Changeset 39018 for trunk

Timestamp:

Nov 2, 2015, 7:49:22 AM (11 years ago)

Author:

eugene

Message:

add -forcedgalaxy mode to mkcmf; fix fitting error in loadgalphot

Location:

trunk/Ohana/src/addstar/src

Files:

• loadgalphot_fit2d.c (modified) (3 diffs)
• mkcmf.c (modified) (7 diffs)

trunk/Ohana/src/addstar/src/loadgalphot_fit2d.c

@@ -183 +183 @@
       if (!finite(*x) || !finite(*y) || !finite(*z)) continue;
       *zf = 0;
-      double X = 1;
+      double Y = 1;
       for (i = ny = 0; ny < fit->nterm; ny++) {
-        double Y = X;
+        double X = Y;
         for (nx = 0; nx < fit->nterm - ny; nx++, i++) {
-          *zf += fit->b[i][0]*Y;
-          Y = Y * (*y);
-        }
-        X = X * (*x);
+          *zf += fit->b[i][0]*X;
+          X = X * (*x);
+        }
+        Y = Y * (*y);
       }
     }
@@ -213 +213 @@
     maxsigma = fit->ClipNsigma * sigma;
 
-    // if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma);
+    if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma);
 
     /* mask outlier points */
@@ -230 +230 @@
       } 
     }
-    // if (VERBOSE) gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask);
+    if (VERBOSE) {
+      gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask);
+    }
   }
     

trunk/Ohana/src/addstar/src/mkcmf.c

@@ -29 +29 @@
 int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
 int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, int Nrad);
+int WriteXGALtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
 int WriteDETFtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
 
@@ -62 +63 @@
   int isStack = FALSE;
   int isForcedWarp = FALSE;
+  int isForcedGalaxy = FALSE;
   int isDiff  = FALSE;
   int WriteXSRC = FALSE;
   int WriteXFIT = FALSE;
   int WriteXRAD = FALSE;
+  int WriteXGAL = FALSE;
   int WriteDETF = FALSE;
   if ((N = get_argument (argc, argv, "-stack"))) {
@@ -76 +79 @@
     WriteXSRC = WriteXFIT = WriteXRAD = TRUE;
     isForcedWarp = TRUE;
+  }
+  if ((N = get_argument (argc, argv, "-forcedgalaxy"))) {
+    remove_argument (N, &argc, argv);
+    WriteXSRC = WriteXFIT = WriteXRAD = TRUE;
+    WriteXGAL = TRUE;
+    isForcedGalaxy = TRUE;
   }
   if ((N = get_argument (argc, argv, "-diff"))) {
@@ -409 +418 @@
     gfits_modify (&header, "PHOT_V",                "%s", 1, "38100");
   }
+  if (isForcedGalaxy) {
+    gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, ""); // does this affect anything? I don't actually think so...
+    gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, 25.0);
+  }
   if (isDiff) {
     gfits_modify (&header, "HIERARCH FPA.ZP",  "%f", 1, zeroPt);
@@ -506 +519 @@
   if (WriteXRAD && isStack) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars, 3);
   if (WriteXRAD && isForcedWarp) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars, 1);
+  if (WriteXGAL) WriteXGALtable (fits, extroot, X, Y, M, Flag, Nstars);
 
   if (WriteDETF) WriteDETFtable (fits, extroot, X, Y, M, Flag, Nstars);
@@ -1541 +1555 @@
 }
 
+# define NMODEL 9
+static char ModelNames[NMODEL][16] = {
+  "PS_MODEL_GAUSS", 
+  "PS_MODEL_PGAUSS", 
+  "PS_MODEL_QGAUSS", 
+  "PS_MODEL_PS1V1", 
+  "PS_MODEL_RGAUSS", 
+  "PS_MODEL_SERSIC", 
+  "PS_MODEL_EXP", 
+  "PS_MODEL_DEV", 
+  "PS_MODEL_TRAIL"};
+
 int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
 
@@ -1677 +1703 @@
     gfits_set_bintable_column (&header, &ftable, "EXT_THETA_ERR"    , EXT_THETA_ERR    , Nstars);
     gfits_set_bintable_column (&header, &ftable, "EXT_PAR_07"       , EXT_PAR_07       , Nstars);
