IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 22, 2007, 7:42:46 AM (19 years ago)
Author:
eugene
Message:

finished psastroModel for boresite, created gpcModel for rough theoretical model, added useModel and fixChips options

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroModelFitBoresite.c

    r15887 r15891  
    22
    33// we now have a set of observed L,M values.  fit these to the boresite model
    4 bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
     4psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po) {
    55
    6     assert (Lo->n > 2);
    7     assert (Lo->n == Mo->n);
    8     assert (Lo->n == Po->n);
     6    assert (Xo->n > 2);
     7    assert (Xo->n == Yo->n);
     8    assert (Xo->n == Po->n);
    99
    1010    // arrays to hold the data to be fitted
    11     psArray *x = psArrayAlloc(2*Lo->n);
    12     psVector *y = psVectorAlloc(2*Lo->n, PS_TYPE_F32);
     11    psArray *x = psArrayAlloc(2*Xo->n);
     12    psVector *y = psVectorAlloc(2*Xo->n, PS_TYPE_F32);
    1313
    1414    int n = 0;
    15     for (int i = 0; i < Lo->n; i++) {
     15    for (int i = 0; i < Xo->n; i++) {
    1616
    1717        psVector *coord = NULL;
    1818
    19         // L coordinate value
     19        // X coordinate value
    2020        coord = psVectorAlloc (2, PS_TYPE_F32);
    2121        coord->data.F32[1] = 0.0;
    2222        coord->data.F32[0] = Po->data.F32[i];
    2323        x->data[n] = coord;
    24         y->data.F32[n] = Lo->data.F32[i];
     24        y->data.F32[n] = Xo->data.F32[i];
    2525        n++;
    2626       
    27         // M coordinate value
     27        // Y coordinate value
    2828        coord = psVectorAlloc (2, PS_TYPE_F32);
    2929        coord->data.F32[1] = 1.0;
    3030        coord->data.F32[0] = Po->data.F32[i];
    3131        x->data[n] = coord;
    32         y->data.F32[n] = Mo->data.F32[i];
     32        y->data.F32[n] = Yo->data.F32[i];
    3333        n++;
    3434    }   
    … …  
    4747    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    4848
    49     // center (Lo) = mean(Lo), RL = range / 2
    50     psVectorStats (stats, Lo, NULL, NULL, 0);
    51     params->data.F32[PAR_L0] = stats->sampleMean;
    52     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
     49    // center (Xo) = mean(Xo), RX = range / 2
     50    psVectorStats (stats, Xo, NULL, NULL, 0);
     51    params->data.F32[PAR_X0] = stats->sampleMean;
     52    params->data.F32[PAR_RX] = (stats->max - stats->min) / 2.0;
    5353
    54     // center (Mo) = mean(Mo), RM = range / 2
    55     psVectorStats (stats, Mo, NULL, NULL, 0);
    56     params->data.F32[PAR_M0] = stats->sampleMean;
    57     params->data.F32[PAR_RM] = (stats->max - stats->min) / 2.0;
     54    // center (Yo) = mean(Yo), RY = range / 2
     55    psVectorStats (stats, Yo, NULL, NULL, 0);
     56    params->data.F32[PAR_Y0] = stats->sampleMean;
     57    params->data.F32[PAR_RY] = (stats->max - stats->min) / 2.0;
    5858
    5959    params->data.F32[PAR_P0] = 0.0;
    … …  
    6363    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    6464   
    65     fprintf (stderr, "guess values:\n");
    66     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    67     fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
    68     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    69     fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
    70     fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
    71     fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
     65    // fprintf (stderr, "guess values:\n");
     66    // fprintf (stderr, "Xo:  %f\n", params->data.F32[PAR_X0]);
     67    // fprintf (stderr, "Yo:  %f\n", params->data.F32[PAR_Y0]);
     68    // fprintf (stderr, "RX:  %f\n", params->data.F32[PAR_RX]);
     69    // fprintf (stderr, "RY:  %f\n", params->data.F32[PAR_RY]);
     70    // fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
     71    // fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
    7272
    7373    // XXX skip the weights for now
    … …  
    7575
    7676    fprintf (stderr, "fitted values:\n");
    77     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    78     fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
    79     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    80     fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
     77    fprintf (stderr, "Xo:  %f\n", params->data.F32[PAR_X0]);
     78    fprintf (stderr, "Yo:  %f\n", params->data.F32[PAR_Y0]);
     79    fprintf (stderr, "RX:  %f\n", params->data.F32[PAR_RX]);
     80    fprintf (stderr, "RY:  %f\n", params->data.F32[PAR_RY]);
    8181    fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
    8282    fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
    8383
    84     return true;
     84    return params;
    8585}
    86 
    87 # if (0)
    88 # undef PAR_L0
    89 # undef PAR_RL
    90 # undef PAR_PHI
    91 # define PAR_L0  0
    92 # define PAR_RL  1
    93 # define PAR_PHI 2
    94 
    95 // we now have a set of observed L,M values.  fit these to the boresite model
    96 bool psastroModelFitBoresite_L (psVector *Lo, psVector *Mo, psVector *Po) {
    97 
    98     assert (Lo->n > 2);
    99     assert (Lo->n == Mo->n);
    100     assert (Lo->n == Po->n);
    101 
    102     // arrays to hold the data to be fitted
    103     psArray *x = psArrayAlloc(Lo->n);
    104     psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
    105 
    106     for (int i = 0; i < Lo->n; i++) {
    107 
    108         psVector *coord = NULL;
    109 
    110         // L coordinate value
    111         coord = psVectorAlloc (1, PS_TYPE_F32);
    112         coord->data.F32[0] = Po->data.F32[i];
    113         x->data[i] = coord;
    114         y->data.F32[i] = Lo->data.F32[i];
    115     }   
    116 
    117     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
    118    
    119     // create the minimization constraints
    120     psMinConstraint *constraint = psMinConstraintAlloc();
    121 
    122     // XXX for now, no parameter masks, skip checkLimits
    123     // constraint->checkLimits = psastroModelBoresiteLimits;
    124 
    125     // make an initial guess:
    126     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    127 
    128     // center (Lo) = mean(Lo), RL = range / 2
    129     psVectorStats (stats, Lo, NULL, NULL, 0);
    130     params->data.F32[PAR_L0] = stats->sampleMean;
    131     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
    132     params->data.F32[PAR_PHI] = 0.0;
    133 
    134     psMinimization *myMin = psMinimizationAlloc (15, 0.001);
    135     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    136    
    137     fprintf (stderr, "guess values:\n");
    138     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    139     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    140     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    141 
    142     // XXX skip the weights for now
    143     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteL);
    144 
    145     fprintf (stderr, "fitted values:\n");
    146     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    147     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    148     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    149 
    150     return true;
    151 }
    152 
    153 // we now have a set of observed L,M values.  fit these to the boresite model
    154 bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
    155 
    156     assert (Lo->n > 2);
    157     assert (Lo->n == Mo->n);
    158     assert (Lo->n == Po->n);
    159 
    160     // arrays to hold the data to be fitted
    161     psArray *x = psArrayAlloc(Lo->n);
    162     psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
    163 
    164     for (int i = 0; i < Lo->n; i++) {
    165 
    166         psVector *coord = NULL;
    167 
    168         // L coordinate value
    169         coord = psVectorAlloc (1, PS_TYPE_F32);
    170         coord->data.F32[0] = Po->data.F32[i];
    171         x->data[i] = coord;
    172         y->data.F32[i] = Mo->data.F32[i];
    173     }   
    174 
    175     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
    176    
    177     // create the minimization constraints
    178     psMinConstraint *constraint = psMinConstraintAlloc();
    179 
    180     // XXX for now, no parameter masks, skip checkLimits
    181     // constraint->checkLimits = psastroModelBoresiteLimits;
    182 
    183     // make an initial guess:
    184     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    185 
    186     // center (Lo) = mean(Lo), RL = range / 2
    187     psVectorStats (stats, Mo, NULL, NULL, 0);
    188     params->data.F32[PAR_L0] = stats->sampleMean;
    189     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
    190     params->data.F32[PAR_PHI] = 0.0;
    191 
    192     psMinimization *myMin = psMinimizationAlloc (15, 0.001);
    193     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    194    
    195     fprintf (stderr, "guess values:\n");
    196     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    197     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    198     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    199 
    200     // XXX skip the weights for now
    201     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteM);
    202 
    203     fprintf (stderr, "fitted values:\n");
    204     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    205     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    206     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    207 
    208     return true;
    209 }
    210 # endif
Note: See TracChangeset for help on using the changeset viewer.