IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27266


Ignore:
Timestamp:
Mar 12, 2010, 9:35:51 AM (16 years ago)
Author:
eugene
Message:

replace bogus radius search code with real Newton-Raphson method; make psf model limits consistent with PS1_V1

Location:
branches/eam_branches/20100225/psModules/src/objects/models
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20100225/psModules/src/objects/models/pmModel_PGAUSS.c

    r26916 r27266  
    236236psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    237237{
    238     psF64 z, f;
     238    psF64 z;
    239239    int Nstep = 0;
    240240    psEllipseShape shape;
    … …  
    271271    z0 = 0.0;
    272272
    273     // perform a type of bisection to find the value
    274     float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
    275     float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
    276     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    277         z = 0.5*(z0 + z1);
    278         f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
    279         if (f > limit) {
    280             z0 = z;
    281             f0 = f;
    282         } else {
    283             z1 = z;
    284             f1 = f;
    285         }
    286         Nstep ++;
    287     }
     273    // starting guess:
     274    z = 0.5*(z0 + z1);
     275    float dz = 1.0;
     276
     277    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     278        // use Newton-Raphson to minimize f(z) - limit = 0
     279        float dqdz = (1.0 + z + z*z/2.0);
     280        float q = (dqdz + z*z*z/6.0);
     281
     282        float f = 1.0 / q;
     283        float dfdz = -dqdz * f / q;
     284
     285        dz = (f - limit) / dfdz;
     286
     287        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     288        z -= dz;
     289    }
     290
    288291    psF64 radius = sigma * sqrt (2.0 * z);
    289292

  • branches/eam_branches/20100225/psModules/src/objects/models/pmModel_PS1_V1.c

    r26916 r27266  
    267267psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    268268{
    269     psF64 z, f;
     269    psF64 z;
    270270    int Nstep = 0;
    271271    psEllipseShape shape;
    … …  
    273273    psF32 *PAR = params->data.F32;
    274274
    275     if (flux <= 0) {
    276         return 1.0;
    277     }
    278     if (PAR[PM_PAR_I0] <= 0) {
    279         return 1.0;
    280     }
    281     if (flux >= PAR[PM_PAR_I0]) {
    282         return 1.0;
    283     }
    284     if (PAR[PM_PAR_7] == 0.0) {
    285         return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    286     }
     275    if (flux <= 0) return 1.0;
     276    if (PAR[PM_PAR_I0] <= 0) return 1.0;
     277    if (flux >= PAR[PM_PAR_I0]) return 1.0;
     278    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    287279
    288280    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    305297    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
    306298
    307     // perform a type of bisection to find the value
    308     float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, ALPHA));
    309     float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, ALPHA));
    310     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    311         z = 0.5*(z0 + z1);
    312         f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    313         if (f > limit) {
    314             z0 = z;
    315             f0 = f;
    316         } else {
    317             z1 = z;
    318             f1 = f;
    319         }
    320         Nstep ++;
     299    // starting guess:
     300    z = 0.5*(z0 + z1);
     301    float dz = 1.0;
     302
     303    // use Newton-Raphson to minimize f(z) - limit = 0
     304    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     305        float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     306        float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
     307
     308        float f = 1.0 / q;
     309        float dfdz = -dqdz * f / q;
     310
     311        dz = (f - limit) / dfdz;
     312
     313        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     314        z -= dz;
    321315    }
    322316    psF64 radius = sigma * sqrt (2.0 * z);
    … …  
    350344    // convert to shape terms (SXX,SYY,SXY)
    351345    if (!pmPSF_FitToModel (out, 0.1)) {
    352         // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    353346        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    354347        return false;
    … …  
    475468}
    476469
    477 
    478470# undef PM_MODEL_FUNC
    479471# undef PM_MODEL_FLUX

  • branches/eam_branches/20100225/psModules/src/objects/models/pmModel_QGAUSS.c

    r26916 r27266  
    3939# define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
    4040
     41# define ALPHA   2.250
     42# define ALPHA_M 1.250
     43
    4144// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    4245// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
    … …  
    4447
    4548// Lax parameter limits
    46 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
     49static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
    4750static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    4851
    4952// Moderate parameter limits
    50 static float *paramsMinModerate = paramsMinLax;
    51 static float *paramsMaxModerate = paramsMaxLax;
     53// Tolerate a small divot (k < 0)
     54static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
     55static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    5256
    5357// Strict parameter limits
    54 static float *paramsMinStrict = paramsMinLax;
    55 static float *paramsMaxStrict = paramsMaxLax;
     58// k = PAR_7 < 0 is very undesirable (big divot in the middle)
     59static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
     60static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
    5661
    5762// Parameter limits to use
    5863static float *paramsMinUse = paramsMinLax;
    5964static float *paramsMaxUse = paramsMaxLax;
    60 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     65static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
    6166
    6267static bool limitsApply = true;         // Apply limits?
    … …  
    8186    assert (z >= 0);
    8287
    83     psF32 zp = pow(z,1.25);
     88    psF32 zp = pow(z,ALPHA_M);
    8489    psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
    8590
    … …  
    9297        // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
    9398        psF32 t = r1*r;
    94         psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
     99        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
    95100
    96101        dPAR[PM_PAR_SKY]  = +1.0;
    … …  
    246251    float f1, f2;
    247252    for (z = DZ; z < 50; z += DZ) {
    248         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     253        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    249254        z += DZ;
    250         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     255        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    251256        norm += f0 + 4*f1 + f2;
    252257        f0 = f2;
    … …  
    263268psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    264269{
    265     psF64 z, f;
     270    psF64 z;
    266271    int Nstep = 0;
    267272    psEllipseShape shape;
    … …  
    269274    psF32 *PAR = params->data.F32;
    270275
    271     if (flux <= 0)
    272         return (1.0);
    273     if (PAR[PM_PAR_I0] <= 0)
    274         return (1.0);
    275     if (flux >= PAR[PM_PAR_I0])
    276         return (1.0);
     276    if (flux <= 0) return 1.0;
     277    if (PAR[PM_PAR_I0] <= 0) return 1.0;
     278    if (flux >= PAR[PM_PAR_I0]) return 1.0;
     279    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
    277280
    278281    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    … …  
    290293
    291294    // choose a z value guaranteed to be beyond our limit
    292     float z0 = pow((1.0 / limit), (1.0 / 2.25));
    293     psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
    294     float z1 = (1.0 / limit) / PAR[PM_PAR_7];
     295    float z0 = 0.0;
     296    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
    295297    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
    296     z1 = PS_MAX (z0, z1);
    297     z0 = 0.0;
    298 
    299     // perform a type of bisection to find the value
    300     float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
    301     float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
    302     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    303         z = 0.5*(z0 + z1);
    304         f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    305         if (f > limit) {
    306             z0 = z;
    307             f0 = f;
    308         } else {
    309             z1 = z;
    310             f1 = f;
    311         }
    312         Nstep ++;
     298    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
     299
     300    // starting guess:
     301    z = 0.5*(z0 + z1);
     302    float dz = 1.0;
     303
     304    // use Newton-Raphson to minimize f(z) - limit = 0
     305    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     306        float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     307        float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
     308
     309        float f = 1.0 / q;
     310        float dfdz = -dqdz * f / q;
     311
     312        dz = (f - limit) / dfdz;
     313
     314        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     315        z -= dz;
    313316    }
    314317    psF64 radius = sigma * sqrt (2.0 * z);
    … …  
    475478# undef PM_MODEL_FIT_STATUS
    476479# undef PM_MODEL_SET_LIMITS
     480# undef ALPHA
     481# undef ALPHA_M

  • branches/eam_branches/20100225/psModules/src/objects/models/pmModel_RGAUSS.c

    r26916 r27266  
    257257psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
    258258{
    259     psF64 z, f;
     259    psF64 z;
    260260    int Nstep = 0;
    261261    psEllipseShape shape;
    … …  
    291291    z0 = 0.0;
    292292
    293     // perform a type of bisection to find the value
    294     float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7]));
    295     float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7]));
    296     while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
    297         z = 0.5*(z0 + z1);
    298         f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
    299         if (f > limit) {
    300             z0 = z;
    301             f0 = f;
    302         } else {
    303             z1 = z;
    304             f1 = f;
    305         }
    306         Nstep ++;
    307     }
     293    // starting guess:
     294    z = 0.5*(z0 + z1);
     295    float dz = 1.0;
     296
     297    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
     298        // use Newton-Raphson to minimize f(z) - limit = 0
     299        float q = (1 + z + pow(z,PAR[PM_PAR_7]));
     300        float dqdz = (1.0 + PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0));
     301
     302        float f = 1.0 / q;
     303        float dfdz = -dqdz * f / q;
     304
     305        dz = (f - limit) / dfdz;
     306
     307        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
     308        z -= dz;
     309    }
     310
    308311    psF64 radius = sigma * sqrt (2.0 * z);
    309312
Note: See TracChangeset for help on using the changeset viewer.