IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 20, 2006, 3:34:56 AM (20 years ago)
Author:
rhl
Message:

Use symbolic param names

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9621 r9675  
    44
    55/******************************************************************************
    6     params->data.F32[0] = So;
    7     params->data.F32[1] = Zo;
    8     params->data.F32[2] = Xo;
    9     params->data.F32[3] = Yo;
    10     params->data.F32[4] = sqrt(2.0) / SigmaX;
    11     params->data.F32[5] = sqrt(2.0) / SigmaY;
    12     params->data.F32[6] = Sxy;
     6    params->data.F32[PM_PAR_SKY] = So;
     7    params->data.F32[PM_PAR_I0] = Zo;
     8    params->data.F32[PM_PAR_XPOS] = Xo;
     9    params->data.F32[PM_PAR_YPOS] = Yo;
     10    params->data.F32[PM_PAR_SXX] = sqrt(2.0) / SigmaX;
     11    params->data.F32[PM_PAR_SYY] = sqrt(2.0) / SigmaY;
     12    params->data.F32[PM_PAR_SXY] = Sxy;
    1313*****************************************************************************/
    1414
    … …  
    1717                        const psVector *x)
    1818{
    19     psF32 X  = x->data.F32[0] - params->data.F32[2];
    20     psF32 Y  = x->data.F32[1] - params->data.F32[3];
    21     psF32 px = params->data.F32[4]*X;
    22     psF32 py = params->data.F32[5]*Y;
    23     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     19    psF32 X  = x->data.F32[0] - params->data.F32[PM_PAR_XPOS];
     20    psF32 Y  = x->data.F32[1] - params->data.F32[PM_PAR_YPOS];
     21    psF32 px = params->data.F32[PM_PAR_SXX]*X;
     22    psF32 py = params->data.F32[PM_PAR_SYY]*Y;
     23    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[PM_PAR_SXY]*X*Y;
    2424    psF32 r  = exp(-z);
    25     psF32 q  = params->data.F32[1]*r;
    26     psF32 f  = q + params->data.F32[0];
     25    psF32 q  = params->data.F32[PM_PAR_I0]*r;
     26    psF32 f  = q + params->data.F32[PM_PAR_SKY];
    2727
    2828    if (deriv != NULL) {
    29         deriv->data.F32[0] = +1.0;
    30         deriv->data.F32[1] = +r;
    31         deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
    32         deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
    33         deriv->data.F32[4] = -2.0*q*px*X;
    34         deriv->data.F32[5] = -2.0*q*py*Y;
    35         deriv->data.F32[6] = -q*X*Y;
     29        deriv->data.F32[PM_PAR_SKY] = +1.0;
     30        deriv->data.F32[PM_PAR_I0] = +r;
     31        deriv->data.F32[PM_PAR_XPOS] = q*(2*px*params->data.F32[PM_PAR_SXX] + params->data.F32[PM_PAR_SXY]*Y);
     32        deriv->data.F32[PM_PAR_YPOS] = q*(2*py*params->data.F32[PM_PAR_SYY] + params->data.F32[PM_PAR_SXY]*X);
     33        deriv->data.F32[PM_PAR_SXX] = -2.0*q*px*X;
     34        deriv->data.F32[PM_PAR_SYY] = -2.0*q*py*Y;
     35        deriv->data.F32[PM_PAR_SXY] = -q*X*Y;
    3636    }
    3737    return(f);
    … …  
    4949    (*params_max)->n = (*params_max)->nalloc;
    5050
    51     beta_lim[0][0].data.F32[0] = 1000;
    52     beta_lim[0][0].data.F32[1] = 3e6;
    53     beta_lim[0][0].data.F32[2] = 5;
    54     beta_lim[0][0].data.F32[3] = 5;
    55     beta_lim[0][0].data.F32[4] = 0.5;
    56     beta_lim[0][0].data.F32[5] = 0.5;
    57     beta_lim[0][0].data.F32[6] = 0.5;
     51    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
     52    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
     53    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
     54    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
     55    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
     56    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
     57    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
    5858
    59     params_min[0][0].data.F32[0] = -1000;
    60     params_min[0][0].data.F32[1] = 0;
    61     params_min[0][0].data.F32[2] = -100;
    62     params_min[0][0].data.F32[3] = -100;
    63     params_min[0][0].data.F32[4] = 0.01;
    64     params_min[0][0].data.F32[5] = 0.01;
    65     params_min[0][0].data.F32[6] = -5.0;
     59    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
     60    params_min[0][0].data.F32[PM_PAR_I0] = 0;
     61    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
     62    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
     63    params_min[0][0].data.F32[PM_PAR_SXX] = 0.01;
     64    params_min[0][0].data.F32[PM_PAR_SYY] = 0.01;
     65    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    6666
    67     params_max[0][0].data.F32[0] = 1e5;
    68     params_max[0][0].data.F32[1] = 1e8;
    69     params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    70     params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    71     params_max[0][0].data.F32[4] = 2.0;
    72     params_max[0][0].data.F32[5] = 2.0;
    73     params_max[0][0].data.F32[6] = +5.0;
     67    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
     68    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
     69    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
     70    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
     71    params_max[0][0].data.F32[PM_PAR_SXX] = 2.0;
     72    params_max[0][0].data.F32[PM_PAR_SYY] = 2.0;
     73    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    7474
    7575    return (TRUE);
    … …  
    8383    psF32     *params  = model->params->data.F32;
    8484
    85     params[0] = moments->Sky;
    86     params[1] = moments->Peak - moments->Sky;
    87     params[2] = moments->x;
    88     params[3] = moments->y;
    89     params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
    90     params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
    91     params[6] = 0.0;
     85    params[PM_PAR_SKY] = moments->Sky;
     86    params[PM_PAR_I0] = moments->Peak - moments->Sky;
     87    params[PM_PAR_XPOS] = moments->x;
     88    params[PM_PAR_YPOS] = moments->y;
     89    params[PM_PAR_SXX] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
     90    params[PM_PAR_SYY] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
     91    params[PM_PAR_SXY] = 0.0;
    9292    return(true);
    9393}
    … …  
    9595psF64 pmModelFlux_GAUSS(const psVector *params)
    9696{
    97     psF64 A1   = PS_SQR(params->data.F32[4]);
    98     psF64 A2   = PS_SQR(params->data.F32[5]);
    99     psF64 A3   = PS_SQR(params->data.F32[6]);
     97    psF64 A1   = PS_SQR(params->data.F32[PM_PAR_SXX]);
     98    psF64 A2   = PS_SQR(params->data.F32[PM_PAR_SYY]);
     99    psF64 A3   = PS_SQR(params->data.F32[PM_PAR_SXY]);
    100100    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    101101    // Area is equivalent to 2 pi sigma^2
    102102
    103     psF64 Flux = params->data.F32[1] * Area;
     103    psF64 Flux = params->data.F32[PM_PAR_I0] * Area;
    104104
    105105    return(Flux);
    … …  
    113113    if (flux <= 0)
    114114        return (1.0);
    115     if (params->data.F32[1] <= 0)
     115    if (params->data.F32[PM_PAR_I0] <= 0)
    116116        return (1.0);
    117     if (flux >= params->data.F32[1])
     117    if (flux >= params->data.F32[PM_PAR_I0])
    118118        return (1.0);
    119119
    120     psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
    121     psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
     120    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[PM_PAR_SXX], 1.0 / params->data.F32[PM_PAR_SYY]);
     121    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[PM_PAR_I0] / flux));
    122122    return (radius);
    123123}
    … …  
    130130    psF32 *in  = modelFLT->params->data.F32;
    131131
    132     out[0] = in[0];
    133     out[1] = in[1];
    134     out[2] = in[2];
    135     out[3] = in[3];
     132    out[PM_PAR_SKY] = in[PM_PAR_SKY];
     133    out[PM_PAR_I0] = in[PM_PAR_I0];
     134    out[PM_PAR_XPOS] = in[PM_PAR_XPOS];
     135    out[PM_PAR_YPOS] = in[PM_PAR_YPOS];
    136136
     137    assert (PM_PAR_YPOS + 1 == 4);  // so starting at 4 is correct
    137138    for (int i = 4; i < 7; i++) {
    138139        psPolynomial2D *poly = psf->params->data[i-4];
    139         out[i] = psPolynomial2DEval(poly, out[2], out[3]);
     140        out[i] = psPolynomial2DEval(poly, out[PM_PAR_XPOS], out[PM_PAR_YPOS]);
    140141    }
    141142    return(true);
    … …  
    153154
    154155    dP = 0;
    155     dP += PS_SQR(dPAR[4] / PAR[4]);
    156     dP += PS_SQR(dPAR[5] / PAR[5]);
     156    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
     157    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    157158    dP = sqrt (dP);
    158159
    159160    status = true;
    160161    status &= (dP < 0.5);
    161     status &= (PAR[1] > 0);
    162     status &= ((dPAR[1]/PAR[1]) < 0.5);
     162    status &= (PAR[PM_PAR_I0] > 0);
     163    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
    163164
    164165    if (status)
Note: See TracChangeset for help on using the changeset viewer.