Changeset 23239 for branches/cnb_branches/cnb_branch_20090215/Ohana

Timestamp:

Mar 9, 2009, 3:41:26 PM (17 years ago)

Author:

beaumont

Message:

sync with trunk

Location:

branches/cnb_branches/cnb_branch_20090215

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/addstar/src/get2mass_full.c (modified) (10 diffs)
• Ohana/src/addstar/src/sky_tessalation.c (modified) (5 diffs)
• Ohana/src/imregister/base/get_fwhm.c (copied) (copied from trunk/Ohana/src/imregister/base/get_fwhm.c ) (1 prop)
• Ohana/src/imregister/detrend/usage.c (copied) (copied from trunk/Ohana/src/imregister/detrend/usage.c ) (1 prop)
• Ohana/src/imregister/imreg/load_probes.c (copied) (copied from trunk/Ohana/src/imregister/imreg/load_probes.c ) (1 prop)
• Ohana/src/libdvo (modified) (1 prop)
• Ohana/src/libdvo/lib (modified) (1 prop)
• Ohana/src/libfits/lib (modified) (1 prop)
• Ohana/src/libkapa/lib (modified) (1 prop)
• Ohana/src/libohana/lib (modified) (1 prop)
• Ohana/src/libtap/lib (modified) (1 prop)

branches/cnb_branches/cnb_branch_20090215/Ohana/src/addstar/src/get2mass_full.c

@@ -1 +1 @@
 # include "addstar.h"
 # include "2mass.h"
+
+// XXX check to see if desired output format is PS1_V1 or later?  (use 16bit version if not?)
 
 // fill in the data for a JHK triplet star.  takes a pointer to the start of the line the
@@ -171 +173 @@
 
   switch (qual) {
-    case 'X': star[0].measure.photFlags |= 0x0000; break;
-    case 'U': star[0].measure.photFlags |= 0x0001; break;
-    case 'F': star[0].measure.photFlags |= 0x0002; break;
-    case 'E': star[0].measure.photFlags |= 0x0003; break;
-    case 'A': star[0].measure.photFlags |= 0x0004; break;
-    case 'B': star[0].measure.photFlags |= 0x0005; break;
-    case 'C': star[0].measure.photFlags |= 0x0006; break;
-    case 'D': star[0].measure.photFlags |= 0x0007; break;
+    case 'A': star[0].measure.photFlags |= 0x00000001; break; // was: 0x0004
+    case 'B': star[0].measure.photFlags |= 0x00000002; break; // was: 0x0005
+    case 'C': star[0].measure.photFlags |= 0x00000004; break; // was: 0x0006
+    case 'D': star[0].measure.photFlags |= 0x00000008; break; // was: 0x0007
+    case 'E': star[0].measure.photFlags |= 0x00000010; break; // was: 0x0003
+    case 'F': star[0].measure.photFlags |= 0x00000020; break; // was: 0x0002
+    case 'U': star[0].measure.photFlags |= 0x00000040; break; // was: 0x0001
+    case 'X': star[0].measure.photFlags |= 0x00000080; break; // was: 0x0000
     default: 
       fprintf (stderr, "error!\n");
@@ -189 +191 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags |= 0x0000; break;
-    case '1': star[0].measure.photFlags |= 0x0010; break;
-    case '2': star[0].measure.photFlags |= 0x0020; break;
-    case '3': star[0].measure.photFlags |= 0x0030; break;
-    case '4': star[0].measure.photFlags |= 0x0040; break;
-    case '6': star[0].measure.photFlags |= 0x0050; break;
-    case '9': star[0].measure.photFlags |= 0x0060; break;
+    case '0': star[0].measure.photFlags |= 0x00000100; break; // was: 0x0000 
+    case '1': star[0].measure.photFlags |= 0x00000200; break; // was: 0x0010 
+    case '2': star[0].measure.photFlags |= 0x00000400; break; // was: 0x0020 
+    case '3': star[0].measure.photFlags |= 0x00000800; break; // was: 0x0030 
+    case '4': star[0].measure.photFlags |= 0x00001000; break; // was: 0x0040 
+    case '6': star[0].measure.photFlags |= 0x00002000; break; // was: 0x0050 
+    case '9': star[0].measure.photFlags |= 0x00004000; break; // was: 0x0060 
     default: 
       fprintf (stderr, "error!\n");
@@ -206 +208 @@
 
   switch (qual) {
-    case 'p': star[0].measure.photFlags |= 0x0000; break;
-    case 'c': star[0].measure.photFlags |= 0x0100; break;
-    case 'd': star[0].measure.photFlags |= 0x0200; break;
-    case 's': star[0].measure.photFlags |= 0x0300; break;
-    case 'b': star[0].measure.photFlags |= 0x0400; break;
-    case '0': star[0].measure.photFlags |= 0x0500; break;
+    case 'p': star[0].measure.photFlags |= 0x00010000; break; // was: 0x0000
+    case 'c': star[0].measure.photFlags |= 0x00020000; break; // was: 0x0100
+    case 'd': star[0].measure.photFlags |= 0x00040000; break; // was: 0x0200
+    case 's': star[0].measure.photFlags |= 0x00080000; break; // was: 0x0300
+    case 'b': star[0].measure.photFlags |= 0x00010000; break; // was: 0x0400
+    case '0': star[0].measure.photFlags |= 0x00020000; break; // was: 0x0500
     default: 
       fprintf (stderr, "error!\n");
@@ -222 +224 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags &= ~0x0008; break;
-    case '1': star[0].measure.photFlags &= ~0x0008; break;
-    default:  star[0].measure.photFlags |=  0x0008; break;
+    case '0': star[0].measure.photFlags &= ~0x00300000; break; // was: ~0x0008
+    case '1': star[0].measure.photFlags |=  0x00100000; break; // was: ~0x0008
+    default:  star[0].measure.photFlags |=  0x00200000; break; // was:  0x0008
   }      
   return (TRUE);
@@ -232 +234 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags &= ~0x0080; break;
-    case '1': star[0].measure.photFlags &= ~0x0080; break;
-    default:  
-      star[0].measure.photFlags |= 0x0080; 
+    case '0': star[0].measure.photFlags &= ~0x00c00000; break; // was: ~0x0080 
+    case '1': star[0].measure.photFlags |=  0x00400000; break; // was: ~0x0080 
+    default:  star[0].measure.photFlags |=  0x00800000;        // was:  0x0080
       star[0].measure.extNsigma = 100.0;
       break;
@@ -245 +246 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags &= ~0x0800; break;
-    case '1': star[0].measure.photFlags &= ~0x0800; break;
-    default:  star[0].measure.photFlags |=  0x0800; break;
+    case '0': star[0].measure.photFlags &= ~0x03000000; break; // was: ~0x0800
+    case '1': star[0].measure.photFlags |=  0x01000000; break; // was: ~0x0800
+    default:  star[0].measure.photFlags |=  0x02000000; break; // was:  0x0800
   }      
   return (TRUE);
