Changeset 32854

Timestamp:

Dec 2, 2011, 11:47:26 AM (15 years ago)

Author:

eugene

Message:

testing the parallax factor; use slightly more accurate solar ephemeris

Location:

branches/eam_branches/ipp-20111122/Ohana/src/relastro

Files:

• Makefile (modified) (2 diffs)
• doc/plx_factor.pro (added)
• src/ParFactor.c (modified) (5 diffs)
• src/UpdateObjects.c (modified) (1 diff)
• src/testparallax.c (modified) (1 diff)

branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile

@@ -18 +18 @@
 
 relastro: $(BIN)/relastro.$(ARCH)
-install: $(DESTBIN)/relastro
+install: $(DESTBIN)/relastro $(DESTBIN)/testparallax
 
 RELASTRO = \
@@ -68 +68 @@
 $(RELASTRO): $(INC)/relastro.h $(KAPA_INCS)
 $(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS)
+
+testparallax: $(BIN)/testparallax.$(ARCH)
+
+TESTPAR = \
+$(SRC)/FitPMandPar.$(ARCH).o         \
+$(SRC)/ParFactor.$(ARCH).o           \
+$(SRC)/fitpoly.$(ARCH).o             \
+$(SRC)/mkpolyterm.$(ARCH).o            \
+$(SRC)/testparallax.$(ARCH).o
+
+$(TESTPAR): $(INC)/relastro.h $(KAPA_INCS)
+$(BIN)/testparallax.$(ARCH): $(TESTPAR) $(KAPA_LIBS)

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ParFactor.c

@@ -1 +1 @@
 # include "relastro.h"
+# define J2000 2451545.       /* Julian date at standard epoch */
 
 # if (0)
@@ -31 +32 @@
   double n, L, g;
 
-# define J2000 2451545.       /* Julian date at standard epoch */
 
   n = jd - J2000;
@@ -46 +46 @@
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) {
 
-# define J2000 2451545.       /* Julian date at standard epoch */
-
   double n = jd - J2000;              // day number
   double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
@@ -60 +58 @@
 
 /* given RA, DEC, Time, calculate the parallax factor */
+// Time is relative to Tmean, Tmean is years relative to J2000
 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) {
 
@@ -66 +65 @@
   /* given a time T in UNIX seconds, determine the solar longitude S */
 
-  jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
+  // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
+  jd = 365.25*(Time + Tmean) + J2000;
+  // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd);
+
   sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius);
 

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

@@ -133 +133 @@
         R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
         D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+
+        // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects
+        // to get Time in years relative to UNIX Tzero (1970/01/01)
         T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
 

branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/testparallax.c

@@ -1 +1 @@
 # include "relastro.h"
+# define NPTS 64
+# define J2000 2451545.       /* Julian date at standard epoch */
+# define T1970 2440587.500000 /* JD at UNIX ref time */
 
 int main (int argc, char **argv) {
+   
+  int i, Npts;
+  char line[1024];
+  double R[NPTS], D[NPTS], Time[NPTS], JD[NPTS];
+  double X[NPTS], Y[NPTS], dX[NPTS], dY[NPTS], pX[NPTS], pY[NPTS];
+  double Tref[NPTS], Tjyrs[NPTS], TrefS, TjyrsS, TrefMean, TjyrsMean, Ro, Do;
+  Coords coords;
+  PMFit fitPAR;
 
-    if (argc != N) {
-        fprintf (stderr, "USAGE: testparallax (file.dat)\n");
-        exit (1);
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: testparallax (file.dat)\n");
+    exit (1);
+  }
+
+  // test parallax program
+  FILE *f = fopen (argv[1], "r");
+  if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); }
+
+  // XXX uncomment to skip first line scan_line(f, line);
+    
+  Npts = 0;
+  for (i = 0; scan_line(f, line) != EOF; i++) {
+    if (Npts == NPTS) {
+      fprintf (stderr, "too many point: use dynamic alloc\n");
+      exit (2);
     }
+    if (line[0] == '#') continue;
+    dparse (&R[Npts],    1, line);
+    dparse (&D[Npts],    2, line);
+    dparse (&JD[Npts],   3, line);
+    dparse (&Time[Npts], 3, line);
 
-    // test parallax program
-    FILE *f = fopen (argv[1], "r");
-    if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); }
+    dparse (&dX[Npts],   4, line);
+    dparse (&dY[Npts],   4, line);
 
