Changeset 39591

Timestamp:

May 28, 2016, 5:07:36 PM (10 years ago)

Author:

eugene

Message:

add parallax factor function; add pm transformation

Location:

trunk/Ohana/src/opihi/cmd.astro

Files:

• Makefile (modified) (1 diff)
• csystem.c (modified) (3 diffs)
• fitplx_irls.c (modified) (1 diff)
• init.c (modified) (2 diffs)
• parallax_factor.c (added)

trunk/Ohana/src/opihi/cmd.astro/Makefile

@@ -60 +60 @@
 $(SRC)/objload.$(ARCH).o           \
 $(SRC)/outline.$(ARCH).o           \
+$(SRC)/parallax_factor.$(ARCH).o   \
 $(SRC)/polar.$(ARCH).o             \
 $(SRC)/precess.$(ARCH).o           \

trunk/Ohana/src/opihi/cmd.astro/csystem.c

@@ -4 +4 @@
 
   /* USAGE: csystem [C/G/E/H] [C/G/E/H] [epoch] */
-  int i;
-  double X, Y, x, y;
-  opihi_flt *xptr, *yptr;
-  Vector *xvec, *yvec;
+  int i, N;
+  double X, Y, x, y, uX, uY, ux, uy;
   CoordTransformSystem input, output;
-  CoordTransform *transform;
+
+  Vector *xvec = NULL;
+  Vector *yvec = NULL;
+  Vector *uxvec = NULL;
+  Vector *uyvec = NULL;
+
+  // transform proper motions at the same time
+  char *uXname = NULL;
+  char *uYname = NULL;
+  int pmBackwards = FALSE;
+  if ((N = get_argument (argc, argv, "-pm"))) {
+    remove_argument (N, &argc, argv);
+    uXname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    uYname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  } 
+  if ((N = get_argument (argc, argv, "-pmback"))) {
+    if (uXname) {
+      gprint (GP_ERR, "cannot mix -pm and -pmback\n");
+      return FALSE;
+    }
+    remove_argument (N, &argc, argv);
+    uXname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    uYname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    pmBackwards = TRUE;
+  } 
 
   if (argc != 5) goto syntax;
@@ -35 +61 @@
   }
 