@@ -255 +256 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags &= ~0x1000; break;
-    case '1': star[0].measure.photFlags &= ~0x1000; break;
-    default:  star[0].measure.photFlags |=  0x1000; break;
+    case '0': star[0].measure.photFlags &= ~0x0c000000; break; // was: ~0x1000
+    case '1': star[0].measure.photFlags |=  0x04000000; break; // was: ~0x1000
+    default:  star[0].measure.photFlags |=  0x08000000; break; // was:  0x1000
   }      
   return (TRUE);
@@ -265 +266 @@
 
   switch (qual) {
-    case '0': star[0].measure.photFlags &= ~0x2000; break;
-    case '1': star[0].measure.photFlags |=  0x2000; break;
+    case '0': star[0].measure.photFlags &= ~0x10000000; break; // was: ~0x2000
+    case '1': star[0].measure.photFlags |=  0x10000000; break; // was:  0x2000
     default:  abort();
   }      
@@ -272 +273 @@
 }
 
+// unused photFlags:
+// 0x0000.8000
+// 0x0004.0000
+// 0x0008.0000
+// 0x2000.0000
+// 0x4000.0000
+// 0x8000.0000

branches/cnb_branches/cnb_branch_20090215/Ohana/src/addstar/src/sky_tessalation.c

@@ -240 +240 @@
   // integer number of dec zones between -90 and +90
 
-  nDEC = 180.0 / CELLSIZE;
-  dDEC = 180.0 / nDEC;
+  // in fact, we place a single image on each pole, so the real range of dec is 180.0 - CELLSIZE:
+
+  nDEC = (180.0 - CELLSIZE) / CELLSIZE;
+  dDEC = (180.0 - CELLSIZE) / nDEC;
+  nDEC += 2;
+
+  // a test
+  // for (dec = 0.0 + 0.5*dDEC; dec < +90.0; dec += dDEC) {
 
   // generate the a collection of rectangles for each ring
