IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38599


Ignore:
Timestamp:
Jul 17, 2015, 8:39:34 PM (11 years ago)
Author:
eugene
Message:

working on relastro

Location:
branches/eam_branches/ipp-20150625/Ohana/src
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d

    r38441 r38599  
    6969FIELD catID,          CAT_ID,      unsigned int,    unique ID for table in which object was first realized
    7070FIELD extID,          EXT_ID,      uint64_t,        external ID for object (eg PSPS objID)
    71 FIELD extIDgc,        EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
     71
     72# replace extIDgc (unused) with uRgal, uDgal:
     73# FIELD extIDgc,      EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
     74FIELD uRgal,          U_RA_GAL,    float,           modeled proper motion based on galactic motion
     75FIELD uDgal,          U_DEC_GAL,   float,           modeled proper motion based on galactic motion
    7276
    7377# this structure should only be used for internal representations

  • branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/measure.d

    r38062 r38599  
    4242FIELD XoffCAM,        X_OFF_CAM,     float,          X offset from correction,     pixels
    4343FIELD YoffCAM,        Y_OFF_CAM,     float,          Y offset from correction,     pixels
     44
     45# XXX I can deprecate these as I am going to apply the correct uR,uD offset
     46# XXX not sure how to use this yet...
    4447FIELD RoffGAL,        R_OFF_GAL,     float,          RA offset from correction,    arcsec
    4548FIELD DoffGAL,        D_OFF_GAL,     float,          DEC offset from correction,   arcsec

  • branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h

    r38590 r38599  
    171171  ID_PROPER            = 0x00000010, // star with large proper motion
    172172  ID_TRANSIENT         = 0x00000020, // identified as a non-periodic (stationary) transient
    173   ID_VARIABLE          = 0x00000040, // identified as a period variable
     173  ID_VARIABLE          = 0x00000040, // identified as a periodic variable
    174174  ID_ASTEROID          = 0x00000080, // identified with a known solar-system object (asteroid or other)
    175   // bits 0x00000100 - 0x00008000 are currently unused
     175  ID_STACK_ASTROM      = 0x00000100, // stack position used for astrometry
     176 
     177  // bits 0x00000200 - 0x00008000 are currently unused
    176178  ID_STAR_FIT_AVE      = 0x00010000, // average position fitted
    177179  ID_STAR_FIT_PM       = 0x00020000, // proper motion fitted

  • branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_catalog_create.c

    r38553 r38599  
    3333    catalog[0].lensobj_catalog  = dvo_catalog_create_subcat (catalog, "cpy", "LENSOBJ");
    3434    catalog[0].starpar_catalog  = dvo_catalog_create_subcat (catalog, "cpz", "STARPAR");
    35     catalog[0].galphot_catalog = dvo_catalog_create_subcat (catalog, "cpq", "GALPHOT");
     35    catalog[0].galphot_catalog  = dvo_catalog_create_subcat (catalog, "cpq", "GALPHOT");
    3636
    3737    // lock the additional split files

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h

    r38441 r38599  
    115115  double Ro, dRo;
    116116  double Do, dDo;
    117 
    118117  double uR, duR;
    119118  double uD, duD;
    120 
    121   double p, dp;
    122 
     119  double  p, dp;
     120
     121  int getChisq;
    123122  double chisq;
    124123  int Nfit;
    125 } PMFit;
     124} FitAstromResult;
     125
     126typedef struct {
     127  double **A;
     128  double **B;
     129  int Nterms;
     130} FitAstromData;
     131
     132typedef struct {
     133  double X, dX;
     134  double Y, dY;
     135  double R, dR;
     136  double D, dD;
     137  double T, dT;
     138  double pX;
     139  double pY;
     140  double C_blue;
     141  double C_red;
     142} FitAstromPoints;
     143
     144typedef struct {
     145  off_t Nave;
     146  off_t Npm;
     147  off_t Npar;
     148  off_t Nskip;
     149  off_t Noffset;
     150
     151  FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats
     152  int Nfit;
     153  int NfitAlloc;
     154
     155  FitAstromPoints *points;
     156  FitAstromPoints *sample;
     157  int Npoints;
     158  int NpointsAlloc;
     159
     160  FitAstromData *fitPos;
     161  FitAstromData *fitPM;
     162  FitAstromData *fitPar;
     163
     164  Coords coords;
     165  time_t T2000;
     166} FitStats;
    126167
    127168typedef struct {
    … …  
    148189  int    Nmeas;
    149190} StatType;
    150 
    151 typedef struct {
    152   off_t Nave;
    153   off_t Npm;
    154   off_t Npar;
    155   off_t Nskip;
    156   off_t Noffset;
    157 } FitStats;
    158191
    159192typedef struct {

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c

    r32695 r38599  
    11# include "relastro.h"
    22
    3 /* do we want an init function which does the alloc and a clear function to free? */
    4 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB) {
     3int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    54
    65  int i;
    76
    8   double **A, **B;
    97  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    10   double chisq, Xf, Yf;
    118
    12   /* do I need to do this as 2 2x2 matrix equations? */
    13   A = array_init (4, 4);
    14   B = array_init (4, 1);
     9  myAssert (data->Nterms == 4, "invalid fit arrays");
    1510
    1611  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    17   for (i = 0; i < Npts; i++) {
     12
     13  for (i = 0; i < Npoints; i++) {
    1814    /* handle case where dX or dY = 0.0 */
    19     wx = 1.0 / SQ(dX[i]);
    20     wy = 1.0 / SQ(dY[i]);
     15    wx = 1.0 / SQ(points[i].dX);
     16    wy = 1.0 / SQ(points[i].dY);
    2117
    2218    Wx += wx;
    2319    Wy += wy;
    2420
    25     Tx += T[i]*wx;
    26     Ty += T[i]*wy;
     21    Tx += points[i].T*wx;
     22    Ty += points[i].T*wy;
    2723   
    28     Tx2 += SQ(T[i])*wx;
    29     Ty2 += SQ(T[i])*wy;
     24    Tx2 += SQ(points[i].T)*wx;
     25    Ty2 += SQ(points[i].T)*wy;
    3026   
    31     Xs += X[i]*wx;
    32     Ys += Y[i]*wy;
     27    Xs += points[i].X*wx;
     28    Ys += points[i].Y*wy;
    3329
    34     XT += X[i]*T[i]*wx;
    35     YT += Y[i]*T[i]*wy;
     30    XT += points[i].X*points[i].T*wx;
     31    YT += points[i].Y*points[i].T*wy;
    3632  }
    3733
    38   A[0][0] = Wx;
    39   A[0][1] = Tx;
     34  data->A[0][0] = Wx;
     35  data->A[0][1] = Tx;
    4036
    41   A[1][0] = Tx;
    42   A[1][1] = Tx2;
     37  data->A[1][0] = Tx;
     38  data->A[1][1] = Tx2;
    4339
    44   A[2][2] = Wy;
    45   A[2][3] = Ty;
     40  data->A[2][2] = Wy;
     41  data->A[2][3] = Ty;
    4642
    47   A[3][2] = Ty;
    48   A[3][3] = Ty2;
     43  data->A[3][2] = Ty;
     44  data->A[3][3] = Ty2;
    4945
    50   B[0][0] = Xs;
    51   B[1][0] = XT;
    52   B[2][0] = Ys;
    53   B[3][0] = YT;
     46  data->B[0][0] = Xs;
     47  data->B[1][0] = XT;
     48  data->B[2][0] = Ys;
     49  data->B[3][0] = YT;
    5450
    55   dgaussjordan (A, B, 4, 1);
     51  dgaussjordan (data->A, data->B, 4, 1);
    5652
    57   fit[0].Ro = B[0][0];
    58   fit[0].uR = B[1][0];
    59   fit[0].Do = B[2][0];
    60   fit[0].uD = B[3][0];
    61   fit[0].p  = 0.0;
     53  fit->Ro = data->B[0][0];
     54  fit->uR = data->B[1][0];
     55  fit->Do = data->B[2][0];
     56  fit->uD = data->B[3][0];
     57  fit->p  = 0.0;
    6258 
    63   fit[0].dRo = sqrt(A[0][0]);
    64   fit[0].duR = sqrt(A[1][1]);
    65   fit[0].dDo = sqrt(A[2][2]);
    66   fit[0].duD = sqrt(A[3][3]);
    67   fit[0].dp  = 0.0;
     59  fit->dRo = sqrt(data->A[0][0]);
     60  fit->duR = sqrt(data->A[1][1]);
     61  fit->dDo = sqrt(data->A[2][2]);
     62  fit->duD = sqrt(data->A[3][3]);
     63  fit->dp  = 0.0;
    6864 
    69   array_free (A, 4);
    70   array_free (B, 4);
     65  fit->Nfit = Npoints;
    7166
    7267  // add up the chi square for the fit
    73   chisq = 0.0;
    74   for (i = 0; i < Npts; i++) {
    75     Xf = fit[0].Ro + fit[0].uR*T[i];
    76     Yf = fit[0].Do + fit[0].uD*T[i];
    77     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
    78     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
    79     if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
     68  if (fit->getChisq) {
     69    double chisq = 0.0;
     70    for (i = 0; i < Npoints; i++) {
     71      double Xf = fit->Ro + fit->uR*points[i].T;
     72      double Yf = fit->Do + fit->uD*points[i].T;
     73      chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
     74      chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
     75    }
     76    // the reduced chisq is divided by (Ndof = 2*Npts - 4)
     77    fit->chisq = chisq / (2.0*Npts - 4.0);
    8078  }
    81   fit[0].Nfit = Npts;
    8279
    83   // the reduced chisq is divided by (Ndof = 2*Npts - 4)
    84   fit[0].chisq = chisq / (2.0*Npts - 4.0);
    8580  return (TRUE);
    8681}
    87 
    88 // XXX this function should (optionally?) iterate and clip outlier detections

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c

    r32695 r38599  
    11# include "relastro.h"
    22
    3 /* do we want an init function which does the alloc and a clear function to free? */
    4 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB) {
     3int FitPMandPAR (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    54
    65  int i;
    76
    8   double **A, **B;
    97  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    108  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
    11   double chisq, Xf, Yf;
    129
    13   A = array_init (5, 5);
    14   B = array_init (5, 1);
     10  myAssert (data->Nterms == 5, "invalid fit arrays");
    1511
    1612  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
    1713  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    18   for (i = 0; i < Npts; i++) {
     14
     15  for (i = 0; i < Npoints; i++) {
    1916    /* handle case where dX or dY = 0.0 */
    20     wx = 1.0 / SQ(dX[i]);
    21     wy = 1.0 / SQ(dY[i]);
     17    wx = 1.0 / SQ(points[i].dX);
     18    wy = 1.0 / SQ(points[i].dY);
    2219
    2320    Wx += wx;
    2421    Wy += wy;
    2522
    26     Tx += T[i]*wx;
    27     Ty += T[i]*wy;
     23    Tx += points[i].T*wx;
     24    Ty += points[i].T*wy;
    2825   
    29     Tx2 += SQ(T[i])*wx;
    30     Ty2 += SQ(T[i])*wy;
     26    Tx2 += SQ(points[i].T)*wx;
     27    Ty2 += SQ(points[i].T)*wy;
    3128   
    32     PR += pR[i]*wx;
    33     PD += pD[i]*wy;
     29    PR += points[i].pR*wx;
     30    PD += points[i].pD*wy;
    3431   
    35     PRT += pR[i]*T[i]*wx;
    36     PDT += pD[i]*T[i]*wy;
     32    PRT += points[i].pR*points[i].T*wx;
     33    PDT += points[i].pD*points[i].T*wy;
    3734   
    38     PRX += pR[i]*X[i]*wx;
    39     PDY += pD[i]*Y[i]*wy;
     35    PRX += points[i].pR*points[i].X*wx;
     36    PDY += points[i].pD*points[i].Y*wy;
    4037   
    41     PR2 += SQ(pR[i])*wx;
    42     PD2 += SQ(pD[i])*wy;
     38    PR2 += SQ(points[i].pR)*wx;
     39    PD2 += SQ(points[i].pD)*wy;
    4340
    44     Xs += X[i]*wx;
    45     Ys += Y[i]*wy;
     41    Xs += points[i].X*wx;
     42    Ys += points[i].Y*wy;
    4643
    47     XT += X[i]*T[i]*wx;
    48     YT += Y[i]*T[i]*wy;
     44    XT += points[i].X*points[i].T*wx;
     45    YT += points[i].Y*points[i].T*wy;
    4946  }
    5047
    51   A[0][0] = Wx;
    52   A[0][1] = Tx;
    53   A[0][4] = PR;
     48  data->A[0][0] = Wx;
     49  data->A[0][1] = Tx;
     50  data->A[0][4] = PR;
    5451
    55   A[1][0] = Tx;
    56   A[1][1] = Tx2;
    57   A[1][4] = PRT;
     52  data->A[1][0] = Tx;
     53  data->A[1][1] = Tx2;
     54  data->A[1][4] = PRT;
    5855
    59   A[2][2] = Wy;
    60   A[2][3] = Ty;
    61   A[2][4] = PD;
     56  data->A[2][2] = Wy;
     57  data->A[2][3] = Ty;
     58  data->A[2][4] = PD;
    6259
    63   A[3][2] = Ty;
    64   A[3][3] = Ty2;
    65   A[3][4] = PDT;
     60  data->A[3][2] = Ty;
     61  data->A[3][3] = Ty2;
     62  data->A[3][4] = PDT;
    6663
    67   A[4][0] = PR;
    68   A[4][1] = PRT;
    69   A[4][2] = PD;
    70   A[4][3] = PDT;
    71   A[4][4] = PR2 + PD2;
     64  data->A[4][0] = PR;
     65  data->A[4][1] = PRT;
     66  data->A[4][2] = PD;
     67  data->A[4][3] = PDT;
     68  data->A[4][4] = PR2 + PD2;
    7269
    73   B[0][0] = Xs;
    74   B[1][0] = XT;
    75   B[2][0] = Ys;
    76   B[3][0] = YT;
    77   B[4][0] = PRX + PDY;
     70  data->B[0][0] = Xs;
     71  data->B[1][0] = XT;
     72  data->B[2][0] = Ys;
     73  data->B[3][0] = YT;
     74  data->B[4][0] = PRX + PDY;
    7875
    79   dgaussjordan (A, B, 5, 1);
     76  dgaussjordan (data->A, data->B, 5, 1);
    8077
    81   fit[0].Ro = B[0][0];
    82   fit[0].uR = B[1][0];
    83   fit[0].Do = B[2][0];
    84   fit[0].uD = B[3][0];
    85   fit[0].p  = B[4][0];
     78  fit->Ro = data->B[0][0];
     79  fit->uR = data->B[1][0];
     80  fit->Do = data->B[2][0];
     81  fit->uD = data->B[3][0];
     82  fit->p  = data->B[4][0];
    8683 
    87   fit[0].dRo = sqrt(A[0][0]);
    88   fit[0].duR = sqrt(A[1][1]);
    89   fit[0].dDo = sqrt(A[2][2]);
    90   fit[0].duD = sqrt(A[3][3]);
    91   fit[0].dp  = sqrt(A[4][4]);
     84  fit->dRo = sqrt(data->A[0][0]);
     85  fit->duR = sqrt(data->A[1][1]);
     86  fit->dDo = sqrt(data->A[2][2]);
     87  fit->duD = sqrt(data->A[3][3]);
     88  fit->dp  = sqrt(data->A[4][4]);
    9289 
    93   array_free (A, 5);
    94   array_free (B, 5);
    95 
    96   /* get the chisq from the matrix values */
     90  fit->Nfit = Npoints;
    9791
    9892  // add up the chi square for the fit
    99   chisq = 0.0;
    100   for (i = 0; i < Npts; i++) {
    101     Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
    102     Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
    103     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
    104     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
    105     if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
     93  if (fit->getChisq) {
     94    double chisq = 0.0;
     95    for (i = 0; i < Npoints; i++) {
     96      double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR;
     97      double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD;
     98      chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
     99      chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
     100    }
     101    // the reduced chisq is divided by (Ndof = 2*Npts - 5)
     102    fit->chisq = chisq / (2.0*Npts - 5.0);
     103  }
    106104
    107   }
    108   fit[0].Nfit = Npts;
    109 
    110   // the reduced chisq is divided by (Ndof = 2*Npts - 5)
    111   fit[0].chisq = chisq / (2.0*Npts - 5.0);
    112105  return (TRUE);
    113106}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/GetAstromError.c

    r36833 r38599  
    11# include "relastro.h"
    22# define WEIGHTED_ERRORS 1
     3
     4// XXX hard-wire the trends identified by CZW
     5static float BrightMo[] = {-15.6, -16.8, -17.0, -16.7, -16.0};
     6static float BrightMs[] = {1.3, 1.3, 1.3, 1.8, 2.0};
    37
    48float GetAstromErrorTiny (MeasureTiny *measure, int mode) {
    … …  
    4044  dM    = measure[0].dM;
    4145  dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
     46
     47  // for GPC1 data, we have a bright end model:
     48  if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) {
     49    int Np = ((int) (measure[0].photcode / 100)) % 100;
     50    myAssert (Np >= 0, "oops");
     51    myAssert (Np <= 4, "oops");
     52
     53    float Minst = measure[0].M - measure[0].dt - 25.0;
     54    float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np])));
     55    dPtotal = hypot(dPtotal, dPbright);
     56  }
    4257
    4358  dPtotal = MAX (dPtotal, MIN_ERROR);
    … …  
    8499  dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
    85100
     101  // for GPC1 data, we have a bright end model:
     102  if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) {
     103    int Np = ((int) (measure[0].photcode / 100)) % 100;
     104    myAssert (Np >= 0, "oops");
     105    myAssert (Np <= 4, "oops");
     106
     107    float Minst = measure[0].M - measure[0].dt - 25.0;
     108    float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np])));
     109    dPtotal = hypot(dPtotal, dPbright);
     110  }
     111
    86112  dPtotal = MAX (dPtotal, MIN_ERROR);
    87113  return (dPtotal);

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ImageOps.c

    r38062 r38599  
    748748    n = measure[0].averef;
    749749
     750    // apply proper-motion from average position to measure epoch:
     751    float dTime = (measure[k].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years
     752
    750753    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    751754    ref[i].R = catalog[c].average[n].R;
    752755    ref[i].D = catalog[c].average[n].D;
     756   
     757    // XXX do this in a better way?
     758    ref[i].R += dTime * catalog[c].average[n].uR / 3600.0 / cos(ref[i].D*RAD_DEG);
     759    ref[i].D += dTime * catalog[c].average[n].uD / 3600.0;
    753760
    754761    // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ParFactor.c

    r37261 r38599  
    22# define J2000 2451545.       /* Julian date at standard epoch */
    33
    4 # if (0)
    5 /* Low precision formulae for the sun, from Almanac p. C24 (1990) */
    6 /* ra and dec are returned as decimal hours and decimal degrees. */
    7 void lpsun (double jd, double *ra, double *dec) {
     4/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
     5// jdoff is days since J2000
     6int sun_ecliptic (double jdoff, double *lambda, double *beta, double *epsilon, double *Radius) {
    87
    9   double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
    10 
    11   n = jd - J2000;
    12   L = 280.460 + 0.9856474 * n;
    13   g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
    14   lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
    15   epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
    16 
    17   // this is the conversion from ecliptic to celestial coords
    18   x = cos(lambda);
    19   y = cos(epsilon)*sin(lambda);
    20   z = sin(epsilon)*sin(lambda);
    21 
    22   *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
    23   *dec = (asin(z))*DEG_IN_RADIAN;
    24 }
    25 # endif
    26 
    27 # if (0)
    28 /* code borrowed from Skycalc : fix this stuff XXX */
    29 /* Low precision formulae for the sun, from Almanac p. C24 (1990) */
    30 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
    31 
    32   double n, L, g;
    33 
    34 
    35   n = jd - J2000;
    36   L = 280.460 + 0.9856474 * n;
    37   g = (357.528 + 0.9856003 * n)*RAD_DEG;
    38   *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
    39   *beta = 0.0;                                    // approx latitude
    40   *epsilon = (23.439 - 0.0000004 * n);            // obliquity of ecliptic in degrees
    41   return TRUE;
    42 }
    43 # endif
    44 
    45 /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
    46 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) {
    47 
    48   double n = jd - J2000;              // day number
     8  double n = jdoff;           // day number
    499  double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
    5010  double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly
    … …  
    5818
    5919/* given RA, DEC, Time, calculate the parallax factor */
    60 // Time is relative to Tmean, Tmean is years relative to J2000
    61 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) {
     20// Time is years since J2000
     21int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) {
    6222
    63   double jd, lambda, beta, epsilon, Radius;
     23  double lambda, beta, epsilon, Radius;
    6424
    65   /* given a Time relative to Tmean, Tmean in years since J2000, determine the solar
    66     longitude S */
     25  /* given a Time in years since J2000, determine the solar longitude S */
    6726
    68   // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
    69   jd = 365.25*(Time + Tmean) + J2000;
    70   // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd);
     27  double jdoff = 365.25*Time;
    7128
    72   sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius);
     29  sun_ecliptic (jdoff, &lambda, &beta, &epsilon, &Radius);
    7330
    7431  double lambda_rad = lambda*RAD_DEG;
    … …  
    9956  return TRUE;
    10057}
     58
     59# if (0)
     60/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
     61/* ra and dec are returned as decimal hours and decimal degrees. */
     62void lpsun (double jd, double *ra, double *dec) {
     63
     64  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
     65
     66  n = jd - J2000;
     67  L = 280.460 + 0.9856474 * n;
     68  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
     69  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
     70  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
     71
     72  // this is the conversion from ecliptic to celestial coords
     73  x = cos(lambda);
     74  y = cos(epsilon)*sin(lambda);
     75  z = sin(epsilon)*sin(lambda);
     76
     77  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
     78  *dec = (asin(z))*DEG_IN_RADIAN;
     79}
     80
     81/* code borrowed from Skycalc : fix this stuff XXX */
     82/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
     83int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
     84
     85  double n, L, g;
     86
     87
     88  n = jd - J2000;
     89  L = 280.460 + 0.9856474 * n;
     90  g = (357.528 + 0.9856003 * n)*RAD_DEG;
     91  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
     92  *beta = 0.0;                                    // approx latitude
     93  *epsilon = (23.439 - 0.0000004 * n);            // obliquity of ecliptic in degrees
     94  return TRUE;
     95}
     96# endif
     97

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c

    r37807 r38599  
    11# include "relastro.h"
    22# define PAR_TOOFEW 5
     3# define NBOOT 100
    34
    45int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
    56int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    6 
    7 static off_t Nmax;
    8 static double *X, *dX;
    9 static double *Y, *dY;
    10 static double *R, *dR;
    11 static double *D, *dD;
    12 static double *pX;
    13 static double *pY;
    14 static double *T;
    15 static double *dT;
    16 static double *C_blue;
    17 static double *C_red;
    18 
    19 static Coords coords;
    20 
    21 static time_t T2000;
    22 
    23 void initFitStats (FitStats *fitStats) {
    24   fitStats->Nave = 0; 
    25   fitStats->Npm = 0;   
    26   fitStats->Npar = 0; 
    27   fitStats->Nskip = 0;
    28   fitStats->Noffset = 0;
    29   return;
    30 }
    31 
    32 void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) {
    33   tgtFitStats->Nave    += srcFitStats->Nave    ; 
    34   tgtFitStats->Npm     += srcFitStats->Npm     ;   
    35   tgtFitStats->Npar    += srcFitStats->Npar    ; 
    36   tgtFitStats->Nskip   += srcFitStats->Nskip   ;
    37   tgtFitStats->Noffset += srcFitStats->Noffset ;
    38   return;
    39 }
    40 
    41 void initObjectData (Catalog *catalog, int Ncatalog) {
    42 
    43   off_t i, j;
    44  
    45   Nmax = 0;
    46   for (i = 0; i < Ncatalog; i++) {
    47     for (j = 0; j < catalog[i].Naverage; j++) {
    48       Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
    49     }
    50   }
    51 
    52   ALLOCATE (R, double, MAX (1, Nmax));
    53   ALLOCATE (D, double, MAX (1, Nmax));
    54   ALLOCATE (T, double, MAX (1, Nmax));
    55   ALLOCATE (X, double, MAX (1, Nmax));
    56   ALLOCATE (Y, double, MAX (1, Nmax));
    57 
    58   ALLOCATE (dR, double, MAX (1, Nmax));
    59   ALLOCATE (dD, double, MAX (1, Nmax));
    60   ALLOCATE (dT, double, MAX (1, Nmax));
    61   ALLOCATE (dX, double, MAX (1, Nmax));
    62   ALLOCATE (dY, double, MAX (1, Nmax));
    63 
    64   ALLOCATE (pX, double, MAX (1, Nmax));
    65   ALLOCATE (pY, double, MAX (1, Nmax));
    66 
    67   ALLOCATE (C_blue, double, MAX (1, Nmax));
    68   ALLOCATE (C_red,  double, MAX (1, Nmax));
    69 
    70   /* project coordinates to a plane centered on the object with units of arcsec */
    71   InitCoords (&coords, "DEC--SIN");
    72   coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    73 
    74   // use J2000 as a reference time
    75   T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
    76 } 
    77 
    78 void freeObjectData () {
    79 
    80   free (R);
    81   free (D);
    82   free (T);
    83   free (X);
    84   free (Y);
    85 
    86   free (dR);
    87   free (dD);
    88   free (dT);
    89   free (dX);
    90   free (dY);
    91 
    92   free (pX);
    93   free (pY);
    94 
    95   free (C_blue);
    96   free (C_red);
    97 } 
    987
    998// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    1009// be called with just MeasureTiny set and Measure == NULL
    10110int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) {
    102 
    103   initObjectData (catalog, Ncatalog);
    10411
    10512  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    … …  
    10916  }
    11017
    111   FitStats sumStatsChips; initFitStats (&sumStatsChips);
    112   FitStats sumStatsStack; initFitStats (&sumStatsStack);
     18  int NmeasureMax = CatalogMaxNmeasure (catalog, Ncatalog);
     19
     20  // allocate summary stats with Nmax = 0, Nboot = 0
     21  FitStats *sumStatsChips = FitStatsInit (0, 0);
     22  FitStats *sumStatsStack = FitStatsInit (0, 0);
     23
     24  FitStats *fitStatsChips = FitStatsInit (NmeasureMax, NBOOT);
     25  FitStats *fitStatsStack = FitStatsInit (NmeasureMax, NBOOT);
    11326
    11427  int i;
    … …  
    11730    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
    11831
    119     FitStats fitStatsChips; initFitStats (&fitStatsChips);
    120     FitStats fitStatsStack; initFitStats (&fitStatsStack);
     32    FitStatsReset (fitStatsChips);
     33    FitStatsReset (fitStatsStack);
    12134
    12235    off_t j;
    … …  
    12942      SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt];
    13043
    131       UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack);
    132       UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop);
     44      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack);
     45      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m, Nloop);
    13346    }
    13447    if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
    13548    if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
    136     sumFitStats (&fitStatsChips, &sumStatsChips);
    137     sumFitStats (&fitStatsStack, &sumStatsStack);
    138   }
    139   freeObjectData ();
     49    FitStatsSum (fitStatsChips, sumStatsChips);
     50    FitStatsSum (fitStatsStack, sumStatsStack);
     51  }
     52
     53  FitStatsFree (fitStatsChips);
     54  FitStatsFree (fitStatsStack);
    14055
    14156  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
    14257  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
     58
     59  FitStatsFree (sumStatsChips);
     60  FitStatsFree (sumStatsStack);
     61
    14362  return (TRUE);
    14463}
    … …  
    14665// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    14766// be called with just MeasureTiny set and Measure == NULL
    148 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) {
     67int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff) {
    14968
    15069  int setRefColor = areImagesMatched();
    … …  
    15473  PMFit fit;    memset (&fit,    0, sizeof(fit));
    15574  PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN;
    156   PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq = NAN;
     75  PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq  = NAN;
    15776  PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;
    15877
    … …  
    16988  int NcBlue = 0;
    17089  int NcRed = 0;
    171   int N = 0;
    17290
    17391  int mode = FIT_MODE; // start with the globally-defined fit mode
    … …  
    17795  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
    17896
    179   // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
    180   off_t k;
    181   for (k = 0; k < average[0].Nmeasure; k++) {
    182 
    183     if (XVERB) {
    184       char *date = ohana_sec_to_date (measure[k].t);
    185       int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
    186       fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
    187       free (date);
    188     }
    189 
    190     // SKIP gpc1 stack data
    191     if (isGPC1stack(measure[k].photcode)) continue;
    192 
    193     // SKIP gpc1 forced-warp data
    194     if (isGPC1warp(measure[k].photcode)) continue;
    195 
    196     // reset the bit to note that a detection was used (or not)
    197     measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    198     if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
    199 
    200     // does the measurement pass the supplied filtering constraints?
    201     // MeasFilterTestTiny does not test psfQF
    202     // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    203     int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    204     if (!keepMeasure) {
    205       continue;
    206     }
    207 
    208     double Ri = getMeanR (&measure[k], average, secfilt);
    209     double Di = getMeanD (&measure[k], average, secfilt);
    210 
    211     // if we are correcting for the Galaxy Motion Model, only should apply it here
    212     // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and
    213     // (b) if we think the measure R,D is already at the image epoch position
    214     if (USE_GALAXY_MODEL && applyGalaxyOffset) {
    215       Ri -= measure[k].RoffGAL / 3600.0;
    216       Di -= measure[k].DoffGAL / 3600.0;
    217     }
    218 
    219     // XXX add in dR,dD GAL here
    220 
    221     // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
    222     if (ExcludeBogus) {
    223       coords.crval1 = average[0].R;
    224       coords.crval2 = average[0].D;
    225       double Xi, Yi;
    226       RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
    227       double radius = hypot(Xi, Yi);
    228       if (radius > ExcludeBogusRadius) {
    229         measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    230         if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
    231         continue;
    232       }
    233       measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    234       if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
    235     }
    236 
    237     // outlier rejection
    238     if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
    239       continue;
    240     }
    241 
    242     R[N] = Ri;
    243     D[N] = Di;
    244 
    245     // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
    246     // T[] is time in years since J2000 (jd = 2451545)
    247     T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
    248 
    249     // dX, dY : error in arcsec --
    250     dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    251     dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
    252 
    253     // allow a given photcode or measurement to be
    254     // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    255     if (isnan(dX[N])) continue;
    256     if (isnan(dY[N])) continue;
    257 
    258     // add systematic error in quadrature, if desired
    259     // only do this after the fit has converged (or you will never improve the poor images)
    260     // if (INCLUDE_SYS_ERR) {
    261     // float dRsys = FromShortPixels(measure[k].dRsys);
    262     // dX[N] = hypot(dX[N], dRsys);
    263     // dY[N] = hypot(dY[N], dRsys);
    264     // }
    265 
    266     // dX[N] = 0.1;
    267     // dY[N] = 0.1;
    268 
    269     dT[N] = measure[k].dt;
    270 
    271     // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    272     // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    273     // non-circular errors (different values for X and Y), then dR and dD will be
    274     // incorrect: they would need to be rotated to take out the position angle
    275     dR[N] = dX[N] / 3600.0;
    276     dD[N] = dY[N] / 3600.0;
    277 
    278     if (setRefColor) {
    279       float colorBlue = getColorBlue (m+k, i);
    280       if (!isnan(colorBlue)) {
    281         C_blue[NcBlue] = colorBlue;
    282         NcBlue++;
    283       }
    284       float colorRed = getColorRed (m+k, i);
    285       if (!isnan(colorRed)) {
    286         C_red[NcRed] = colorRed;
    287         NcRed++;
    288       }
    289     }
    290 
    291     measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    292     if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
    293 
    294     N++;
    295   } // loop over measurements : average[0].Nmeasure
    296 
    297   if (N < 1) {
     97  // select the measurements to be used in this analysis
     98  FitAstromObject *points = UpdateObjects_SelectMeasures (average, secfilt, measure, measureBig, cat, measOff, &Npoints);
     99
     100  // if there are no exposure detections, use the stack position
     101  if (Npoints < 1) {
    298102    if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
    299103      average[0].R  = average[0].Rstk;
    … …  
    301105      average[0].dR = average[0].dRstk;
    302106      average[0].dD = average[0].dDstk;
     107      average[0].flags |= ID_STACK_ASTROM;
    303108    }
    304109    return FALSE;
    305110  }
    306111
    307   // if we have too few good detections for the desired fit, or too limited a
    308   // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
    309   // for even the average position, consider including the lower-quality detections?
    310 
    311   // find Tmin & Tmax from the list of accepted measurements
    312   double Tmean = 0.0;
    313   double Tmin = T[0];
    314   double Tmax = T[0];
    315   for (k = 0; k < N; k++) {
    316     Tmin = MIN(Tmin, T[k]);
    317     Tmax = MAX(Tmax, T[k]);
    318     Tmean += T[k];
    319   }
    320   double Trange = Tmax - Tmin;
    321 
    322   if (RELASTRO_OP == OP_HIGH_SPEED) {
    323     Tmean = 0.5*(Tmax - Tmin);
    324   } else {
    325     Tmean /= (float) N;
    326   }
    327 
    328   /* we need to do the fit in a locally linear space; choose a ref coordinate */
    329   coords.crval1 = R[0];
    330   coords.crval2 = D[0];
     112  double Tmean, Trange, parRange;
     113  UpdateObjects_Project (points, Npoints, &Tmean, &Trange, &parRange);
    331114
    332115  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    333116
    334   // project all of the R,D coordinates to a plane centered on this coordinate. set
    335   // the times to be relative to Tmean (this is required for parallax as well)
    336   for (k = 0; k < N; k++) {
    337     RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
    338     T[k] -= Tmean;
    339     if (XVERB) {
    340       fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
    341     }
    342   }       
     117  // if we have too few good detections for the desired fit, or too limited a baseline,
     118  // use a fit with fewer parameters.
    343119
    344120  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
    … …  
    346122    if (Trange < PM_DT_MIN) {
    347123      mode = FIT_AVERAGE;
    348       goto skipPM;
    349     }
    350     if (N <= PM_TOOFEW) {
     124      goto justPosition;
     125    }
     126    if (Npoints <= PM_TOOFEW) {
    351127      mode = FIT_AVERAGE;
    352       goto skipPM;
    353     }
    354 
    355     FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
     128      goto justPosition;
     129    }
     130
     131    fitStats->Nfit = 0;
     132    for (k = 0; k < fitStats->NfitAlloc; k++) {
     133      BootstrapResample (sample, points, Npoints);
     134      if (!FitPM (&fitStats->fit[k], fitdata, sample, Npoints)) continue;
     135      fitStats->Nfit ++;
     136    }
     137    // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
    356138
    357139    if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
    … …  
    374156  }
    375157 
    376 skipPM:
    377158  // fit the parallax + proper-motion model
    378159  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
    379160  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
    380161  if (mode == FIT_PM_AND_PAR) {
    381     if (Trange < PM_DT_MIN) {
    382       mode = FIT_PM_ONLY;
    383       goto skipPAR;
    384     }
    385162    if (N <= PAR_TOOFEW) {
    386163      mode = FIT_PM_ONLY;
    387       goto skipPAR;
    388     }
    389     float pXmin = +2.0;
    390     float pXmax = -2.0;
    391     float pYmin = +2.0;
    392     float pYmax = -2.0;
    393     for (k = 0; k < N; k++) {
    394       ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
    395       pXmin = MIN (pXmin, pX[k]);
    396       pXmax = MAX (pXmax, pX[k]);
    397       pYmin = MIN (pYmin, pY[k]);
    398       pYmax = MAX (pYmax, pY[k]);
    399     }
    400     float dXRange = pXmax - pXmin;
    401     float dYRange = pYmax - pYmin;
    402     float parRange = hypot (dXRange, dYRange);
    403        
     164      goto justPosition;
     165    }
    404166    if (parRange < PAR_FACTOR_MIN) {
    405167      mode = FIT_PM_ONLY;
    406       goto skipPAR;
    407     }
    408 
    409     FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
     168      goto justPosition;
     169    }
     170
     171    fitStats->Nfit = 0;
     172    for (k = 0; k < fitStats->NfitAlloc; k++) {
     173      BootstrapResample (sample, points, Npoints);
     174      FitPMandPar (&fitStats->fit[k], fitsdata, sample, Npoints);
     175      fitStats->Nfit ++;
     176    }
     177    // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
     178
    410179    if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);
    411180
    … …  
    414183    fitStats->Npar ++;
    415184
    416     if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     185    if (fabs(fitPAR.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    417186
    418187    // XXX a hard-wired hack...
    … …  
    422191  }       
    423192
    424 skipPAR:
     193justPosition:
    425194  {
    426     // ALWAYS fit the average model
    427     StatType statsR, statsD;
    428     liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    429     liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
     195    // use bootstrap resampling to check the error distribution
     196    // if we only have one point, this is silly...
     197   
     198    fitStats->Nfit = 0;
     199    for (k = 0; k < fitStats->NfitAlloc; k++) {
     200      BootstrapResample (sample, points, Npoints);
     201      FitPosPMfixed (&fitStats->fit[k], fitsdata, sample, Npoints);
     202      fitStats->Nfit ++;
     203    }
     204
     205    if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
    430206
    431207    // project Ro, Do back to RA,DEC
    432     XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords);
    433     if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
    434 
    435     fitAve.dRo = statsR.sigma;
    436     fitAve.dDo = statsD.sigma;
    437 
    438     fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
    439     fitAve.Nfit = N;
    440 
    441     fitAve.uR = fitAve.duR = 0.0;
    442     fitAve.uD = fitAve.duD = 0.0;
    443     fitAve.p  = fitAve.dp  = 0.0;
     208    XY_to_RD (&fitAve.Ro, &fitAve.Do, fitAve.Ro, fitAve.Do, &coords);
    444209    average[0].flags |= ID_STAR_FIT_AVE;
    445210    fitStats->Nave ++;
     211
     212    // XXX fitAve.dRo = statsR.sigma;
     213    // XXX fitAve.dDo = statsD.sigma;
     214    // XXX
     215    // XXX fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
     216    // XXX fitAve.Nfit = N;
     217    // XXX
     218    // XXX fitAve.uR = fitAve.duR = 0.0;
     219    // XXX fitAve.uD = fitAve.duD = 0.0;
     220    // XXX fitAve.p  = fitAve.dp  = 0.0;
     221  }
     222
     223  // XXX update the bit flags of which points were used
     224  for (k = 0; k < Npoints; k++) {
     225    int Nm = points[k].measure;
     226    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
     227    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
    446228  }
    447229
    448230  if (setRefColor) {
     231    // need to reassign here if isfinite()
    449232    float colorMedian;
    450233    dsort (C_blue, NcBlue);
    … …  
    696479}
    697480
    698 
     481int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) {
     482  int i;
     483
     484  // I need to draw Npoints random entries from 'points' with replacement:
     485  for (i = 0; i < Npoints; i++) {
     486    int N = Npoints * drand48();
     487    sample[i] = points[N];
     488  }
     489  return TRUE;
     490}
     491
     492FitAstromObject *UpdateObjects_SelectMeasures (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int *npoints) {
     493
     494  // XXX move this above, but just for re-working now:
     495  int Npoints = 0;
     496  int NPOINTS = average->Nmeasure;
     497  FitAstromObject *points = NULL;
     498  ALLOCATE (points, FitAstromObject, NPOINTS);
     499
     500  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
     501  off_t k;
     502  for (k = 0; k < average[0].Nmeasure; k++) {
     503
     504    if (XVERB) {
     505      char *date = ohana_sec_to_date (measure[k].t);
     506      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
     507      fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
     508      free (date);
     509    }
     510
     511    // SKIP gpc1 stack data
     512    if (isGPC1stack(measure[k].photcode)) continue;
     513
     514    // SKIP gpc1 forced-warp data
     515    if (isGPC1warp(measure[k].photcode)) continue;
     516
     517    // reset the bit to note that a detection was used (or not)
     518    measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
     519    if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
     520
     521    // does the measurement pass the supplied filtering constraints?
     522    // MeasFilterTestTiny does not test psfQF
     523    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
     524    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
     525    if (!keepMeasure) {
     526      continue;
     527    }
     528
     529    double Ri = measure[k].R;
     530    double Di = measure[k].D;
     531
     532    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
     533    if (ExcludeBogus) {
     534      coords.crval1 = average[0].R;
     535      coords.crval2 = average[0].D;
     536      double Xi, Yi;
     537      RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
     538      double radius = hypot(Xi, Yi);
     539      if (radius > ExcludeBogusRadius) {
     540        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
     541        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
     542        continue;
     543      }
     544      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     545      if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
     546    }
     547
     548    // outlier rejection
     549    if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
     550      continue;
     551    }
     552
     553    FitAstromObjectInit (&points[Npoints]);
     554
     555    points[Npoints].R = Ri;
     556    points[Npoints].D = Di;
     557
     558    // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
     559    // T[] is time in years since J2000 (jd = 2451545)
     560    points[Npoints].T = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
     561
     562    // add measured systematic error in quadrature?  only do this after the fit has
     563    // converged (or you will never improve the poor images)
     564
     565    // dX,dY are the X and Y direction errors in arcseconds.  dR, dD are the errors in
     566    // those directions in degrees.  IF we have non-circular errors (different values for
     567    // X and Y), then dR and dD will be incorrect: they would need to be rotated to take
     568    // out the position angle
     569
     570    // dX, dY : error in arcsec:
     571    points[Npoints].dX = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
     572    points[Npoints].dY = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
     573
     574    // allow a given photcode or measurement to be
     575    // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
     576    if (isnan(points[Npoints].dX)) continue;
     577    if (isnan(points[Npoints].dY)) continue;
     578
     579    points[Npoints].dT = measure[k].dt;
     580
     581    if (setRefColor) {
     582      points[Npoints].C_blue = getColorBlue (measOff+k, cat);
     583      points[Npoints].C_red  = getColorRed (measOff+k, cat);
     584    }
     585
     586    points[Npoints].measure = k;
     587    Npoints++;
     588  } // loop over measurements : average[0].Nmeasure
     589
     590  *npoints = Npoints;
     591  return points;
     592}
     593
     594int FitAstromObject_Project (FitAstromObject *points, int Npoints, double *Tmean, double *Trange, double *parRange) {
     595
     596  int k;
     597
     598  // find Tmin & Tmax from the list of accepted measurements
     599  double Tmin  = points[0].T;
     600  double Tmax  = points[0].T;
     601  double pXmin = +2.0;
     602  double pXmax = -2.0;
     603  double pYmin = +2.0;
     604  double pYmax = -2.0;
     605
     606  *Tmean = 0.0;
     607  for (k = 0; k < Npoints; k++) {
     608    Tmin = MIN(Tmin, points[k].T);
     609    Tmax = MAX(Tmax, points[k].T);
     610    Tmean += points[k].T;
     611
     612    // at this point, T is in years since J2000
     613    ParFactor (&points[k].pX, &points[k].pY, points[k].R, points[k].D, points[k].T);
     614    pXmin = MIN (pXmin, points[k].pX);
     615    pXmax = MAX (pXmax, points[k].pX);
     616    pYmin = MIN (pYmin, points[k].pY);
     617    pYmax = MAX (pYmax, points[k].pY);
     618  }
     619  *Trange = Tmax - Tmin;
     620
     621  // mean epoch
     622  *Tmean /= (float) Npoints;
     623
     624  // for HIGH_SPEED, just use the center of the range
     625  if (RELASTRO_OP == OP_HIGH_SPEED) {
     626    *Tmean = 0.5*(Tmax - Tmin);
     627  }
     628
     629  double dXRange = pXmax - pXmin;
     630  double dYRange = pYmax - pYmin;
     631  *parRange = hypot (dXRange, dYRange);
     632
     633  /* we need to do the fit in a locally linear space; choose a ref coordinate */
     634  coords.crval1 = points[0].R;
     635  coords.crval2 = points[0].D;
     636
     637  // project all of the R,D coordinates to a plane centered on this coordinate. set
     638  // the times to be relative to Tmean
     639  for (k = 0; k < Npoints; k++) {
     640    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &coords);
     641    points[k].T -= Tmean;
     642  }       
     643  return TRUE;
     644}
     645
     646int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
     647
     648  int i, j;
     649
     650  Nmax = 0;
     651  for (i = 0; i < Ncatalog; i++) {
     652    for (j = 0; j < catalog[i].Naverage; j++) {
     653      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
     654    }
     655  }
     656
     657  return Nmax;
     658}
    699659
    700660/* fitting proper-motion and parallax:
Note: See TracChangeset for help on using the changeset viewer.