+
+    // add model type info in header:
+    gfits_modify (&header, "MTNUM", "%d", 1, NMODEL);
+    for (i = 0; i < NMODEL; i++) {
+      char field[64];
+      snprintf (field, 64, "MTNAM%02d", i);
+      gfits_modify (&header, field,   "%s", 1, ModelNames[i]);
+      snprintf (field, 64, "MTVAL%02d", i);
+      gfits_modify (&header, field,   "%d", 1, i);
+    }
+  }
+  /**/
+
+  gfits_fwrite_Theader (fits, &header);
+  gfits_fwrite_table  (fits, &ftable);
+  gfits_free_header (&header);
+  gfits_free_table (&ftable);
+
+  return TRUE;
+}
+
+# define NXGRID 5
+# define NYGRID 5
+# define BINSTEP 0.1
+
+int WriteXGALtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
+
+  int i;
+  Header header;
+  FTable ftable;
+
+  char extdata[80];
+  snprintf (extdata, 80, "%s.xgal", extroot);
+
+  /* Example Code to create a bintable */
+  {
+    gfits_create_table_header (&header, "BINTABLE", extdata);
+
+    int Nbin = NXGRID * NYGRID;
+    char fmt[16];
+    snprintf (fmt, 16, "%dE", Nbin);
+
+    // define the table layout
+    gfits_define_bintable_column (&header, "J",   "IPP_IDET"         , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned
+    gfits_define_bintable_column (&header, "J",   "MODEL_TYPE"       , "no comment", NULL, 1.0, 0.0); 
+    gfits_define_bintable_column (&header, "E",   "X_FIT"            , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "Y_FIT"            , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "NPIX"             , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, fmt,   "GAL_FLUX"         , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, fmt,   "GAL_FLUX_ERR"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, fmt,   "GAL_CHISQ"        , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MAJOR_MIN"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MAJOR_MAX"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MAJOR_DEL"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MINOR_MIN"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MINOR_MAX"     , "no comment", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&header, "E",   "FR_MINOR_DEL"     , "no comment", NULL, 1.0, 0.0);
+
+    // generate the output array that carries the data
+    gfits_create_table (&header, &ftable);
+
+    // create intermediate storage arrays
+    int    *IPP_IDET           ; ALLOCATE (IPP_IDET         ,  int   , Nstars);
+    int    *MODEL_TYPE         ; ALLOCATE (MODEL_TYPE       ,  int   , Nstars);
+    float  *X_FIT              ; ALLOCATE (X_FIT            ,  float , Nstars);
+    float  *Y_FIT              ; ALLOCATE (Y_FIT            ,  float , Nstars);
+    float  *NPIX               ; ALLOCATE (NPIX             ,  float , Nstars);
+    float  *GAL_FLUX           ; ALLOCATE (GAL_FLUX         ,  float , Nstars*Nbin);
+    float  *GAL_FLUX_ERR       ; ALLOCATE (GAL_FLUX_ERR     ,  float , Nstars*Nbin);
+    float  *GAL_CHISQ          ; ALLOCATE (GAL_CHISQ        ,  float , Nstars*Nbin);
+    float  *FR_MAJOR_MIN       ; ALLOCATE (FR_MAJOR_MIN     ,  float , Nstars);
+    float  *FR_MAJOR_MAX       ; ALLOCATE (FR_MAJOR_MAX     ,  float , Nstars);
+    float  *FR_MAJOR_DEL       ; ALLOCATE (FR_MAJOR_DEL     ,  float , Nstars);
+    float  *FR_MINOR_MIN       ; ALLOCATE (FR_MINOR_MIN     ,  float , Nstars);
+    float  *FR_MINOR_MAX       ; ALLOCATE (FR_MINOR_MAX     ,  float , Nstars);
+    float  *FR_MINOR_DEL       ; ALLOCATE (FR_MINOR_DEL     ,  float , Nstars);
+
+    // assign the storage arrays
+    for (i = 0; i < Nstars; i++) {
+      float flux = pow (10.0, -0.4*M[i]);
+
+      double ra, dec;
+      XY_to_RD (&ra, &dec, X[i], Y[i], &coords);
+
+      IPP_IDET         [i] = i;
+      MODEL_TYPE       [i] = 7;
+      X_FIT            [i] = X[i];
+      Y_FIT            [i] = Y[i];
+      NPIX             [i] = 10.0;
+      FR_MAJOR_MIN     [i] = FX*2.7;
+      FR_MAJOR_DEL     [i] = FX*BINSTEP;
+      FR_MAJOR_MAX     [i] = FR_MAJOR_MIN[i] + FR_MAJOR_DEL[i]*(NXGRID - 1); // this is probabaly an error in the real file, careful
+      FR_MINOR_MIN     [i] = FY*2.7; 
+      FR_MINOR_DEL     [i] = FY*BINSTEP; 
+      FR_MINOR_MAX     [i] = FR_MINOR_MIN[i] + FR_MINOR_DEL[i]*(NYGRID - 1); 
+
+      int xs = NXGRID / 2;
+      int ys = NYGRID / 2;
+
+      // is it major (fast) * minor (slow) or vice versa?
+      int Npt = 0;
+      int ix, iy;
+      for (ix = -xs; ix <= xs; ix++) {
+        for (iy = -ys; iy <= ys; iy++) {
+          float scale = (1.0 + ix*ix + 1.2*iy*iy); // a parabola with min at ix,iy = 0,0
+          GAL_FLUX     [Nbin*i + (ix + xs) + NXGRID*(iy + ys)] = flux / scale;
+          GAL_FLUX_ERR [Nbin*i + (ix + xs) + NXGRID*(iy + ys)] = sqrt(flux / scale);
+          GAL_CHISQ    [Nbin*i + (ix + xs) + NXGRID*(iy + ys)] = 10.0 * scale;
+          Npt ++;
+        }
+      }
+    }
+
+    // set the table data
+    gfits_set_bintable_column (&header, &ftable, "IPP_IDET"         , IPP_IDET         , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "MODEL_TYPE"       , MODEL_TYPE       , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "X_FIT"            , X_FIT            , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "Y_FIT"            , Y_FIT            , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "NPIX"             , NPIX             , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "GAL_FLUX"         , GAL_FLUX         , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "GAL_FLUX_ERR"     , GAL_FLUX_ERR     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "GAL_CHISQ"        , GAL_CHISQ         , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MAJOR_MIN"     , FR_MAJOR_MIN     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MAJOR_MAX"     , FR_MAJOR_MAX     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MAJOR_DEL"     , FR_MAJOR_DEL     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MINOR_MIN"     , FR_MINOR_MIN     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MINOR_MAX"     , FR_MINOR_MAX     , Nstars);
+    gfits_set_bintable_column (&header, &ftable, "FR_MINOR_DEL"     , FR_MINOR_DEL     , Nstars);
   }
   /**/