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_EXP.c

    r35768 r37066  
    3333#include "pmMoments.h"
    3434#include "pmModelFuncs.h"
     35#include "pmModelClass.h"
    3536#include "pmModel.h"
    3637#include "pmModelUtils.h"
    37 #include "pmModelClass.h"
    3838#include "pmSourceMasks.h"
    3939#include "pmSourceExtendedPars.h"
    4040#include "pmSourceDiffStats.h"
    4141#include "pmSourceSatstar.h"
     42#include "pmSourceLensing.h"
    4243#include "pmSource.h"
    4344#include "pmSourceFitModel.h"
    … …  
    4546#include "pmPSFtry.h"
    4647#include "pmDetections.h"
     48#include "pmModel_CentralPixel.h"
    4749
    4850#include "pmModel_EXP.h"
    … …  
    5456# define PM_MODEL_LIMITS          pmModelLimits_EXP
    5557# define PM_MODEL_RADIUS          pmModelRadius_EXP
     58# define PM_MODEL_SET_FWHM        pmModelSetFWHM_EXP
    5659# define PM_MODEL_FROM_PSF        pmModelFromPSF_EXP
    5760# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_EXP
    … …  
    6265// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    6366// values need to be pixel coords
     67//
     68
     69// Notes on changing kappa value from 1.70056 to 1.678
     70// I'm using a functional form f(x,y) = Io exp(-kappa (r / r_e)). 
     71// The article by Graham & Driver (2005) uses a form Ie exp(-bn [(r / r_e) -1])
     72// which is equal to Ie exp(-bn (r / r_e)) exp(bn). 
     73// Thus, my Io = Ie exp(bn) and my kappa is their bn.
     74// My value of kappa is 1.700, their value for bn is 1.678., so I am off by a small amount there (1.5%). 
     75
     76
     77#define KAPPA_EXP 1.678
     78#define OLD_KAPP_EXP 1.70056
     79
    6480
    6581// Lax parameter limits
    6682static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 };
    67 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     83static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    6884
    6985// Moderate parameter limits
    … …  
    7894static float *paramsMinUse = paramsMinLax;
    7995static float *paramsMaxUse = paramsMaxLax;
    80 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5};
     96static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5};
    8197
    8298static bool limitsApply = true;         // Apply limits?
    8399
    84 # include "pmModel_SERSIC.CP.h"
     100// # include "pmModel_SERSIC.CP.h"
     101
     102// the problems I'm having with the SERSIC-like functions are:
     103// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius)
     104// 2) getting the central pixel right
     105// 3) getting the derivaties right.
    85106
    86107psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    101122    psAssert (z >= 0, "do not allow negative z values in model");
    102123
    103     float index = 1.0;
    104     float par7 = 0.5;
    105     float bn = 1.9992*index - 0.3271;
    106     float Io = exp(bn);
    107 
    108     psF32 f2 = bn*sqrt(z);
    109     psF32 f1 = Io*exp(-f2);
    110 
     124    // for EXP, we can hard-wire kappa(1):
     125    // float index = 1.0;
     126    float kappa = KAPPA_EXP;
     127
     128    // sqrt(z) is r
     129    float q = kappa*sqrt(z);
     130    psF32 f0 = exp(-q);
     131
     132    assert (isfinite(q));
     133
     134    // only worry about the central 4 pixels at most
    111135    psF32 radius = hypot(X, Y);
    112     if (radius < 1.0) {
    113 
    114         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    115 
    116         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    117         psEllipseAxes axes;
    118         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    119 
    120         // get the central pixel flux from the lookup table
    121         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    122         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    123         float yPix = (index - centralPixelYo) / centralPixeldY;
    124         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    125 
    126         // the integral of a Sersic has an analytical form as follows:
    127         float logGamma = lgamma(2.0*index);
    128         float bnFactor = pow(bn, 2.0*index);
    129         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    130 
    131         // XXX interpolate to get the value
    132         // XXX for the moment, just integerize
    133         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    134         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    135        
    136         float px1 = 1.0 / PAR[PM_PAR_SXX];
    137         float py1 = 1.0 / PAR[PM_PAR_SYY];
    138         float z10 = PS_SQR(px1);
    139         float z01 = PS_SQR(py1);
    140 
    141         // which pixels do we need for this interpolation?
    142         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    143         if ((X >= 0) && (Y >= 0)) {
    144             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    145             float V00 = Vcenter;
    146             float V10 = Io*exp(-bn*pow(z10,par7));
    147             float V01 = Io*exp(-bn*pow(z01,par7));
    148             float V11 = Io*exp(-bn*pow(z11,par7));
    149             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    150         }
    151         if ((X < 0) && (Y >= 0)) {
    152             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    153             float V00 = Io*exp(-bn*pow(z10,par7));
    154             float V10 = Vcenter;
    155             float V01 = Io*exp(-bn*pow(z11,par7));
    156             float V11 = Io*exp(-bn*pow(z01,par7));
    157             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    158         }
    159         if ((X >= 0) && (Y < 0)) {
    160             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    161             float V00 = Io*exp(-bn*pow(z01,par7));
    162             float V10 = Io*exp(-bn*pow(z11,par7));
    163             float V01 = Vcenter;
    164             float V11 = Io*exp(-bn*pow(z10,par7));
    165             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    166         }
    167         if ((X < 0) && (Y < 0)) {
    168             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    169             float V00 = Io*exp(-bn*pow(z11,par7));
    170             float V10 = Io*exp(-bn*pow(z10,par7));
    171             float V01 = Io*exp(-bn*pow(z01,par7));
    172             float V11 = Vcenter;
    173             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    174         }
    175     }   
    176 
    177     psF32 z0 = PAR[PM_PAR_I0]*f1;
    178     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    179 
    180     assert (isfinite(f2));
     136    if (radius <= 1.5) {
     137        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
     138    }
     139    assert (isfinite(f0));
     140
     141    psF32 f1 = PAR[PM_PAR_I0]*f0;
     142    psF32 f = PAR[PM_PAR_SKY] + f1;
     143
    181144    assert (isfinite(f1));
    182     assert (isfinite(z0));
    183     assert (isfinite(f0));
     145    assert (isfinite(f));
    184146
    185147    if (deriv != NULL) {
    … …  
    187149
    188150        dPAR[PM_PAR_SKY]  = +1.0;
    189         dPAR[PM_PAR_I0]   = +f1;
    190 
    191         // gradient is infinite for z = 0; saturate at z = 0.01
    192         // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
    193         psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
    194 
    195         // XXX dampen SXX and SYY as in GAUSS?
    196         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    197         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    198         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
    199         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    200         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    201     }
    202     return (f0);
     151        dPAR[PM_PAR_I0]   = +f0;
     152
     153        if (z > 0.01) {
     154          float z1 = 0.5*f1*kappa/sqrt(z);
     155          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     156          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     157          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     158          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     159          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     160        } else {
     161          // gradient -> 0 for z -> 0, but has undef form
     162          float z1 = 0.5*f1*kappa;
     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        }
     169    }
     170    return (f);
    203171}
    204172
    … …  
    284252bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    285253{
     254    // for the moment, we are going to require moments and KronFlux
     255    if (!source->moments) return false;
     256    pmMoments *moments = source->moments;
     257
     258    if (!isfinite(moments->KronFlux)) return false;
     259    if (!isfinite(moments->Mrf)) return false;
     260    if (moments->Mrf < 0.0) return false;
     261
    286262    psF32 *PAR  = model->params->data.F32;
    287263
    … …  
    289265    PAR[PM_PAR_SKY]  = 0.0;
    290266
    291     // set the shape parameters
    292     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
    293       return false;
    294     }
    295 
    296     // set the model normalization
    297     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    298       return false;
    299     }
     267    psEllipseMoments emoments;
     268    emoments.x2 = moments->Mxx;
     269    emoments.xy = moments->Mxy;
     270    emoments.y2 = moments->Myy;
     271
     272    // force the axis ratio to be < 20.0
     273    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
     274
     275    if (!isfinite(axes.major)) return false;
     276    if (!isfinite(axes.minor)) return false;
     277    if (!isfinite(axes.theta)) return false;
     278
     279    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width
     280    float scale = moments->Mrf / axes.major;
     281    axes.major *= scale;
     282    axes.minor *= scale;
     283
     284    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
     285
     286    // psEllipseAxes axes;
     287    // use the code in SetShape here to avoid doing this 2x
     288    // pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     289
     290    // float norm = pmSersicNorm (4);  // hardwire
     291    float norm = 0.34578;
     292    float normFlux = 2.0 * M_PI * axes.major * axes.minor * norm;
     293    PAR[PM_PAR_I0] = moments->KronFlux / normFlux;
    300294
    301295    // set the model position
    … …  
    306300    return(true);
    307301}
    308 
    309302// An exponential model is equivalent to a Sersic with index = 1.0
    310303psF64 PM_MODEL_FLUX (const psVector *params)
    … …  
    314307    psEllipseAxes axes;
    315308    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    316     float AspectRatio = axes.minor / axes.major;
    317 
    318     float index = 1.0;
    319     float bn = 1.9992*index - 0.3271;
    320 
    321     // the integral of a Sersic has an analytical form as follows:
    322     float logGamma = lgamma(2.0*index);
    323     float bnFactor = pow(bn, 2.0*index);
    324     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    325    
    326     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    327 
    328     return(Flux);
     309
     310    // static value for EXP:
     311    float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr
     312
     313    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     314
     315    return(flux);
    329316}
    330317
    … …  
    345332    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    346333
    347     // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
    348     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    349     psF64 radius = axes.major * sqrt (2.0) * zn;
     334    // static value for EXP:
     335    float kappa = KAPPA_EXP;
     336
     337    // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
     338    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     339    psF64 radius = axes.major * zn;
    350340
    351341    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",
    352342              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
    353343    return (radius);
     344}
     345
     346psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
     347  return (NAN);
    354348}
    355349
    … …  
    375369    // the 2D PSF model fits polarization terms (E0,E1,E2)
    376370    // convert to shape terms (SXX,SYY,SXY)
    377     bool useReff = pmModelUseReff (modelPSF->type);
     371    bool useReff = modelPSF->class->useReff;
    378372    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    379373        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    … …  
    429423    // convert to shape terms (SXX,SYY,SXY)
    430424    // XXX user-defined value for limit?
    431     bool useReff = pmModelUseReff (model->type);
     425    bool useReff = model->class->useReff;
    432426    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    433427        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
    … …  
    501495    return;
    502496}
     497
     498# if (0)
     499void bilin_inter_function () {
     500        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     501        psEllipseAxes axes;
     502        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     503
     504        // get the central pixel flux from the lookup table
     505        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     506        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     507        float yPix = (index - centralPixelYo) / centralPixeldY;
     508        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     509
     510        // the integral of a Sersic has an analytical form as follows:
     511        float logGamma = lgamma(2.0*index);
     512        float bnFactor = pow(bn, 2.0*index);
     513        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     514
     515        // XXX interpolate to get the value
     516        // XXX for the moment, just integerize
     517        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     518        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     519       
     520        float px1 = 1.0 / PAR[PM_PAR_SXX];
     521        float py1 = 1.0 / PAR[PM_PAR_SYY];
     522        float z10 = PS_SQR(px1);
     523        float z01 = PS_SQR(py1);
     524
     525        // which pixels do we need for this interpolation?
     526        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     527        if ((X >= 0) && (Y >= 0)) {
     528            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     529            float V00 = Vcenter;
     530            float V10 = Io*exp(-bn*pow(z10,par7));
     531            float V01 = Io*exp(-bn*pow(z01,par7));
     532            float V11 = Io*exp(-bn*pow(z11,par7));
     533            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     534        }
     535        if ((X < 0) && (Y >= 0)) {
     536            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     537            float V00 = Io*exp(-bn*pow(z10,par7));
     538            float V10 = Vcenter;
     539            float V01 = Io*exp(-bn*pow(z11,par7));
     540            float V11 = Io*exp(-bn*pow(z01,par7));
     541            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     542        }
     543        if ((X >= 0) && (Y < 0)) {
     544            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     545            float V00 = Io*exp(-bn*pow(z01,par7));
     546            float V10 = Io*exp(-bn*pow(z11,par7));
     547            float V01 = Vcenter;
     548            float V11 = Io*exp(-bn*pow(z10,par7));
     549            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     550        }
     551        if ((X < 0) && (Y < 0)) {
     552            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     553            float V00 = Io*exp(-bn*pow(z11,par7));
     554            float V10 = Io*exp(-bn*pow(z10,par7));
     555            float V01 = Io*exp(-bn*pow(z01,par7));
     556            float V11 = Vcenter;
     557            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     558        }
     559}
     560# endif
Note: See TracChangeset for help on using the changeset viewer.