Changeset 34826

Timestamp:

Dec 15, 2012, 7:16:10 AM (14 years ago)

Author:

eugene

Message:

add projection ops for boundary tree cells

Location:

branches/eam_branches/ipp-20121130/Ohana/src

Files:

• addstar/src/findskycell.c (modified) (3 diffs)
• libdvo/include/dvo.h (modified) (1 diff)
• libdvo/src/BoundaryTree.c (modified) (1 diff)

branches/eam_branches/ipp-20121130/Ohana/src/addstar/src/findskycell.c

@@ -202 +202 @@
     tree.ra[zone][band] = ra;
     tree.dec[zone][band] = dec;
+    tree.Xo[zone][band] = x;
+    tree.Yo[zone][band] = y;
     tree.cell[zone][band] = i;
 
@@ -342 +344 @@
     // the proj cell is divided into Nx, Ny bits.
     // (ra,dec) for (Ro,Do) -> (x,y).  given (Xo,Yo),(dX,dY) I can find ix,iy
-    // I currently track Ro,Do (tree->RA_origin[zone], tree->...)
-
+    // I currently track Ro,Do (tree->ra[zone][band], tree->dec[zone][band])
+    // I need to have the skycell center in pixels (Xo,Yo) and the scale dX,dY
+    
+    // convert R,D to X,Y with hard-wired projection and scale, orientation?
+    
+    TreeCellProjection (&x, &y, ra, dec, tree, zone, band);
+
+    xi = x / tree->NX_SUB;
+    yi = y / tree->NY_SUB;
+
+    
 
   }
@@ -349 +360 @@
   exit (0);
 }
+

branches/eam_branches/ipp-20121130/Ohana/src/libdvo/include/dvo.h

@@ -322 +322 @@
   int *NBAND;
 
-  double   **ra;
-  double  **dec;
-  int    **cell;
-  char  ***name;
+  double   **ra; // RA of projection cell center
+  double  **dec; // DEC of projection cell center
+  int    **cell; // zone,band -> proj cell sequence
+  char  ***name; // projection cell name
+  
+  int NX_SUB, NY_SUB;
+
+  int **Xo;
+  int **Yo;
 } BoundaryTree;
 

branches/eam_branches/ipp-20121130/Ohana/src/libdvo/src/BoundaryTree.c

@@ -318 +318 @@
 }
 
+
+// projection = TAN
+// need Ro, Do, Xo, Yo, dPix
+
+int TreeCellProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band) {
+
+    double Xo = tree->Yo[zone][band];
+    double Yo = tree->Xo[zone][band];
+    double Ro = tree->ra[zone][band];
+    double Do = tree->dec[zone][band];
+    double dPix = tree->dPix;
+
+    // this block only depends on Ro, Do
+
+    sdp  = sin(RAD_DEG*Do);
+    cdp  = cos(RAD_DEG*Do);
+    salp = sin(RAD_DEG*(ra - Ro));
+    calp = cos(RAD_DEG*(ra - Ro));
+    sdel = sin(RAD_DEG*dec);
+    cdel = cos(RAD_DEG*dec);
+    
+    stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
+    sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
+    cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
+
+    // defines the TAN projection (one of zenithal projections available, libdvo/src/coordops.c
+    // R = cot (theta) = cos(theta) / sin(theta)
+    if (stht == 0) {
+        Rc = hypot(sphi, cphi);
+        L = 180.0 * sphi / Rc;
+        M = 180.0 * cphi / Rc;
+    } else {
+        L = +DEG_RAD * sphi / stht;
+        M = -DEG_RAD * cphi / stht;
+    }
+
+    // scale, rotation, parity:
+    // rotation == 0.0 (pc1_1 == pc2_2 == 1.0, pc1_2 = pc2_1 = 0.0)
+
+    // if there were rotation or parity:
+    // Ro = (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+    // Xo = (coords[0].pc2_2*L - coords[0].pc1_2*M) / Ro;
+    // Yo = (coords[0].pc1_1*M - coords[0].pc2_1*L) / Ro;
+
+    Xo = L;
+    Yo = M;
+
+    // scale is dPix
+
+    *x = Xo / dPix + Xo;
+    *y = Yo / dPix + Yo;
+    
+    return TRUE;
+}
+