Changeset 12047

Timestamp:

Feb 25, 2007, 4:22:22 PM (19 years ago)

Author:

eugene

Message:

adding pm / par fiting functions

Location:

branches/dvo-mods-2007-02/Ohana/src/relastro/src

Files:

• FitPM.c (added)
• UpdateObjects.c (modified) (6 diffs)
• dvo_astrom_ops.c (added)

branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c

@@ -2 +2 @@
 
 static int Nmax;
-static double *listR, *listD;
-static double *dlistR, *dlistD;
+static double *listX, *dlistX;
+static double *listY, *dlistY;
+static double *listR, *dlistR;
+static double *listD, *dlistD;
+static double *listT, *dlistT;
 
 void initObjectData (Catalog *catalog, int Ncatalog) {
@@ -25 +28 @@
   float chisq;
   StatType statsR, statsD;
+  Coords coords;
+
+  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");
 
   for (i = 0; i < Ncatalog; i++) {
@@ -39 +53 @@
         listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
         listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        listT[N] = catalog[i].measure[m].t;
 
         /* the astrometric errors are not being carried yet (but should be!) */
@@ -44 +59 @@
         dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
         dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+        dlistT[N] = 1.0; /* XXX hard-wire this for now */
+
         N++;
       }
@@ -52 +69 @@
       } 
 
+      /* we need to do the fit in a locally linear space; choose a ref coordinate */
+      coords.crval1 = listR[0];
+      coords.crval2 = listD[0];
+      
+      /* project all of the R,D coordinates to a plane centered on this coordinate */
+      /* calculate pR[i], pD[i] for each point */
+      for (k = 0; k < N; k++) {
+          RD_to_XY (&listX[k], &listY[k], listR[k], listD[k], &coords);
+          ParFactor (&pX[k], &pY[k], listR[k], listD[k]);
+      }   
+
       /* in here, we should be fitting the parallax and proper-motion components */
+      /* XXX if I fit R,D,T as above, I will be getting the wrong answers near the pole */
+      if (FitProperMotion) {
+          fit = FitPM (listX, dlistX, listY, dlistY, listT, pR, pD, N);
+      }
 
-      liststats (listR, dlistR, N, &statsR);
-      liststats (listR, dlistR, N, &statsD);
+      if (AverageCoords) {
+          liststats (listR, dlistR, N, &statsR);
+          liststats (listD, dlistD, N, &statsD);
 
-      catalog[i].average[j].R = statsR.mean;
-      catalog[i].average[j].dR = statsR.sigma;
+          catalog[i].average[j].R = statsR.mean;
+          catalog[i].average[j].dR = statsR.sigma;
 
-      catalog[i].average[j].D = statsD.mean;
-      catalog[i].average[j].dD = statsD.sigma;
+          catalog[i].average[j].D = statsD.mean;
+          catalog[i].average[j].dD = statsD.sigma;
 
-      chisq = 0.5*(statsR.chisq + statsD.chisq);
+          chisq = 0.5*(statsR.chisq + statsD.chisq);
+      }
+
       catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
-
     }
   }
@@ -71 +105 @@
 }
 
+/* fitting proper-motion and parallax:
+
+ given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
+
+ let x,y be the coordinate in the linearized frame with y parallel to DEC lines
+
+ L,B are the ecliptic longitude and latitude of the object, 
+ dL and dB are the offsets in the L and B directions
+
+ dL = sin(t - topp)
+ dB = cos(t - topp)*sin(B)
+
+ these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
+ for the source in the x,y frame is:
+
+ x = Ro + uR * (t - to) + p * pR 
+ y = Do + uD * (t - to) + p * pD
+
+ the unknowns in these equations are Ro, uR, Do, uD, and p
+
+ XXX think through the concepts for the pole a bit better.  all objects near the pole 
+ move the same way with the same phase.  choose a projection center and define dL,dB relative 
+ to that center point coordinate system?
+
+*/