-  transform = InitTransform (input, output);
+  CoordTransform *transform = InitTransform (input, output);
   if (transform == NULL) {
     gprint (GP_ERR, "transform %c to %c is not yet defined\n", argv[1][0], argv[2][0]);
@@ -42 +68 @@
     
   if (SelectScalar (argv[3], &X)) {
-    if (!SelectScalar (argv[4], &Y)) return (FALSE);
+    if (!SelectScalar (argv[4], &Y)) goto escape;
       
+    // apply transform at the current coordinate to transform motions:
+    if (uXname) {
+      if (!SelectScalar (uXname, &uX)) goto escape;
+      if (!SelectScalar (uYname, &uY)) goto escape;
+
+      if (pmBackwards) {
+        TransformProperMotionBackwards (&ux, &uy, uX, uY, X, Y, transform);
+      } else {
+        TransformProperMotionForewards (&ux, &uy, uX, uY, X, Y, transform);
+      }
+    }
     ApplyTransform (&x, &y, X, Y, transform);
 
-    gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
+    if (uXname) {
+      gprint (GP_LOG, "%10.6f %10.6f @ %10.6f %10.6f\n", x, y, ux, uy);
+    } else {
+      gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
+    }
+    free (transform);
     return (TRUE);
   }
 
   /* find vectors */
-  if ((xvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((yvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((xvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((yvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) goto escape;
 
   if (xvec[0].Nelements != yvec[0].Nelements) {
     gprint (GP_ERR, "vectors %s and %s not the same length\n", argv[3], argv[4]);
-    return (FALSE);
+    goto escape;
   }
   
+  // apply transform at the current coordinate to transform motions:
+  if (uXname) {
+    if ((uxvec = SelectVector (uXname, OLDVECTOR, TRUE)) == NULL) goto escape;
+    if ((uyvec = SelectVector (uYname, OLDVECTOR, TRUE)) == NULL) goto escape;
+    
+    if (xvec[0].Nelements != uxvec[0].Nelements) {
+      gprint (GP_ERR, "vectors %s and %s not the same length\n", argv[3], uXname);
+      goto escape;
+    }
+    if (xvec[0].Nelements != uyvec[0].Nelements) {
+      gprint (GP_ERR, "vectors %s and %s not the same length\n", argv[3], uYname);
+      goto escape;
+    }
+
+    CastVector (uxvec, OPIHI_FLT);
+    CastVector (uyvec, OPIHI_FLT);
+  }
+
   CastVector (xvec, OPIHI_FLT);
   CastVector (yvec, OPIHI_FLT);
 
-  xptr = xvec[0].elements.Flt;
-  yptr = yvec[0].elements.Flt;
+  opihi_flt *xptr = xvec[0].elements.Flt;
+  opihi_flt *yptr = yvec[0].elements.Flt;
+
+  opihi_flt *uxptr = uxvec ? xvec[0].elements.Flt : NULL;
+  opihi_flt *uyptr = uyvec ? yvec[0].elements.Flt : NULL;
 
   for (i = 0; i < xvec[0].Nelements; i++, xptr++, yptr++) {
+    if (uXname) {
+      if (pmBackwards) {
+        TransformProperMotionBackwards (uxptr, uyptr, *uxptr, *uyptr, *xptr, *yptr, transform);
+      } else {
+        TransformProperMotionForewards (uxptr, uyptr, *uxptr, *uyptr, *xptr, *yptr, transform);
+      }
+      uxptr ++;
+      uyptr ++;
+    }
     ApplyTransform (xptr, yptr, *xptr, *yptr, transform);
   }
 
+  FREE (transform);
   return (TRUE);
 
+escape:
+  FREE (transform);
+  return FALSE;
+
 syntax:
-  gprint (GP_ERR, "USAGE: csystems [C/G/E/H] [C/G/E/H] X Y\n");
-  return (FALSE);
+  gprint (GP_ERR, "USAGE   : csystems [C/G/E/H] [C/G/E/H] X Y [-pm uX uY] [-pmback ux uy]\n");
+  gprint (GP_ERR, " -pm    : transform motions in source system to target system\n");
+  gprint (GP_ERR, " -pmback: transform motions in target system to source system\n");
+  gprint (GP_ERR, " e.g.   : csystem C G R D -pm uR uD -> R D will contain glon, glat; uR, uD will contain uL, uB\n");
+  gprint (GP_ERR, " e.g.   : csystem C G R D -pmback uL uB -> R D will contain glon, glat; uL, uB will contain uR, uD\n");
+  return FALSE;
 }

trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

@@ -123 +123 @@
   }
   // FitPMandPar_IRLS sets the values of Wx,Wy based on the fit distance
+
+  if (VERBOSE && Nresample) {
+    fprintf (stderr, "-- solution before resample --\n");
+    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f (%f, %f)\n", Rmean, Dmean, fit.dRo, fit.dDo, fit.Ro, fit.Do);
+    fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD);
+    fprintf (stderr, "par: %f +/- %f\n", fit.p, fit.dp);
+    fprintf (stderr, "chisq: %f Nfit %d\n", fit.chisq, fit.Nfit);
+  }
 
   // update the mask based on the input mask and the outlier limits.

trunk/Ohana/src/opihi/cmd.astro/init.c

@@ -48 +48 @@
 int objload                 PROTO((int, char **));
 int outline                 PROTO((int, char **));
+int parallax_factor         PROTO((int, char **));
 int polar                   PROTO((int, char **));
 int precess                 PROTO((int, char **));
@@ -117 +118 @@
   {1, "objload",     objload,      "plot obj data on Ximage "},
   {1, "outline",     outline,      "fit outline region"},
+  {1, "parallax_factor", parallax_factor, "generate parallax_factors"},
   {1, "polar",       polar,        "convert polar image to cartesian"},
   {1, "precess",     precess,      "precess coordinates"},