IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 18, 2009, 6:24:53 AM (17 years ago)
Author:
eugene
Message:

updates for petrosian analysis study

Location:
branches/eam_branches/20090715/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy

  • branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c

    r25032 r25105  
    55// model parameters
    66enum {PAR_PHI, PAR_EPSILON, PAR_RMIN};
     7psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
    78
    89bool psphotEllipticalContour (pmPetrosian *petrosian) {
    … …  
    1011    // use LMM to fit theta vs radius to an ellipse
    1112    psVector *theta = petrosian->theta;
    12     psVector *radius = petrosian->isophotalRadius;
     13    psVector *radius = petrosian->isophotalRadii;
    1314
    1415    // find Rmin and Rmax for the initial guess
    … …  
    2021    psArray *x = psArrayAlloc(2*radius->n);
    2122    psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32);
     23    psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32);
    2224
     25    int n = 0;
    2326    for (int i = 0; i < radius->n; i++) {
    2427
    … …  
    2831        coord = psVectorAlloc (2, PS_TYPE_F32);
    2932        coord->data.F32[1] = 0.0;
    30         coord->data.F32[0] = radius->data.F32[i]*cos(theta->data.F32[i]);
     33        coord->data.F32[0] = theta->data.F32[i];
    3134        x->data[n] = coord;
    32         y->data.F32[n] = theta->data.F32[i];
     35        y->data.F32[n] = radius->data.F32[i]*cos(theta->data.F32[i]);
     36        yErr->data.F32[n] = 1000.0;
    3337        n++;
    3438
    … …  
    3640        coord = psVectorAlloc (2, PS_TYPE_F32);
    3741        coord->data.F32[1] = 1.0;
    38         coord->data.F32[0] = radius->data.F32[i]*sin(theta->data.F32[i]);
     42        coord->data.F32[0] = theta->data.F32[i];
    3943        x->data[n] = coord;
    40         y->data.F32[n] = theta->data.F32[i];
     44        y->data.F32[n] = radius->data.F32[i]*sin(theta->data.F32[i]);
     45        yErr->data.F32[n] = 1000.0;
    4146        n++;
    4247
    … …  
    4954
    5055    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
     56   
     57    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 7);
    5158   
    5259    // create the minimization constraints
    … …  
    6471   
    6572    // XXX skip the weights for now
    66     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psphotEllipticalContourFunc);
     73    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc);
    6774
    6875    fprintf (stderr, "# fitted values:\n");
    69     fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]);
     76    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
    7077    fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
    7178    fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
    … …  
    7582    petrosian->axes.minor = params->data.F32[PAR_RMIN];
    7683    petrosian->axes.theta = params->data.F32[PAR_PHI];
     84
     85    // show the results
     86    psphotPetrosianVisualEllipticalContour (petrosian);
     87
     88    psFree (x);
     89    psFree (y);
     90    psFree (yErr);
    7791
    7892    return true;
    … …  
    88102psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
    89103
     104    static int pass = 0;
     105
    90106    psF32 *par = params->data.F32;
     107
     108    float alpha = coord->data.F32[0];
    91109
    92110    float cs_alpha = cos(alpha);
    93111    float sn_alpha = sin(alpha);
    94112
    95     float alpha = coord->data.F32[0];
    96113    float cs_phi = cos(alpha - par[PAR_PHI]);
    97114    float sn_phi = sin(alpha - par[PAR_PHI]);
    … …  
    103120
    104121    // value is X
    105     if (coord->data.F32[1] == 0) {
     122    // if (coord->data.F32[1] == 0) {
     123    if (pass == 0) {
     124        pass = 1;
    106125
    107126        float value = par[PAR_RMIN]*cs_alpha*r;
    … …  
    109128        if (deriv) {
    110129            psF32 *dPAR = deriv->data.F32;
    111             dpar[PAR_RMIN]    = r*cs_alpha;
    112             dpar[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE;
    113             dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*cs_alpha*drdP;
     130            dPAR[PAR_RMIN]    = r*cs_alpha;
     131            dPAR[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE;
     132            dPAR[PAR_PHI]     = 4.0*par[PAR_RMIN]*cs_alpha*drdP;
    114133        }
    115134        return (value);
    … …  
    117136
    118137    // value is Y
    119     if (coord->data.F32[1] == 1) {
     138    // if (coord->data.F32[1] == 1) {
     139    if (pass == 1) {
     140        pass = 0;
    120141
    121142        float value = par[PAR_RMIN]*sn_alpha*r;
    … …  
    123144        if (deriv) {
    124145            psF32 *dPAR = deriv->data.F32;
    125             dpar[PAR_RMIN]    = r*sn_alpha;
    126             dpar[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE;
    127             dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*sn_alpha*drdP;
     146            dPAR[PAR_RMIN]    = r*sn_alpha;
     147            dPAR[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE;
     148            dPAR[PAR_PHI]     = 4.0*par[PAR_RMIN]*sn_alpha*drdP;
    128149        }
    129150        return (value);
Note: See TracChangeset for help on using the changeset viewer.