Changeset 39591 for trunk/Ohana/src/opihi/cmd.astro/csystem.c

Timestamp:

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

Author:

eugene

Message:

add parallax factor function; add pm transformation

File:

• trunk/Ohana/src/opihi/cmd.astro/csystem.c (modified) (3 diffs)

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;
 }