IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects/models/pmModel_SERSIC.c

    r35768 r37066  
    2020   * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
    2121     construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
    22      Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
     22     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using:
    2323     shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
    2424     psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
    … …  
    4343#include "pmMoments.h"
    4444#include "pmModelFuncs.h"
     45#include "pmModelClass.h"
    4546#include "pmModel.h"
    4647#include "pmModelUtils.h"
    47 #include "pmModelClass.h"
    4848#include "pmSourceMasks.h"
    4949#include "pmSourceExtendedPars.h"
    5050#include "pmSourceDiffStats.h"
    5151#include "pmSourceSatstar.h"
     52#include "pmSourceLensing.h"
    5253#include "pmSource.h"
    5354#include "pmSourceFitModel.h"
    … …  
    5556#include "pmPSFtry.h"
    5657#include "pmDetections.h"
     58#include "pmModel_CentralPixel.h"
    5759
    5860#include "pmModel_SERSIC.h"
    … …  
    6466# define PM_MODEL_LIMITS          pmModelLimits_SERSIC
    6567# define PM_MODEL_RADIUS          pmModelRadius_SERSIC
     68# define PM_MODEL_SET_FWHM        pmModelSetFWHM_SERSIC
    6669# define PM_MODEL_FROM_PSF        pmModelFromPSF_SERSIC
    6770# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC
    … …  
    7477
    7578// Lax parameter limits
    76 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
    77 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
     79static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.0625 };
     80static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 };
    7881
    7982// Moderate parameter limits
    … …  
    8891static float *paramsMinUse = paramsMinLax;
    8992static float *paramsMaxUse = paramsMaxLax;
    90 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     93static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0};
    9194
    9295static bool limitsApply = true;         // Apply limits?
    9396
    94 # include "pmModel_SERSIC.CP.h"
     97// # include "pmModel_SERSIC.CP.h"
    9598
    9699psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    111114    psAssert (z >= 0, "do not allow negative z values in model");
    112115
    113     float index = 0.5 / PAR[PM_PAR_7];
    114     float par7 = PAR[PM_PAR_7];
    115     float bn = 1.9992*index - 0.3271;
    116     float Io = exp(bn);
    117 
    118     psF32 f2 = bn*pow(z,par7);
    119     psF32 f1 = Io*exp(-f2);
     116    float Sindex = 0.5 / PAR[PM_PAR_7];
     117    float kappa = pmSersicKappa (Sindex);
     118
     119    float q = kappa*pow(z,PAR[PM_PAR_7]);
     120    psF32 f0 = exp(-q);
     121
     122    assert (isfinite(q));
    120123
    121124    psF32 radius = hypot(X, Y);
    122     if (radius < 1.0) {
    123 
    124         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    125 
    126         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    127         psEllipseAxes axes;
    128         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    129 
    130         // get the central pixel flux from the lookup table
    131         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    132         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    133         float yPix = (index - centralPixelYo) / centralPixeldY;
    134         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    135 
    136         // the integral of a Sersic has an analytical form as follows:
    137         float logGamma = lgamma(2.0*index);
    138         float bnFactor = pow(bn, 2.0*index);
    139         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    140 
    141         // XXX interpolate to get the value
    142         // XXX for the moment, just integerize
    143         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    144         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    145        
    146         float px1 = 1.0 / PAR[PM_PAR_SXX];
    147         float py1 = 1.0 / PAR[PM_PAR_SYY];
    148         float z10 = PS_SQR(px1);
    149         float z01 = PS_SQR(py1);
    150 
    151         // which pixels do we need for this interpolation?
    152         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    153         if ((X >= 0) && (Y >= 0)) {
    154             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    155             float V00 = Vcenter;
    156             float V10 = Io*exp(-bn*pow(z10,par7));
    157             float V01 = Io*exp(-bn*pow(z01,par7));
    158             float V11 = Io*exp(-bn*pow(z11,par7));
    159             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    160         }
    161         if ((X < 0) && (Y >= 0)) {
    162             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    163             float V00 = Io*exp(-bn*pow(z10,par7));
    164             float V10 = Vcenter;
    165             float V01 = Io*exp(-bn*pow(z11,par7));
    166             float V11 = Io*exp(-bn*pow(z01,par7));
    167             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    168         }
    169         if ((X >= 0) && (Y < 0)) {
    170             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    171             float V00 = Io*exp(-bn*pow(z01,par7));
    172             float V10 = Io*exp(-bn*pow(z11,par7));
    173             float V01 = Vcenter;
    174             float V11 = Io*exp(-bn*pow(z10,par7));
    175             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    176         }
    177         if ((X < 0) && (Y < 0)) {
    178             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    179             float V00 = Io*exp(-bn*pow(z11,par7));
    180             float V10 = Io*exp(-bn*pow(z10,par7));
    181             float V01 = Io*exp(-bn*pow(z01,par7));
    182             float V11 = Vcenter;
    183             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    184         }
    185     }   
    186 
    187     psF32 z0 = PAR[PM_PAR_I0]*f1;
    188     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    189 
    190     if (!isfinite(z0)) {
    191         fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
     125    if (radius <= 1.5) {
     126        // Nsub ~ 10*index^2 + 1
     127        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, true); // SERSIC model uses Reff
     128        int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1;
     129        Nsub = PS_MIN (Nsub, 121);
     130        Nsub = PS_MAX (Nsub, 11);
     131        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub);
     132    }
     133    if (!isfinite(f0)) {
     134        fprintf(stderr, "f0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, q);
    192135        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
    193136            PAR[5], PAR[6], PAR[7]);
    194137    }
    195 
    196     assert (isfinite(f2));
     138    assert (isfinite(f0));
     139
     140    psF32 f1 = PAR[PM_PAR_I0]*f0;
     141    psF32 f = PAR[PM_PAR_SKY] + f1;
     142
    197143    assert (isfinite(f1));
    198     assert (isfinite(z0));
    199     assert (isfinite(f0));
     144    assert (isfinite(f));
    200145
    201146    if (deriv != NULL) {
    … …  
    203148
    204149        dPAR[PM_PAR_SKY]  = +1.0;
    205         dPAR[PM_PAR_I0]   = +f1;
    206 
    207         // gradient is infinite for z = 0; saturate at z = 0.01
    208         psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
    209 
    210         dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
    211         dPAR[PM_PAR_7]   *= 3.0;
    212 
    213         assert (isfinite(z1));
     150        dPAR[PM_PAR_I0]   = +f0;
     151
     152        if (z > 0.01) {
     153          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0);
     154          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     155          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     156          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     157          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     158          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     159          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z);
     160        } else {
     161          // gradient -> 0 for z -> 0, but has undef form
     162          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]);
     163          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     164          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     165          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     166          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     167          dPAR[PM_PAR_SXY]  = -1.0*z1;
     168          // dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001);
     169          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001); // factor of 16 to reduce the gain
     170        }
    214171        assert (isfinite(dPAR[PM_PAR_7]));
    215 
    216         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    217         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    218         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    219         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    220         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    221     }
    222     return (f0);
     172    }
     173    return (f);
    223174}
    224175
    … …  
    370321    psEllipseAxes axes;
    371322    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    372     float AspectRatio = axes.minor / axes.major;
    373 
    374     float index = 0.5 / PAR[PM_PAR_7];
    375     float bn = 1.9992*index - 0.3271;
    376 
    377     // the integral of a Sersic has an analytical form as follows:
    378     float logGamma = lgamma(2.0*index);
    379     float bnFactor = pow(bn, 2.0*index);
    380     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    381    
    382     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    383 
    384     return(Flux);
     323
     324    float Sindex = 0.5 / PAR[PM_PAR_7];
     325    float norm = pmSersicNorm (Sindex);
     326
     327    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     328
     329    return(flux);
    385330}
    386331
    … …  
    401346    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    402347
     348    float Sindex = 0.5 / PAR[PM_PAR_7];
     349    float kappa = pmSersicKappa (Sindex);
     350
    403351    // f = Io exp(-z^n) -> z^n = ln(Io/f)
    404     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    405     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     352    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     353    psF64 radius = axes.major * pow(zn, Sindex);
    406354
    407355    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f",
    408356              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]);
    409357    return (radius);
     358}
     359
     360psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
     361  return (NAN);
    410362}
    411363
    … …  
    431383    // the 2D PSF model fits polarization terms (E0,E1,E2)
    432384    // convert to shape terms (SXX,SYY,SXY)
    433     bool useReff = pmModelUseReff (modelPSF->type);
     385    bool useReff = modelPSF->class->useReff;
    434386    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    435387        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    … …  
    485437    // convert to shape terms (SXX,SYY,SXY)
    486438    // XXX user-defined value for limit?
    487     bool useReff = pmModelUseReff (model->type);
     439    bool useReff = model->class->useReff;
    488440    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    489441        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
Note: See TracChangeset for help on using the changeset viewer.