IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12731


Ignore:
Timestamp:
Apr 3, 2007, 10:54:48 AM (19 years ago)
Author:
eugene
Message:

fixed proper motion measurement, updated options, etc

Location:
trunk/Ohana/src/relastro
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/include/relastro.h

    r12332 r12731  
    9595int    STAR_BAD;
    9696int    MEAS_BAD;
    97 int    STAR_TOOFEW;
     97int    POS_TOOFEW;
     98int    PM_TOOFEW;
     99double PM_DT_MIN;
    98100int    IMAGE_TOOFEW;
    99101double IMAGE_GOOD_FRACTION;
    … …  
    272274int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
    273275int ParFactor (double *pR, double *pD, double R, double D, time_t T);
    274 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts);
    275 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts);
     276int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
     277int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
    276278
    277279Image *getMosaicForImage (int N);

  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r12332 r12731  
    2525
    2626  GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
    27   GetConfig (config, "STAR_TOOFEW",            "%d",  0, &STAR_TOOFEW);
     27  GetConfig (config, "PM_DT_MIN",              "%lf",  0, &PM_DT_MIN);
     28  GetConfig (config, "PM_TOOFEW",              "%d",  0, &PM_TOOFEW);
     29  GetConfig (config, "POS_TOOFEW",             "%d",  0, &POS_TOOFEW);
    2830  GetConfig (config, "IMAGE_TOOFEW",           "%d",  0, &IMAGE_TOOFEW);
    2931  GetConfig (config, "IMAGE_GOOD_FRACTION",    "%lf", 0, &IMAGE_GOOD_FRACTION);
    … …  
    3133  GetConfig (config, "GSCFILE",                "%s",  0, GSCFILE);
    3234  GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
    33   ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
    34   ScanConfig (config, "CATFORMAT",              "%s",  0, CATFORMAT);
    35   ScanConfig (config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
     35  ScanConfig(config, "CATMODE",                "%s",  0, CATMODE);
     36  ScanConfig(config, "CATFORMAT",              "%s",  0, CATFORMAT);
     37  ScanConfig(config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
    3638
    3739  sprintf (ImageCat, "%s/Images.dat", CATDIR);
    … …  
    4446  }
    4547
    46   GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
    47 
    48   GetConfig (config, "RELPHOT_GRID_X",         "%d",  0, &RELPHOT_GRID_X);
    49   GetConfig (config, "RELPHOT_GRID_Y",         "%d",  0, &RELPHOT_GRID_Y);
    50   GetConfig (config, "RELPHOT_GRID_BINNING",   "%d",  0, &RELPHOT_GRID_BINNING);
    51   GetConfig (config, "CAMERA_CONFIG",          "%s",  0, CameraConfig);
     48  // GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
     49  // GetConfig (config, "RELPHOT_GRID_X",         "%d",  0, &RELPHOT_GRID_X);
     50  // GetConfig (config, "RELPHOT_GRID_Y",         "%d",  0, &RELPHOT_GRID_Y);
     51  // GetConfig (config, "RELPHOT_GRID_BINNING",   "%d",  0, &RELPHOT_GRID_BINNING);
     52  // GetConfig (config, "CAMERA_CONFIG",          "%s",  0, CameraConfig);
    5253
    5354  if (*CATMODE == 0) strcpy (CATMODE, "RAW");

  • trunk/Ohana/src/relastro/src/FitPM.c

    r12332 r12731  
    22
    33/* 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, time_t *T, int Npts) {
     4int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
    55
    66  int i;
    … …  
    88  double **A, **B;
    99  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
     10  double chisq, Xf, Yf;
    1011
    1112  /* do I need to do this as 2 2x2 matrix equations? */
    … …  
    6970  array_free (B, 4);
    7071
     72  // 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  }
     80  fit[0].Nfit = Npts;
     81
     82  // the reduced chisq is divided by (Ndof = 2*Npts - 4)
     83  fit[0].chisq = chisq / (2.0*Npts - 4.0);
    7184  return (TRUE);
    7285}

  • trunk/Ohana/src/relastro/src/FitPMandPar.c

    r12332 r12731  
    22
    33/* 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, time_t *T, double *pR, double *pD, int Npts) {
     4int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
    55
    66  int i;

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r12332 r12731  
    88static double *pX;
    99static double *pY;
    10 static time_t *T;
     10static double *T;
    1111static double *dT;
    1212
    … …  
    2424  ALLOCATE (R, double, MAX (1, Nmax));
    2525  ALLOCATE (D, double, MAX (1, Nmax));
    26   ALLOCATE (T, time_t, MAX (1, Nmax));
     26  ALLOCATE (T, double, MAX (1, Nmax));
    2727  ALLOCATE (X, double, MAX (1, Nmax));
    2828  ALLOCATE (Y, double, MAX (1, Nmax));
    … …  
    4444  Coords coords;
    4545  PMFit fit;
     46  time_t To;
     47  int mode, Nave, Npm, Npar, Nskip;
     48  double Tmin, Tmax;
    4649
    4750  initObjectData (catalog, Ncatalog);
    4851
     52  /* project coordinates to a plane centered on the object with units of arcsec */
    4953  coords.crval1 = 0;
    5054  coords.crval2 = 0;
    … …  
    5761  strcpy (coords.ctype, "RA---SIN");
    5862
     63  // use J2000 as a reference time
     64  To = date_to_sec ("2000/01/01");
     65  Nave = Npar = Npm = 0;
     66
    5967  for (i = 0; i < Ncatalog; i++) {
     68
     69    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure);
     70
     71    Nskip = 0;
    6072    for (j = 0; j < catalog[i].Naverage; j++) {
    6173
    6274      /* calculate the average value of R,D for a single star */
    63       if (catalog[i].average[j].code & STAR_BAD) continue; 
     75      if (catalog[i].average[j].code & STAR_BAD) {
     76        Nskip ++;
     77        continue; 
     78      }
    6479
    6580      N = 0;
    6681      m = catalog[i].average[j].offset;
     82
     83      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
     84      mode = FIT_MODE;
     85
    6786      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
    6887        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    … …  
    7089        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    7190        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    72         T[N] = catalog[i].measure[m].t;
     91        T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years
     92
     93        Tmin = MIN(Tmin, T[N]);
     94        Tmax = MAX(Tmax, T[N]);
    7395
    7496        /* the astrometric errors are not being carried yet (but should be!) */
    … …  
    81103      }
    82104
     105      if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE;
     106      if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
     107
    83108      // XXX This criterion needs to be better considered: adjust to match Ndof
    84       if (N < STAR_TOOFEW) { /* too few measurements */
     109      if (N < POS_TOOFEW) { /* too few measurements */
    85110        catalog[i].average[j].code |= ID_STAR_FEW;
    86111        continue;
    … …  
    95120      /* project all of the R,D coordinates to a plane centered on this coordinate */
    96121      for (k = 0; k < N; k++) {
    97           RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
     122        RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
     123        dX[k] =  dR[k];
     124        dY[k] =  dD[k];
     125        // fprintf (stderr, "%d %f %f %f  %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]);
    98126      }   
    99127
    100128      /* fit the model components as needed */
    101       switch (FIT_MODE) {
     129      switch (mode) {
    102130        case FIT_AVERAGE:
    103131          liststats (R, dR, N, &statsR);
    … …  
    105133
    106134          fit.Ro = statsR.mean;
    107           fit.dRo = statsR.sigma;
     135          fit.dRo = 3600.0*statsR.sigma;
    108136
    109137          fit.Do = statsD.mean;
    110           fit.dDo = statsD.sigma;
     138          fit.dDo = 3600.0*statsD.sigma;
    111139
    112140          fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
    113141          fit.Nfit = N;
     142          Nave ++;
    114143          break;
    115144
    116145        case FIT_PM_ONLY:
    117146          FitPM (&fit, X, dX, Y, dY, T, N);
     147          // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
     148          // project Ro, Do back to RA,DEC
     149          XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
     150          // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
     151          // continue;
     152          Npm ++;
    118153          break;
    119154
    … …  
    123158          }
    124159          FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
     160          Npar ++;
    125161          break;
    126162        default:
    … …  
    130166
    131167      if (0 && (j < 100)) {
    132           fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
    133                    catalog[i].average[j].R,
    134                    catalog[i].average[j].D,
    135                    fit.Ro, fit.Do,
    136                    3600*(catalog[i].average[j].R - fit.Ro),
    137                    3600*(catalog[i].average[j].D - fit.Do));
     168        fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
     169                 catalog[i].average[j].R,
     170                 catalog[i].average[j].D,
     171                 fit.Ro, fit.Do,
     172                 3600*(catalog[i].average[j].R - fit.Ro),
     173                 3600*(catalog[i].average[j].D - fit.Do));
    138174      }
    139175
    … …  
    146182      }     
    147183
    148       catalog[i].average[j].R   = fit.Ro;
    149       catalog[i].average[j].D   = fit.Do;
    150       catalog[i].average[j].dR  = fit.dRo;
    151       catalog[i].average[j].dD  = fit.dDo;
    152 
    153       catalog[i].average[j].uR  = fit.uR;
    154       catalog[i].average[j].uD  = fit.uD;
    155       catalog[i].average[j].duR = fit.duR;
    156       catalog[i].average[j].duD = fit.duD;
    157 
    158       catalog[i].average[j].P   = fit.p;
    159       catalog[i].average[j].dP  = fit.dp;
     184      catalog[i].average[j].R   = fit.Ro; // RA in degrees
     185      catalog[i].average[j].D   = fit.Do; // DEC in degrees
     186      catalog[i].average[j].dR  = fit.dRo; // RA scatter in arcsec
     187      catalog[i].average[j].dD  = fit.dDo; // DEC scatter in arcsec
     188
     189      catalog[i].average[j].uR  = fit.uR; // RA proper motion in arcsec/year
     190      catalog[i].average[j].uD  = fit.uD; // DEC proper motion in arcsec/year
     191      catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year
     192      catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year
     193
     194      catalog[i].average[j].P   = fit.p; // parallax in arcsec
     195      catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec
    160196
    161197      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
    162 
    163198    }
     199
     200    if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip);
    164201  }
     202
     203  if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar);
    165204  return (TRUE);
    166205}
    … …  
    168207/* fitting proper-motion and parallax:
    169208
    170  given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
    171 
    172  let x,y be the coordinate in the linearized frame with y parallel to DEC lines
    173 
    174  L,B are the ecliptic longitude and latitude of the object,
    175  dL and dB are the offsets in the L and B directions
    176 
    177  dL = sin(t - topp)
    178  dB = cos(t - topp)*sin(B)
    179 
    180  these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
    181  for the source in the x,y frame is:
    182 
    183  x = Ro + uR * (t - to) + p * pR
    184  y = Do + uD * (t - to) + p * pD
    185 
    186  the unknowns in these equations are Ro, uR, Do, uD, and p
    187 
    188  XXX think through the concepts for the pole a bit better.  all objects near the pole
    189  move the same way with the same phase.  choose a projection center and define dL,dB relative
    190  to that center point coordinate system?
     209given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
     210
     211let x,y be the coordinate in the linearized frame with y parallel to DEC lines
     212
     213L,B are the ecliptic longitude and latitude of the object,
     214dL and dB are the offsets in the L and B directions
     215
     216dL = sin(t - topp)
     217dB = cos(t - topp)*sin(B)
     218
     219these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
     220for the source in the x,y frame is:
     221
     222x = Ro + uR * (t - to) + p * pR
     223y = Do + uD * (t - to) + p * pD
     224
     225the unknowns in these equations are Ro, uR, Do, uD, and p
     226
     227XXX think through the concepts for the pole a bit better.  all objects near the pole
     228move the same way with the same phase.  choose a projection center and define dL,dB relative
     229to that center point coordinate system?
    191230
    192231*/

  • trunk/Ohana/src/relastro/src/initialize.c

    r12332 r12731  
    4848    }
    4949    fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
    50     fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n",
    51              RELPHOT_GRID_X,
    52              RELPHOT_GRID_Y,
    53              RELPHOT_GRID_BINNING,
    54              (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
    5550
    5651    fprintf (stderr, "MAG_LIM                %lf\n", MAG_LIM);

  • trunk/Ohana/src/relastro/src/load_catalogs.c

    r12332 r12731  
    33Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
    44
    5   int i, Nstar;
     5  int i, j, k, m, Nstar;
    66  Catalog *catalog, *pcatalog, tcatalog;
    77
    … …  
    2020    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    2121    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
     22    pcatalog[0].Nsecfilt  = GetPhotcodeNsecfilt ();
    2223
    2324    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
    … …  
    3233      dvo_catalog_unlock (&tcatalog);
    3334      dvo_catalog_free (&tcatalog);
     35    } else {
     36      if (RESET) {
     37        for (j = 0; j < catalog[i].Naverage; j++) {
     38          catalog[i].average[j].code = 0;
     39          m = catalog[i].average[j].offset;
     40          for (k = 0; k < catalog[i].average[j].Nm; k++) {
     41            catalog[i].measure[m+k].flags = 0;
     42          }
     43        }
     44      }
    3445    }
    3546  }

  • trunk/Ohana/src/relastro/src/save_catalogs.c

    r12332 r12731  
    88  for (i = 0; i < Ncatalog; i++) {
    99
     10    if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename);
    1011    dvo_catalog_save (&catalog[i], VERBOSE);
    1112    dvo_catalog_unlock (&catalog[i]);

  • trunk/Ohana/src/relastro/src/select_images.c

    r12332 r12731  
    6666
    6767  if (VERBOSE) fprintf (stderr, "finding images\n");
     68  BuildChipMatch (timage, Ntimage);
    6869
    6970  nimage = 0;
    … …  
    9091    }
    9192   
     93    if (!FindMosaicForImage (timage, Ntimage, i)) continue;
     94
    9295    /* define image corners */
    9396    Xi[0] = 0;            Yi[0] = 0;
Note: See TracChangeset for help on using the changeset viewer.