IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 3, 2013, 2:37:22 PM (13 years ago)
Author:
eugene
Message:

deprecate KiiOpen,KiiClose (now KapaOpen,etc); major rework of psEllipse translations : use common functions pmModelAxesToParams and pmModelParamsToAxes ; use new convergence method in pmPCM_MinimizeChisq; add convergence crerition options to psMinimization; threaded versions of pmPSFtryFitEXT and pmPSFtryFitPSF

Location:
trunk/psModules
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

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

    r35560 r35768  
    123123
    124124        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    125         psEllipseShape shape;
    126 
    127         shape.sx  = PAR[PM_PAR_SXX];
    128         shape.sy  = PAR[PM_PAR_SYY];
    129         shape.sxy = PAR[PM_PAR_SXY];
    130 
    131         // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
    132         psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     125        psEllipseAxes axes;
     126        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    133127
    134128        // get the central pixel flux from the lookup table
    … …  
    238232        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    239233        // angle and let f2,f1 fight it out
    240         q2 = 0.5*sqrtf(q1);
     234        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     235        q2 = 2.0*0.5*sqrtf(q1);
    241236    }
    242237
    … …  
    303298
    304299    // set the shape parameters
    305     // XXX adjust this?
    306     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     300    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
    307301      return false;
    308302    }
    … …  
    330324psF64 PM_MODEL_FLUX (const psVector *params)
    331325{
    332     psEllipseShape shape;
    333 
    334326    psF32 *PAR = params->data.F32;
    335327
    336     shape.sx  = PAR[PM_PAR_SXX];
    337     shape.sy  = PAR[PM_PAR_SYY];
    338     shape.sxy = PAR[PM_PAR_SXY];
    339 
    340     // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
    341     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     328    psEllipseAxes axes;
     329    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    342330    float AspectRatio = axes.minor / axes.major;
    343331
    … …  
    359347psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    360348{
    361     psEllipseShape shape;
    362 
    363349    psF32 *PAR = params->data.F32;
    364350
    … …  
    370356        return (1.0);
    371357
    372     shape.sx  = PAR[PM_PAR_SXX];
    373     shape.sy  = PAR[PM_PAR_SYY];
    374     shape.sxy = PAR[PM_PAR_SXY];
    375 
    376     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     358    psEllipseAxes axes;
     359    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    377360
    378361    // f = Io exp(-z^n) -> z^n = ln(Io/f)
    … …  
    382365    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",
    383366              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
    384 
    385367    return (radius);
    386368}
    … …  
    407389    // the 2D PSF model fits polarization terms (E0,E1,E2)
    408390    // convert to shape terms (SXX,SYY,SXY)
    409     if (!pmPSF_FitToModel (out, 0.1)) {
     391    bool useReff = pmModelUseReff (modelPSF->type);
     392    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    410393        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    411394        return false;
    … …  
    460443    // convert to shape terms (SXX,SYY,SXY)
    461444    // XXX user-defined value for limit?
    462     if (!pmPSF_FitToModel (PAR, 0.1)) {
     445    bool useReff = pmModelUseReff (model->type);
     446    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    463447        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    464448        return false;

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

    r35560 r35768  
    115115
    116116        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    117         psEllipseShape shape;
    118 
    119         shape.sx  = PAR[PM_PAR_SXX];
    120         shape.sy  = PAR[PM_PAR_SYY];
    121         shape.sxy = PAR[PM_PAR_SXY];
    122 
    123         // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
    124         psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     117        psEllipseAxes axes;
     118        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    125119
    126120        // get the central pixel flux from the lookup table
    … …  
    230224        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    231225        // angle and let f2,f1 fight it out
    232         q2 = 0.5*sqrtf(q1);
     226        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     227        q2 = 2.0*0.5*sqrtf(q1);
    233228    }
    234229
    … …  
    295290
    296291    // set the shape parameters
    297     // XXX adjust this?
    298     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     292    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
    299293      return false;
    300294    }
    … …  
    316310psF64 PM_MODEL_FLUX (const psVector *params)
    317311{
    318     psEllipseShape shape;
    319 
    320312    psF32 *PAR = params->data.F32;
    321313
    322     shape.sx  = PAR[PM_PAR_SXX];
    323     shape.sy  = PAR[PM_PAR_SYY];
    324     shape.sxy = PAR[PM_PAR_SXY];
    325 
    326     // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
    327     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     314    psEllipseAxes axes;
     315    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    328316    float AspectRatio = axes.minor / axes.major;
    329317
    … …  
    345333psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    346334{
    347     psEllipseShape shape;
    348 
    349335    psF32 *PAR = params->data.F32;
    350336
    … …  
    356342        return (1.0);
    357343
    358     shape.sx  = PAR[PM_PAR_SXX];
    359     shape.sy  = PAR[PM_PAR_SYY];
    360     shape.sxy = PAR[PM_PAR_SXY];
    361 
    362     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     344    psEllipseAxes axes;
     345    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    363346
    364347    // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
    … …  
    392375    // the 2D PSF model fits polarization terms (E0,E1,E2)
    393376    // convert to shape terms (SXX,SYY,SXY)
    394     if (!pmPSF_FitToModel (out, 0.1)) {
     377    bool useReff = pmModelUseReff (modelPSF->type);
     378    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    395379        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    396380        return false;
    … …  
    445429    // convert to shape terms (SXX,SYY,SXY)
    446430    // XXX user-defined value for limit?
    447     if (!pmPSF_FitToModel (PAR, 0.1)) {
     431    bool useReff = pmModelUseReff (model->type);
     432    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    448433        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    449434        return false;

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

    r35560 r35768  
    129129    float q2 = NAN;
    130130    if (nParam == PM_PAR_SXY) {
    131         float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    132         float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     131        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     132        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
     133        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
    133134        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    134135        q1 = (q1 < 0.0) ? 0.0 : q1;
    … …  
    200201
    201202    // set the shape parameters
    202     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     203    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
    203204      return false;
    204205    }
    … …  
    219220psF64 PM_MODEL_FLUX (const psVector *params)
    220221{
    221 
    222     psEllipseShape shape;
    223 
    224222    psF32 *PAR = params->data.F32;
    225223
    226     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    227     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    228     shape.sxy = PAR[PM_PAR_SXY];
     224    psEllipseAxes axes;
     225    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    229226
    230227    // Area is equivalent to 2 pi sigma^2
    231     // axes ratio < 20
    232     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    233228    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    234229
    … …  
    242237psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    243238{
    244     psEllipseShape shape;
    245 
    246239    psF32 *PAR = params->data.F32;
    247240
    … …  
    253246        return (1.0);
    254247
    255     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    256     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    257     shape.sxy = PAR[PM_PAR_SXY];
    258 
    259     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     248    psEllipseAxes axes;
     249    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     250
    260251    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
    261252    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",
    … …  
    285276    }
    286277
    287     // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
    288     // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
    289 
    290278    // the 2D PSF model fits polarization terms (E0,E1,E2)
    291279    // convert to shape terms (SXX,SYY,SXY)
    292     // XXX user-defined value for limit?
    293     if (!pmPSF_FitToModel (out, 0.1)) {
    294         // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
     280    bool useReff = pmModelUseReff (modelPSF->type);
     281    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    295282        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    296283        return false;
    … …  
    343330    // the 2D PSF model fits polarization terms (E0,E1,E2)
    344331    // convert to shape terms (SXX,SYY,SXY)
    345     // XXX user-defined value for limit?
    346     if (!pmPSF_FitToModel (PAR, 0.1)) {
     332    bool useReff = pmModelUseReff (model->type);
     333    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    347334        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    348335        return false;

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

    r35560 r35768  
    129129    float q2 = NAN;
    130130    if (nParam == PM_PAR_SXY) {
    131         float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    132         float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     131        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     132        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
     133        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
    133134        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    134135        q1 = (q1 < 0.0) ? 0.0 : q1;
    … …  
    201202
    202203    // set the shape parameters
    203     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     204    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
    204205      return false;
    205206    }
    … …  
    222223{
    223224    float z, norm;
    224     psEllipseShape shape;
    225225
    226226    psF32 *PAR = params->data.F32;
    227227
    228     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    229     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    230     shape.sxy = PAR[PM_PAR_SXY];
    231 
    232     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     228    psEllipseAxes axes;
     229    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     230
    233231    float AspectRatio = axes.minor / axes.major;
    234232
    … …  
    262260{
    263261    psF64 z;
    264     psEllipseShape shape;
    265262
    266263    psF32 *PAR = params->data.F32;
    … …  
    273270        return (1.0);
    274271
    275     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    276     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    277     shape.sxy = PAR[PM_PAR_SXY];
    278 
    279     // this estimates the radius assuming f(z) is roughly exp(-z)
    280     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     272    psEllipseAxes axes;
     273    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     274
    281275    psF64 sigma = axes.major;
    282276
    … …  
    347341    }
    348342
    349     // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
    350     // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
    351 
    352343    // the 2D PSF model fits polarization terms (E0,E1,E2)
    353344    // convert to shape terms (SXX,SYY,SXY)
    354     if (!pmPSF_FitToModel (out, 0.1)) {
     345    bool useReff = pmModelUseReff (modelPSF->type);
     346    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    355347        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    356348        return false;
    … …  
    403395    // the 2D PSF model fits polarization terms (E0,E1,E2)
    404396    // convert to shape terms (SXX,SYY,SXY)
    405     // XXX user-defined value for limit?
    406     if (!pmPSF_FitToModel (PAR, 0.1)) {
     397    bool useReff = pmModelUseReff (model->type);
     398    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    407399        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    408400        return false;

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

    r35560 r35768  
    11/******************************************************************************
    2  * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
    3  * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
    4  * models use a psVector to represent the set of parameters, with the sequence used to specify
    5  * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
    6  * specifics of the model.  All models which are used as a PSF representations share a few
    7  * parameters, for which # define names are listed in pmModel.h:
     2 * this file defines the PS1_V1 source shape model.  Note that these model functions are
     3 * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
     4 * their own.  The models use a psVector to represent the set of parameters, with the
     5 * sequence used to specify the meaning of the parameter.  The meaning of the parameters
     6 * may thus vary depending on the specifics of the model.  All models which are used as a
     7 * PSF representations share a few parameters, for which # define names are listed in
     8 * pmModel.h:
    89
    910   power-law with fitted linear term
    … …  
    148149    float q2 = NAN;
    149150    if (nParam == PM_PAR_SXY) {
    150         float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    151         float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     151        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
     152        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
    152153        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    153154        q1 = (q1 < 0.0) ? 0.0 : q1;
    … …  
    220221
    221222    // set the shape parameters
    222     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     223    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
    223224      return false;
    224225    }
    … …  
    244245{
    245246    float z, norm;
    246     psEllipseShape shape;
    247247
    248248    psF32 *PAR = params->data.F32;
    249249
    250     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    251     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    252     shape.sxy = PAR[PM_PAR_SXY];
    253 
    254     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     250    psEllipseAxes axes;
     251    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    255252    float AspectRatio = axes.minor / axes.major;
    256253
    … …  
    284281{
    285282    psF64 z;
    286     psEllipseShape shape;
    287283
    288284    psF32 *PAR = params->data.F32;
    … …  
    292288    if (flux >= PAR[PM_PAR_I0]) return 1.0;
    293289
    294     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    295     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    296     shape.sxy = PAR[PM_PAR_SXY];
    297 
    298     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     290    psEllipseAxes axes;
     291    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    299292    psF64 sigma = axes.major;
    300293
    … …  
    363356    // the 2D PSF model fits polarization terms (E0,E1,E2)
    364357    // convert to shape terms (SXX,SYY,SXY)
    365     if (!pmPSF_FitToModel (out, 0.1)) {
     358    bool useReff = pmModelUseReff (modelPSF->type);
     359    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    366360        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    367361        return false;
    … …  
    416410    // convert to shape terms (SXX,SYY,SXY)
    417411    // XXX user-defined value for limit?
    418     if (!pmPSF_FitToModel (PAR, 0.1)) {
     412    bool useReff = pmModelUseReff (model->type);
     413    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    419414        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    420415        return false;

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

    r35560 r35768  
    11/******************************************************************************
    2  * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
    3  * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
    4  * statements of their own.  The models use a psVector to represent the set of parameters, with
    5  * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
    6  * may thus vary depending on the specifics of the model.  All models which are used a PSF
    7  * representations share a few parameters, for which # define names are listed in pmModel.h:
     2 * this file defines the QGAUSS source shape model.  Note that these model functions are
     3 * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
     4 * their own.  The models use a psVector to represent the set of parameters, with the
     5 * sequence used to specify the meaning of the parameter.  The meaning of the parameters
     6 * may thus vary depending on the specifics of the model.  All models which are used as a
     7 * PSF representations share a few parameters, for which # define names are listed in
     8 * pmModel.h:
    89
    910   power-law with fitted linear term
    … …  
    1415   * PM_PAR_XPOS 2  - X center of object
    1516   * PM_PAR_YPOS 3  - Y center of object
    16    * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
    17    * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     17   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
     18   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    1819   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    1920   * PM_PAR_7   7   - amplitude of the linear component (k)
    … …  
    138139# define AR_MAX 20.0
    139140# define AR_RATIO 0.99
    140 
    141141bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
    142142{
    … …  
    149149    float q2 = NAN;
    150150    if (nParam == PM_PAR_SXY) {
    151         float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    152         float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     151        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     152        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
     153        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
    153154        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    154155        q1 = (q1 < 0.0) ? 0.0 : q1;
    … …  
    203204          return true;
    204205      }
    205     default:
     206      default:
    206207        psAbort("invalid choice for limits");
    207208    }
    … …  
    221222
    222223    // set the shape parameters
    223     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     224    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
    224225      return false;
    225226    }
    … …  
    245246{
    246247    float z, norm;
    247     psEllipseShape shape;
    248248
    249249    psF32 *PAR = params->data.F32;
    250250
    251     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    252     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    253     shape.sxy = PAR[PM_PAR_SXY];
    254 
    255     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     251    psEllipseAxes axes;
     252    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    256253    float AspectRatio = axes.minor / axes.major;
    257254
    … …  
    285282{
    286283    psF64 z;
    287     psEllipseShape shape;
    288284
    289285    psF32 *PAR = params->data.F32;
    … …  
    293289    if (flux >= PAR[PM_PAR_I0]) return 1.0;
    294290
    295     // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    296 
    297     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    298     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    299     shape.sxy = PAR[PM_PAR_SXY];
    300 
    301     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     291    psEllipseAxes axes;
     292    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    302293    psF64 sigma = axes.major;
    303294
    … …  
    307298        return ( sigma * sqrt (2.0 * z) );
    308299    }
    309 
    310300    psF64 limit = flux / PAR[PM_PAR_I0];
    311301
    … …  
    367357    // the 2D PSF model fits polarization terms (E0,E1,E2)
    368358    // convert to shape terms (SXX,SYY,SXY)
    369     if (!pmPSF_FitToModel (out, 0.1)) {
     359    bool useReff = pmModelUseReff (modelPSF->type);
     360    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    370361        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    371362        return false;
    … …  
    424415    // the 2D PSF model fits polarization terms (E0,E1,E2)
    425416    // convert to shape terms (SXX,SYY,SXY)
    426     // XXX user-defined value for limit?
    427     if (!pmPSF_FitToModel (PAR, 0.1)) {
     417    bool useReff = pmModelUseReff (model->type);
     418    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    428419        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    429420        return false;

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

    r35560 r35768  
    138138    float q2 = NAN;
    139139    if (nParam == PM_PAR_SXY) {
    140         float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
    141         float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
     140        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     141        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
     142        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
    142143        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
    143144        q1 = (q1 < 0.0) ? 0.0 : q1;
    … …  
    210211
    211212    // set the shape parameters
    212     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
     213    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
    213214      return false;
    214215    }
    … …  
    234235{
    235236    float z, norm;
    236     psEllipseShape shape;
    237237
    238238    psF32 *PAR = params->data.F32;
    239239
    240     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    241     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    242     shape.sxy = PAR[PM_PAR_SXY];
    243 
    244     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     240    psEllipseAxes axes;
     241    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    245242    float AspectRatio = axes.minor / axes.major;
    246243
    … …  
    274271{
    275272    psF64 z;
    276     psEllipseShape shape;
    277273
    278274    psF32 *PAR = params->data.F32;
    … …  
    285281        return (1.0);
    286282
    287     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    288     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
    289     shape.sxy = PAR[PM_PAR_SXY];
    290 
    291     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     283    psEllipseAxes axes;
     284    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
    292285    psF64 sigma = axes.major;
    293286
    … …  
    357350    // the 2D PSF model fits polarization terms (E0,E1,E2)
    358351    // convert to shape terms (SXX,SYY,SXY)
    359     if (!pmPSF_FitToModel (out, 0.1)) {
     352    bool useReff = pmModelUseReff (modelPSF->type);
     353    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    360354        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    361355        return false;
    … …  
    409403    // the 2D PSF model fits polarization terms (E0,E1,E2)
    410404    // convert to shape terms (SXX,SYY,SXY)
    411     // XXX user-defined value for limit?
    412     if (!pmPSF_FitToModel (PAR, 0.1)) {
     405    bool useReff = pmModelUseReff (model->type);
     406    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    413407        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    414408        return false;

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

    r35560 r35768  
    125125
    126126        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    127         psEllipseShape shape;
    128 
    129         shape.sx  = PAR[PM_PAR_SXX];
    130         shape.sy  = PAR[PM_PAR_SYY];
    131         shape.sxy = PAR[PM_PAR_SXY];
    132 
    133         // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
    134         psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     127        psEllipseAxes axes;
     128        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    135129
    136130        // get the central pixel flux from the lookup table
    … …  
    250244        // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
    251245        // angle and let f2,f1 fight it out
    252         q2 = 0.5*sqrtf(q1);
     246        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
     247        q2 = 2.0*0.5*sqrtf(q1);
    253248    }
    254249
    … …  
    347342    axes.major = Rmajor;
    348343    axes.minor = Rminor;
    349     psEllipseShape shape = psEllipseAxesToShape (axes);
    350 
    351     if (!isfinite(shape.sx))  return false;
    352     if (!isfinite(shape.sy))  return false;
    353     if (!isfinite(shape.sxy)) return false;
     344
     345    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
    354346
    355347    float bn = 1.9992*index - 0.3271;
    … …  
    357349    float Io = exp(0.5*bn);
    358350
    359     float Sxx = PS_MAX(0.5, shape.sx);
    360     float Syy = PS_MAX(0.5, shape.sy);
    361 
    362     PAR[PM_PAR_SXX]  = Sxx;
    363     PAR[PM_PAR_SYY]  = Syy;
    364     PAR[PM_PAR_SXY]  = shape.sxy;
    365 
    366351    // set the model normalization (adjust for Sersic best guess)
    367352    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    … …  
    381366psF64 PM_MODEL_FLUX (const psVector *params)
    382367{
    383     psEllipseShape shape;
    384 
    385368    psF32 *PAR = params->data.F32;
    386369
    387     shape.sx  = PAR[PM_PAR_SXX];
    388     shape.sy  = PAR[PM_PAR_SYY];
    389     shape.sxy = PAR[PM_PAR_SXY];
    390 
    391     // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
    392     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     370    psEllipseAxes axes;
     371    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    393372    float AspectRatio = axes.minor / axes.major;
    394373
    … …  
    410389psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    411390{
    412     psEllipseShape shape;
    413 
    414391    psF32 *PAR = params->data.F32;
    415392
    … …  
    421398        return (1.0);
    422399
    423     shape.sx  = PAR[PM_PAR_SXX];
    424     shape.sy  = PAR[PM_PAR_SYY];
    425     shape.sxy = PAR[PM_PAR_SXY];
    426 
    427     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     400    psEllipseAxes axes;
     401    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    428402
    429403    // f = Io exp(-z^n) -> z^n = ln(Io/f)
    … …  
    457431    // the 2D PSF model fits polarization terms (E0,E1,E2)
    458432    // convert to shape terms (SXX,SYY,SXY)
    459     if (!pmPSF_FitToModel (out, 0.1)) {
     433    bool useReff = pmModelUseReff (modelPSF->type);
     434    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    460435        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    461436        return false;
    … …  
    510485    // convert to shape terms (SXX,SYY,SXY)
    511486    // XXX user-defined value for limit?
    512     if (!pmPSF_FitToModel (PAR, 0.1)) {
     487    bool useReff = pmModelUseReff (model->type);
     488    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    513489        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    514490        return false;

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

    r35577 r35768  
    350350    PAR[PM_PAR_SKY]  = 0.0;
    351351
    352     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
    353     psEllipseShape psfShape;
    354     psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
    355     psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
    356     psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
    357     psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
     352    psF32 *psfPAR  = source->modelPSF->params->data.F32;
     353    bool useReff = pmModelUseReff (source->modelPSF->type);
     354
     355    psEllipseAxes psfAxes;
     356    pmModelParamsToAxes (&psfAxes, psfPAR[PM_PAR_SXX], psfPAR[PM_PAR_SXY], psfPAR[PM_PAR_SYY], useReff);
    358357
    359358    psEllipseMoments emoments;
    … …  
    369368    if (!isfinite(axes.theta)) return false;
    370369
    371     float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
    372     //    if (size > psfAxes.major) { size -= psfAxes.major; }
    373     //else { size = psfAxes.major; }
     370    float size = NAN;
     371    if (!isfinite(source->moments->Mrf)) {
     372      size = axes.major;
     373    } else {
     374      size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
     375    }
    374376
    375377    float theta, peak;
Note: See TracChangeset for help on using the changeset viewer.