Changeset 38606

Timestamp:

Jul 20, 2015, 8:37:33 AM (11 years ago)

Author:

eugene

Message:

working on relastro

Location:

branches/eam_branches/ipp-20150625/Ohana/src

Files:

• addstar/src/loadgalphot_readstars.c (modified) (4 diffs)
• libautocode/def/average-ps1-v5.d (modified) (1 diff)
• libdvo/src/dvo_catalog.c (modified) (1 diff)
• libdvo/src/dvo_convert_PS1_V5.c (modified) (3 diffs)
• relastro/src/UpdateObjects.c (modified) (17 diffs)

branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_readstars.c

@@ -15 +15 @@
 
 static float *chisqFit = NULL;
+static float ZeroPt = 0;
 
 GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars) {
@@ -67 +68 @@
     namehash[i][2] = -1;
   }
+
+  ZeroPt = GetZeroPoint();
 
   char name[256], tmpname[256];
@@ -228 +231 @@
     stars[i].galphot.modelType = MODEL_TYPE_gal[i];
     stars[i].galphot.detID     = ID_gal[i];
+    stars[i].galphot.phtocode  = options->photcode;
 
     // I have a grid of measurements with (Flux, dFlux, Chisq) at each point
@@ -377 +381 @@
 
   // use bilinear interpolation to choose Flux @ (Xmin,Ymin)?
-  galphot->mag    = -2.5*log10(flux[iMin]); // correct for exptime?
+  galphot->mag    = -2.5*log10(flux[iMin]) + ZeroPt; // correct for exptime?
   galphot->magErr = sqrt(fluxErr[iMin]) / flux[iMin];
   galphot->chisq  = chisqMin;

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

@@ -69 +69 @@
 FIELD catID,          CAT_ID,      unsigned int,    unique ID for table in which object was first realized
 FIELD extID,          EXT_ID,      uint64_t,        external ID for object (eg PSPS objID)
-FIELD extIDgc,        EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
+
+# replace extIDgc (unused) with uRgal, uDgal:
+# FIELD extIDgc,      EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
+FIELD uRgal,          U_RA_GAL,    float,           modeled proper motion based on galactic motion
+FIELD uDgal,          U_DEC_GAL,   float,           modeled proper motion based on galactic motion
 
 # 4 double, 2 uint64_t, 17 float, 12 int, 6 short

branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_catalog.c

@@ -200 +200 @@
   average->catID           = 0;
   average->extID           = 0;
-  average->extIDgc         = 0;
+
+  average->uRgal           = NAN;
+  average->uDgal           = NAN;
 }
 

branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c

@@ -218 +218 @@
     out[i].catID         = in[i].catID;
     out[i].extID         = in[i].extID;
-    out[i].extIDgc       = in[i].extIDgc;
+
+    out[i].uRgal         = in[i].uRgal;
+    out[i].uDgal         = in[i].uDgal;
   }
   return (out);
@@ -286 +288 @@
     out[i].catID          = in[i].catID;
     out[i].extID          = in[i].extID;
-    out[i].extIDgc        = in[i].extIDgc;
+
+    out[i].uRgal         = in[i].uRgal;
+    out[i].uDgal         = in[i].uDgal;
   }
   return (out);
@@ -1318 +1322 @@
     out[i].catID          = in[i].catID;
     out[i].extID          = in[i].extID;
-    out[i].extIDgc        = in[i].extIDgc;
   }
   return (out);

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

@@ -42 +42 @@
       UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m);
     }
-    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);
-    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);
+    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);
+    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);
     FitStatsSum (fitStatsChips, sumStatsChips);
     FitStatsSum (fitStatsStack, sumStatsStack);
@@ -51 +51 @@
   FitStatsFree (fitStatsStack);
 
-  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);
-  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);
+  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);
+  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);
 
   FitStatsFree (sumStatsChips);
