Changeset 8200

Timestamp:

Aug 6, 2006, 11:51:06 AM (20 years ago)

Author:

eugene

Message:

added euler angle rotations

File:

• trunk/Ohana/src/opihi/cmd.astro/drizzle.c (modified) (12 diffs)

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

@@ -5 +5 @@
 static double XO, XX, XY;
 static double YO, YX, YY;
-int ZERO;
+static int ZERO;
+
+static int ROTATE;
+static double rot_phi, rot_alpha, rot_delta;
+static double rot_cdp, rot_sdp;
 
 int map_output_to_input (int Npix, double df);
@@ -13 +17 @@
 
 Coords coords_in, coords_out;
-Buffer *in, *out, *wt;
+Buffer *in, *out, *wt, *mask;
 
 int drizzle (int argc, char **argv) {
@@ -23 +27 @@
   if ((N = get_argument (argc, argv, "-zero"))) {
     ZERO = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  ROTATE = FALSE;
+  if ((N = get_argument (argc, argv, "-roll"))) {
+    /* -roll phi alpha_pole delta_pole */
+    /* XXX need to clarify the meaning of phi, alpha, delta */
+    ROTATE = TRUE;
+    remove_argument (N, &argc, argv);
+    rot_phi   = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    rot_alpha = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    rot_delta = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+
+    rot_cdp = cos(RAD_DEG*rot_delta);
+    rot_sdp = sin(RAD_DEG*rot_delta);
+  }
+
+  mask = NULL;
+  if ((N = get_argument (argc, argv, "-mask"))) {
+    remove_argument (N, &argc, argv);
+    if ((mask = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
   }
@@ -72 +100 @@
 
   int i, j, Ni, No, Nx, Ny, nx, ny;
-  float *Vin, *Vout, *Vwt;
+  float *Vin, *Vout, *Vwt, *Vmk;
   double x, y, X, Y;
 
@@ -80 +108 @@
   Vout = (float *) out[0].matrix.buffer;
   Vwt  = (float *) wt[0].matrix.buffer;
+  Vmk  = NULL;
+  Vmk = (mask == NULL) ? NULL : (float *) mask[0].matrix.buffer;
 
   nx = in[0].header.Naxis[0];
@@ -110 +140 @@
           Ni = (int)x + ((int)y)*nx;
 
+          if (Vmk && Vmk[Ni]) continue;
           Vout[No] += Vin[Ni];
           Vwt[No] ++;
@@ -123 +154 @@
 
   int i, j, Ni, No, Nx, Ny, nx, ny;
-  float *Vin, *Vout, *Vwt;
+  float *Vin, *Vout, *Vwt, *Vmk;
   double x, y, X, Y;
 
@@ -131 +162 @@
   Vout = (float *) out[0].matrix.buffer;
   Vwt  = (float *) wt[0].matrix.buffer;
+  Vmk  = NULL;
+  Vmk = (mask == NULL) ? NULL : (float *) mask[0].matrix.buffer;
 
   Nx = in[0].header.Naxis[0];
@@ -161 +194 @@
           No = (int)x + ((int)y)*nx;
 
+          if (Vmk && Vmk[Ni]) continue;
           Vout[No] += Vin[Ni];
           Vwt[No] ++;
@@ -170 +204 @@
 }
 
+int rotate_coords (double *phi, double *theta, double alpha, double delta) {
+
+  double sda, cda, cd, sd, sth;
+  double x, y;
+  
+  sda = sin(RAD_DEG*(alpha - rot_alpha));
+  cda = cos(RAD_DEG*(alpha - rot_alpha));
+  sd = sin(RAD_DEG*delta);
+  cd = cos(RAD_DEG*delta);
+  
+  sth = -cd*sda*rot_cdp + sd*rot_sdp;
+  y   = +cd*sda*rot_sdp + sd*rot_cdp;
+  x   = +cd*cda;
+
+  *theta = DEG_RAD*asin(sth);
+  *phi   = DEG_RAD*atan2(y,x) + rot_phi;
+  
+  while (*phi <   0.0) *phi += 360.0;
+  while (*phi > 360.0) *phi -= 360.0;
+  return (TRUE);
+}
+
 /* find the linear astrometric fix between images at this location */
 int set_linear_terms (Coords *in, Coords *out, int i, int j, int Npix) {
@@ -177 +233 @@
   double Xin, Yin, Xout, Yout;
   double Sx2, Sy2, Sxy, SXx, SXy, SYx, SYy;
-  double N, r, d;
+  double N, r, d, phi, theta;
 
   Xin = Yin = 0;
@@ -201 +257 @@
 
     if (!XY_to_RD (&r, &d, Xin, Yin, in)) return (FALSE);
+    if (ROTATE) { 
+      rotate_coords (&phi, &theta, r, d);
+      r = phi;
+      d = theta;
+    }
     if (!RD_to_XY (&Xout, &Yout, r, d, out)) return (FALSE);