-    char buffer[1024];
-    scan_line(f, line);
-    
-    while (scan_line(f, line) != EOF) {
-        dparse (&jd,   1, line);
-        dparse (&time, 2, line);
-        dparse (&ra,   3, line);
-        dparse (&dec,  4, line);
+    // trent's file
+    // dparse (&JD[Npts],   1, line);
+    // dparse (&Time[Npts], 2, line);
+    // dparse (&R[Npts],    3, line);
+    // dparse (&D[Npts],    4, line);
+    Npts ++;
+  }
 
-    }
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  coords.crval1 = 0;
+  coords.crval2 = 0;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "RA---SIN");
 
-    // project to local coords
+  // use one point as a local reference
+  coords.crval1 = R[0];
+  coords.crval2 = D[0];
 
-    // run fitter
-    FitPMandPar
+  TrefS = TjyrsS = 0.0;
 
+  // project to local coords
+  for (i = 0; i < Npts; i++) {
+    RD_to_XY (&X[i], &Y[i], R[i], D[i], &coords);
 
+    Tjyrs[i] = Time[i] / 365.25;
+    Tref[i] = Time[i] / 365.25;
+    // Tunix[i] = (JD[i] - T1970) / 365.25; // time relative to T1970 in years
+    // Tjyrs[i] = (JD[i] - J2000) / 365.25; // time relative to J2000 in years
+
+    TrefS   += Tref[i];
+    TjyrsS   += Tjyrs[i];
+
+    // nominal for PS1
+    // dX[i] = 0.020;
+    // dY[i] = 0.020;
+  }
+  TrefMean = TrefS / Npts;
+  TjyrsMean = TjyrsS / Npts;
+
+  for (i = 0; i < Npts; i++) {
+    Tref[i] -= TrefMean;
+    ParFactor (&pX[i], &pY[i], R[i], D[i], Tjyrs[i], 0.0);
+    fprintf (stderr, "%f %f : %f %f : %f %f : %f %f\n", R[i], D[i], X[i], Y[i], Tref[i], Tjyrs[i], pX[i], pY[i]);
+  }
+
+  // run fitter
+  // FitPMandPar (&fitPAR, X, dX, Y, dY, Tref, pX, pY, Npts, TRUE);
+  FitPMandPar (&fitPAR, X, dX, Y, dY, Tjyrs, pX, pY, Npts, TRUE);
+  XY_to_RD (&Ro, &Do, fitPAR.Ro, fitPAR.Do, &coords);
+
+  fprintf (stderr, "Rx   : %f\n", fitPAR.Ro);
+  fprintf (stderr, "Dx   : %f\n", fitPAR.Do);
+  fprintf (stderr, "Ro   : %f\n", Ro);
+  fprintf (stderr, "Do   : %f\n", Do);
+  fprintf (stderr, "dRo  : %f\n", fitPAR.dRo);
+  fprintf (stderr, "dDo  : %f\n", fitPAR.dDo);
+  fprintf (stderr, "uR   : %f\n", fitPAR.uR);
+  fprintf (stderr, "uD   : %f\n", fitPAR.uD);
+  fprintf (stderr, "duR  : %f\n", fitPAR.duR);
+  fprintf (stderr, "duD  : %f\n", fitPAR.duD);
+  fprintf (stderr, "p    : %f\n", fitPAR.p);
+  fprintf (stderr, "dp   : %f\n", fitPAR.dp);
+  fprintf (stderr, "uTot : %f\n", hypot(fitPAR.uR,fitPAR.uD));
+  fprintf (stderr, "PA   : %f\n", DEG_RAD*atan2(fitPAR.uR,fitPAR.uD) + 360.0);
+  
+  exit (0);
 }
+