@@ -63 +63 @@
 // be called with just MeasureTiny set and Measure == NULL
 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff) {
+
+  int k;
 
   int setRefColor = areImagesMatched();
@@ -82 +84 @@
   // an object with no measurements is externally supplied
   if (average[0].Nmeasure == 0) return TRUE;
-
-  int NcBlue = 0;
-  int NcRed = 0;
 
   int mode = FIT_MODE; // start with the globally-defined fit mode
@@ -108 +107 @@
 
   double Tmean, Trange, parRange;
-  UpdateObjects_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange);
+  FitAstromPoints_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange);
 
   // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -131 +130 @@
     for (k = 0; k < fitStats->NfitAlloc; k++) {
       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
-      // XXX init fitStats->fit[k]? set fitStats->fit[k].getChisq?
-      if (!FitPM (&fitStats->fit[k], fitStats->fitPM, fitStats->sample, fitStats->Npoints)) continue;
+      if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
       fitStats->Nfit ++;
     }
@@ -143 +141 @@
 
     // project Ro, Do back to RA,DEC
-    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
+    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
     if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
     if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
@@ -175 +173 @@
     for (k = 0; k < fitStats->NfitAlloc; k++) {
       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
-      FitPMandPar (&fitStats->fit[k], fitStats->fitPar, fitStats->sample, fitStats->Npoints);
+      FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
       fitStats->Nfit ++;
     }
@@ -186 +184 @@
     FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
 
-    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &coords);
+    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
     average[0].flags |= ID_STAR_FIT_PAR;
     fitStats->Npar ++;
@@ -206 +204 @@
     for (k = 0; k < fitStats->NfitAlloc; k++) {
       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
-      FitPosPMfixed (&fitStats->fit[k], fitStats->fitPos, fitStats->sample, fitStats->Npoints);
+      FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
       fitStats->Nfit ++;
     }
@@ -214 +212 @@
 
     // project Ro, Do back to RA,DEC
-    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &coords);
+    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
     average[0].flags |= ID_STAR_FIT_AVE;
     fitStats->Nave ++;
@@ -286 +284 @@
 
   // what is the offset relative to the mean fit position?
-  coords.crval1 = average[0].R;
-  coords.crval2 = average[0].D;
-  if (isnan(coords.crval1)) {
+  fitStats->coords.crval1 = average[0].R;
+  fitStats->coords.crval2 = average[0].D;
+  if (isnan(fitStats->coords.crval1)) {
     return (FALSE);
   }
-  if (isnan(coords.crval2)) {
+  if (isnan(fitStats->coords.crval2)) {
     return (FALSE);
   }
 
   double dXoff, dYoff;
-  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
+  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &fitStats->coords);
   float dPos = hypot (dXoff, dYoff);
   if (dPos > MaxMeanOffset) {
     if (fitStats->Noffset < 100) {
-      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
     }
     fitStats->Noffset ++;
@@ -333 +331 @@
   average[0].ChiSqPar   = fitPar.chisq;
 
-  average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
+  average[0].Tmean      = (Tmean * 86400 * 365.25) + fitStats->T2000;
   average[0].Trange     = (Trange * 86400 * 365.25);
   average[0].Npos       = fit.Nfit;
@@ -456 +454 @@
 
   // what is the offset relative to the mean fit position?
-  coords.crval1 = average[0].R;
-  coords.crval2 = average[0].D;
+  fitStats->coords.crval1 = average[0].R;
+  fitStats->coords.crval2 = average[0].D;
 
   double dXoff, dYoff;
-  RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords);
+  RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &fitStats->coords);
   float dPos = hypot (dXoff, dYoff);
   if (dPos > MaxMeanOffset) {
     if (fitStats->Noffset < 100) {
-      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
     }
     fitStats->Noffset ++;
@@ -490 +488 @@
   for (k = 0; k < average[0].Nmeasure; k++) {
 
-    if (XVERB) {
+    if (0) {
       char *date = ohana_sec_to_date (measure[k].t);
       int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
@@ -520 +518 @@
     // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
     if (ExcludeBogus) {
-      coords.crval1 = average[0].R;
-      coords.crval2 = average[0].D;
+      fitStats->coords.crval1 = average[0].R;
+      fitStats->coords.crval2 = average[0].D;
       double Xi, Yi;
-      RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
+      RD_to_XY (&Xi, &Yi, Ri, Di, &fitStats->coords);
       double radius = hypot(Xi, Yi);
       if (radius > ExcludeBogusRadius) {
@@ -546 +544 @@
     // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
     // T[] is time in years since J2000 (jd = 2451545)
-    points[Npoints].T = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
+    points[Npoints].T = (measure[k].t - fit->T2000) / (86400*365.25) ; // time relative to J2000 in years
 
     // add measured systematic error in quadrature?  only do this after the fit has