IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:02:53 PM (15 years ago)
Author:
eugene
Message:
  • add concept of saddlePoints to peaks (not actually used in the end)
  • add tmp flags to mark sources for analysis or not in psphotStack
  • autocode the pmSourceIO_CMF_PS1_* functions
  • use 1D gauss approx for convolution in PCM fitting
  • added pmSourceExtFitPars (not actually used)
  • in model guess, use 1st radial moments to define size (if it exists)
  • include PSF_INST_MAG, AP_MAG, KRON_MAG in xfit output
  • fix the position for extended source fits (avoid instability)
  • Sersic-like models (incl. Exp and Dev) use Reff, not sigma; conversion tools need to respect this
  • only use a single pass on the centroid (unwindowed, but limited to 1.5 sigma radius) - this avoids moving the centroid because of nearby neighbors
  • use symmetrical averaging (geometric mean) to calculated 1st radial moment (and avoid neighbor biases), do not use symm. averaging for the flux
  • fix the integration of the sersic, pgauss, and related model functions.
  • fix the central pixel to have the full flux for sersic-like models (interpolated value)
Location:
trunk/psModules/src/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        55*.la
        66*.lo
         7pmSourceIO_CMF_PS1_V1.c
         8pmSourceIO_CMF_PS1_V2.c
         9pmSourceIO_CMF_PS1_V3.c

  • trunk/psModules/src/objects/models/pmModel_EXP.c

    r31451 r32347  
    8181static bool limitsApply = true;         // Apply limits?
    8282
     83# include "pmModel_SERSIC.CP.h"
     84
    8385psF32 PM_MODEL_FUNC (psVector *deriv,
    8486                     const psVector *params,
    … …  
    9395    psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    9496
    95     // XXX if the elliptical contour is defined in valid way, this step should not be required.
    96     // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
    97     // for negative values of z
    98     // XXX use an assert here to force the elliptical parameters to be correctly determined
    99     // if (z < 0) z = 0;
    100     assert (z >= 0);
    101 
    102     psF32 f2 = sqrt(z);
    103     psF32 f1 = exp(-f2);
     97    // If the elliptical contour is defined in a valid way, we should never trigger this
     98    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
     99    // NaN values for negative values of z
     100    psAssert (z >= 0, "do not allow negative z values in model");
     101
     102    float index = 1.0;
     103    float par7 = 0.5;
     104    float bn = 1.9992*index - 0.3271;
     105    float Io = exp(bn);
     106
     107    psF32 f2 = bn*sqrt(z);
     108    psF32 f1 = Io*exp(-f2);
     109
     110    psF32 radius = hypot(X, Y);
     111    if (radius < 1.0) {
     112
     113        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
     114
     115        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     116        psEllipseShape shape;
     117
     118        shape.sx  = PAR[PM_PAR_SXX];
     119        shape.sy  = PAR[PM_PAR_SYY];
     120        shape.sxy = PAR[PM_PAR_SXY];
     121
     122        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     123        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     124
     125        // get the central pixel flux from the lookup table
     126        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     127        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     128        float yPix = (index - centralPixelYo) / centralPixeldY;
     129        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     130
     131        // the integral of a Sersic has an analytical form as follows:
     132        float logGamma = lgamma(2.0*index);
     133        float bnFactor = pow(bn, 2.0*index);
     134        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     135
     136        // XXX interpolate to get the value
     137        // XXX for the moment, just integerize
     138        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     139        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     140       
     141        float px1 = 1.0 / PAR[PM_PAR_SXX];
     142        float py1 = 1.0 / PAR[PM_PAR_SYY];
     143        float z10 = PS_SQR(px1);
     144        float z01 = PS_SQR(py1);
     145
     146        // which pixels do we need for this interpolation?
     147        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     148        if ((X >= 0) && (Y >= 0)) {
     149            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     150            float V00 = Vcenter;
     151            float V10 = Io*exp(-bn*pow(z10,par7));
     152            float V01 = Io*exp(-bn*pow(z01,par7));
     153            float V11 = Io*exp(-bn*pow(z11,par7));
     154            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     155        }
     156        if ((X < 0) && (Y >= 0)) {
     157            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     158            float V00 = Io*exp(-bn*pow(z10,par7));
     159            float V10 = Vcenter;
     160            float V01 = Io*exp(-bn*pow(z11,par7));
     161            float V11 = Io*exp(-bn*pow(z01,par7));
     162            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     163        }
     164        if ((X >= 0) && (Y < 0)) {
     165            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     166            float V00 = Io*exp(-bn*pow(z01,par7));
     167            float V10 = Io*exp(-bn*pow(z11,par7));
     168            float V01 = Vcenter;
     169            float V11 = Io*exp(-bn*pow(z10,par7));
     170            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     171        }
     172        if ((X < 0) && (Y < 0)) {
     173            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     174            float V00 = Io*exp(-bn*pow(z11,par7));
     175            float V10 = Io*exp(-bn*pow(z10,par7));
     176            float V01 = Io*exp(-bn*pow(z01,par7));
     177            float V11 = Vcenter;
     178            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     179        }
     180    }   
     181
    104182    psF32 z0 = PAR[PM_PAR_I0]*f1;
    105183    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    118196        // gradient is infinite for z = 0; saturate at z = 0.01
    119197        // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
    120         psF32 z1 = (z < 0.01) ? 0.5*z0/sqrt(0.01) : 0.5*z0/sqrt(z);
     198        psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
    121199
    122200        // XXX dampen SXX and SYY as in GAUSS?
    … …  
    216294
    217295    // set the shape parameters
     296    // XXX adjust this?
    218297    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
    219298      return false;
    … …  
    233312}
    234313
     314// An exponential model is equivalent to a Sersic with index = 1.0
    235315psF64 PM_MODEL_FLUX (const psVector *params)
    236316{
    237     float z, norm;
    238317    psEllipseShape shape;
    239318
    240319    psF32 *PAR = params->data.F32;
    241320
    242     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    243     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     321    shape.sx  = PAR[PM_PAR_SXX];
     322    shape.sy  = PAR[PM_PAR_SYY];
    244323    shape.sxy = PAR[PM_PAR_SXY];
    245324
    246     // Area is equivalent to 2 pi sigma^2
     325    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
    247326    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    248     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    249 
    250     // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity]
    251     norm = 0.0;
    252 
    253     # define DZ 0.25
    254 
    255     float f0 = 1.0;
    256     float f1, f2;
    257     for (z = DZ; z < 150; z += DZ) {
    258         f1 = exp(-sqrt(z));
    259         z += DZ;
    260         f2 = exp(-sqrt(z));
    261         norm += f0 + 4*f1 + f2;
    262         f0 = f2;
    263     }
    264     norm *= DZ / 3.0;
    265 
    266     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     327    float AspectRatio = axes.minor / axes.major;
     328
     329    float index = 1.0;
     330    float bn = 1.9992*index - 0.3271;
     331
     332    // the integral of a Sersic has an analytical form as follows:
     333    float logGamma = lgamma(2.0*index);
     334    float bnFactor = pow(bn, 2.0*index);
     335    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     336   
     337    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    267338
    268339    return(Flux);
    … …  
    284355        return (1.0);
    285356
    286     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    287     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     357    shape.sx  = PAR[PM_PAR_SXX];
     358    shape.sy  = PAR[PM_PAR_SYY];
    288359    shape.sxy = PAR[PM_PAR_SXY];
    289360
Note: See TracChangeset for help on using the changeset viewer.