-  for (dec = -90.0 + 0.5*dDEC; dec < +90.0; dec += dDEC) {
+  for (dec = -90.0; dec < +90.0 + 0.5*dDEC; dec += dDEC) {
 
     ring = sky_rectangle_ring (dec, dDEC, &Nring);
+    if (!ring) continue;
 
     // subdivide each image (Nx x Ny subcells)
@@ -529 +536 @@
 SkyRectangle *sky_rectangle_ring (float dec, float dDEC, int *nring) {
 
-  int i, nRA, NX, NY;
-  float dRA, decLower;
+  int i, NX, NY, nRA;
   SkyRectangle *ring;
+  float theta, dRA;
 
   // 'dec' is a guess at the center of the cell; in fact, we need to choose decLower and
@@ -537 +544 @@
 
   // we can determine the 'lower' bound (bound closest to the equator):
-  decLower = (dec > 0.0) ? dec - 0.5*dDEC : dec + 0.5*dDEC;
-
-  // Subdivide the 'lower' bound into an integer number of segments:
-  nRA = cos(dec*RAD_DEG) * 360.0 / CELLSIZE; // CELLSIZE is a projection size
-  dRA = 360.0 / nRA;                         // dRA is a size in RA degrees
+  float decLower = (dec > 0.0) ? dec - 0.5*dDEC : dec + 0.5*dDEC;
+
+  // solve for actual cellsize (\theta):  tan(\delta_{n+1} - \theta/2) = tan(\delta_n + \theta/2)cos(\alpha_n / 2)
+  float decUpper = (dec > 0.0) ? dec + dDEC : dec - dDEC;
+
+  if (fabs(dec) + 0.5*dDEC > 90.0) {
+    // onPole = TRUE;
+    theta = dDEC;
+    nRA = 1;
+    dRA = theta / cos(decLower*RAD_DEG); // make a square at the pole
+  } else {
+    // onPole = FALSE;
+    // Subdivide the 'lower' bound into an integer number of segments:
+    nRA = cos(RAD_DEG*decLower) * 360.0 / CELLSIZE; // CELLSIZE is a projection size
+    dRA = 360.0 / nRA;                         // dRA is a size in RA degrees == \alpha_n
+
+    // tan(decUpper - theta/2) = tan(dec + theta/2) cos(dRA / 2);
+
+    // we solve this equation for theta (fairly ugly: expand the tangents into sin/cos, expand the 
+    // sum-of-angle sine and cosine, multiply through, convert via half-angle formulae and write 
+    // as a quadratic expression in sine(theta/2)
+  
+    float sd1 = sin(RAD_DEG*decUpper);
+    float cd1 = cos(RAD_DEG*decUpper);
+    float sd2 = sin(RAD_DEG*dec);
+    float cd2 = cos(RAD_DEG*dec);
+    float   k = cos(RAD_DEG*dRA/2.0);
+
+    float c1 =  (sd1*cd2 + sd2*cd1)*(1.0 - k);
+    float c2 =  (sd1*cd2 - sd2*cd1)*(1.0 + k);
+    float c3 = -(sd1*sd2 + cd1*cd2)*(1.0 + k); 
+
+    float A = SQ(c3) + SQ(c2);
+    float B = 2*c1*c3;
+    float C = SQ(c1) - SQ(c2);
+
+    float arg = SQ(B) - 4.0*A*C;
+
+    float root;
+
+    if (dec >= 0.0) {
+      root = (-B + sqrt (arg)) / (2.0*A);
+      theta = +DEG_RAD*asin(root);
+    } else {
+      root = (-B - sqrt (arg)) / (2.0*A);
+      theta = -DEG_RAD*asin(root);
+    }
+
+    // the negative solution yields a negative cellsize 
+    // float root2 = (-B - sqrt (arg)) / (2.0*A);
+    // float theta2 = DEG_RAD*asin(root2);
+
+    // test lines:
+    // float r1 = tan(RAD_DEG*(decUpper - 0.5*theta1));
+    // float r2 = tan(RAD_DEG*(dec + 0.5*theta1));
+    // fprintf (stdout, "%f %f  %f  %f  %f %f  %f %f  %f %f %f\n", dec, decUpper, dRA, arg, root1, root2, theta1, theta2, r1, r2, k*r2);
+  }
+  fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
 
   // I think we need to return the value of dec for the next ring, but I am not sure...
@@ -559 +619 @@
   
     // range values are in projected degrees
-    NX = cos(dec*RAD_DEG) * dRA  * 3600.0 / SCALE;
-    NY =                    dDEC * 3600.0 / SCALE;
+    NX = cos(decLower*RAD_DEG) * dRA   * 3600.0 / SCALE;
+    NY =                         theta * 3600.0 / SCALE;
 
     // crpix1,crpix2 is the projection center
@@ -576 +636 @@
 
 
-    fprintf (stderr, "%f %f  : %f %f\n", 
-             ring[i].coords.crval1, ring[i].coords.crval2, 
-             ring[i].coords.crpix1, ring[i].coords.crpix2);
+    // fprintf (stderr, "%f %f  : %f %f\n", 
+    // ring[i].coords.crval1, ring[i].coords.crval2, 
+    // ring[i].coords.crpix1, ring[i].coords.crpix2);
   }
 

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libdvo

@@ -1 @@
-lib
+

@@ -1 @@
-lib
+

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libdvo/lib

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libfits/lib

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libkapa/lib

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libohana/lib

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

branches/cnb_branches/cnb_branch_20090215/Ohana/src/libtap/lib

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-lib
+

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so

@@ -1 @@
-*.linux
-*.lin64
-*.sol
-*.sun
-*.sid
-*.hp
-*.irix
-*.darwin.dylib
-*.darwin_x86.dylib
+*.a
+*.so