IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20186


Ignore:
Timestamp:
Oct 16, 2008, 9:15:20 AM (18 years ago)
Author:
eugene
Message:

handle chip IDs which are not a simple sequence; add keywords to output grid fits image to identify camera and format; fix a bug in the mosaic grid coordinates; handle nan values in measurement mags; add stats on rejections; output plot and grid names based on outroot

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/src/GridOps.c

    r15509 r20186  
    2323  int *Fx, *Fy;  /* chip flip */
    2424  int *Ox, *Oy;  /* chip offset */
     25  int *valid;
    2526  char **ccdname;
    2627} camera;   
     28
     29static int *ccdnum;
     30static char *config;
    2731
    2832# ifdef GRID_V1
    … …  
    4953void initGrid (int dX, int dY) {
    5054
    51   int i;
    52   char field[64], *config, line[256];
     55  int i, N, ccdnum_max;
     56  char *p, field[64], line[256];
     57  int *Fx, *Fy;  /* chip flip */
     58  int *Ox, *Oy;  /* chip offset */
     59  char **ccdname;
    5360
    5461  /* load camera config file */
    … …  
    6673  ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
    6774
    68   ALLOCATE (camera.Fx, int, camera.Nchip);
    69   ALLOCATE (camera.Fy, int, camera.Nchip);
    70   ALLOCATE (camera.Ox, int, camera.Nchip);
    71   ALLOCATE (camera.Oy, int, camera.Nchip);
    72   ALLOCATE (camera.ccdname, char *, camera.Nchip);
     75  // temporary storage
     76  ALLOCATE (Fx,      int,    camera.Nchip);
     77  ALLOCATE (Fy,      int,    camera.Nchip);
     78  ALLOCATE (Ox,      int,    camera.Nchip);
     79  ALLOCATE (Oy,      int,    camera.Nchip);
     80  ALLOCATE (ccdname, char *, camera.Nchip);
     81  ALLOCATE (ccdnum,  int,    camera.Nchip);
    7382
    7483  /* load per-chip parameters */
    7584  for (i = 0; i < camera.Nchip; i++) {
    76     ALLOCATE (camera.ccdname[i], char, 256);
     85    ALLOCATE (ccdname[i], char, 256);
    7786    sprintf (field, "CCD.%d", i);
    7887    ScanConfig (config, field, "%s", 1, line);
    79     sscanf (line, "%s %d %d %d %d", camera.ccdname[i], &camera.Ox[i], &camera.Oy[i], &camera.Fx[i], &camera.Fy[i]);
    80   }
    81   free (config);
     88    // XXX get error status!
     89
     90    sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]);
     91   
     92    p = ccdname[i];
     93    while (!isdigit(*p) && *p) p++;
     94    if (*p == 0) continue;
     95    ccdnum[i] = atoi (p);
     96  }
     97
     98  /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */
     99  ccdnum_max = 0;
     100  for (i = 0; i < camera.Nchip; i++) {
     101    ccdnum_max = MAX(ccdnum_max, ccdnum[i]);
     102  }
     103  ccdnum_max ++;
     104
     105  if (ccdnum_max < camera.Nchip) {
     106    fprintf (stderr, "problem with camera config: duplicate ccd IDs\n");
     107    exit (1);
     108  }
     109 
     110  if (ccdnum_max > 0x1000) {
     111    fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max);
     112    exit (1);
     113  }
     114
     115  ALLOCATE (camera.Fx,      int,    ccdnum_max);
     116  ALLOCATE (camera.Fy,      int,    ccdnum_max);
     117  ALLOCATE (camera.Ox,      int,    ccdnum_max);
     118  ALLOCATE (camera.Oy,      int,    ccdnum_max);
     119  ALLOCATE (camera.valid,   int,    ccdnum_max);
     120  ALLOCATE (camera.ccdname, char *, ccdnum_max);
     121
     122  for (i = 0; i < ccdnum_max; i++) {
     123    camera.valid[i] = FALSE;
     124  }
     125
     126  for (i = 0; i < camera.Nchip; i++) {
     127    N = ccdnum[i];
     128    camera.Fx[N] = Fx[i];
     129    camera.Fy[N] = Fy[i];
     130    camera.Ox[N] = Ox[i];
     131    camera.Oy[N] = Oy[i];
     132    camera.ccdname[N] = ccdname[i];
     133    camera.valid[N] = TRUE;
     134  }
    82135
    83136  /* define mosaic 2d correction grid:
    … …  
    215268  bin[cat][meas] = i;
    216269  Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
    217   Ymeas[cat][meas] = Y + camera.Oy[ccdnum]*camera.Ny;
     270  Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
    218271  clist[i][Nlist[i]] = cat;
    219272  mlist[i][Nlist[i]] = meas;
    … …  
    248301void setMgrid (Catalog *catalog) {
    249302
    250   int i, j, m, c, n, N, Nmax;
     303  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys;
    251304  double *list, *dlist;
    252305  float Msys, Mrel, Mcal, Mmos;
    … …  
    261314  ALLOCATE (list, double, Nmax);
    262315  ALLOCATE (dlist, double, Nmax);
     316
     317  Nbad = Ncal = Nmos = Nrel = Nsys = 0;
    263318
    264319  for (i = 0; i < Ngrid; i++) {
    … …  
    270325      c = clist[i][j];
    271326     
    272       if (catalog[c].measure[m].dbFlags & MEAS_BAD) continue;
     327      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
     328        Nbad ++;
     329        continue;
     330      }
    273331      Mcal = getMcal  (m, c);
    274       if (isnan(Mcal)) continue;
     332      if (isnan(Mcal)) {
     333        Ncal ++;
     334        continue;
     335      }
    275336      Mmos = getMmos  (m, c);
    276       if (isnan(Mmos)) continue;
     337      if (isnan(Mmos)) {
     338        Nmos ++;
     339        continue;
     340      }
    277341      Mrel  = getMrel  (catalog, m, c);
    278       if (isnan(Mrel)) continue;
     342      if (isnan(Mrel)) {
     343        Nrel ++;
     344        continue;
     345      }
    279346     
    280347      n = catalog[c].measure[m].averef;
    281348      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
     349      if (isnan(Msys)) {
     350        Nsys++;
     351        continue;
     352      }
    282353      list[N] = Msys - Mrel - Mcal - Mmos;
    283354      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
    … …  
    290361    gridN[i] = N;
    291362  }
     363
     364  fprintf (stderr, "grid cells having too few measurements (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ncal, Nrel, Nsys);
     365
    292366  free (list);
    293367  free (dlist);
    … …  
    347421  graphdata.ymin = PlotdMmin;
    348422  graphdata.ymax = PlotdMmax;
    349   plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "XdM.png");
    350   plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
    351   plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
     423  plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "%s.XdM.png", OUTROOT);
     424  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT);
     425  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT);
    352426
    353427  plot_defaults (&graphdata);
    354   plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
     428  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT);
    355429
    356430  free (ylist);
    … …  
    363437void dump_grid () {
    364438
    365   int i, j, Nimage;
     439    int i, j, Nimage, Nbytes, Nformat;
    366440  int *imlist;
    367441  FILE *f;
    … …  
    369443  Matrix matrix;
    370444  Mosaic *refmosaic;
     445  char *filename;
     446  char formatline[32], key[32], value[64];
    371447   
     448  Nbytes = strlen (OUTROOT) + 6;
     449  ALLOCATE (filename, char, Nbytes);
     450  snprintf (filename, Nbytes, "%s.fits", OUTROOT);
     451
    372452  /* select reference mosaic image */
    373453  imlist = SelectRefMosaic (&refmosaic, &Nimage);
    374454
    375455  /* we are writing to this file */
    376   f = fopen ("mosaic.fits", "w");
     456  f = fopen (filename, "w");
    377457  if (f == (FILE *) NULL) {
    378     fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
     458    fprintf (stderr, "cannot open %s for output\n", filename);
     459    free (filename);
    379460    return;
    380461  }
    … …  
    388469  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
    389470  gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
     471
     472  // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system
     473  // Note that config must have been loaded (and not freed) above.
     474  ScanConfig (config, "NFORMAT", "%d", 1, &Nformat);
     475  for (i = 1; i <= Nformat; i++) {
     476      ScanConfig (config, "FORMAT", "%s", i, formatline);
     477      sscanf (formatline, "%s %s", key, value);
     478      gfits_modify (&header, key, "%s", 1, value);
     479  }
     480
    390481  gfits_fwrite_header (f, &header);
    391482  gfits_fwrite_matrix (f, &matrix);
    … …  
    483574  /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
    484575  for (i = 0; i < camera.Nchip; i++) {
    485     int ix, iy, x, y, X, Y, bin;
    486 
    487     gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[i]);
     576    int N, ix, iy, x, y, X, Y, bin;
     577
     578    N = ccdnum[i];
     579
     580    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
    488581    gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
    489582    gfits_modify (&theader, "NX", "%d", 1, camera.Nx);
    … …  
    501594        /* normalize X & Y */
    502595        x = X;
    503         if (camera.Fx[i]) x = RELPHOT_GRID_X - X - 1;
     596        if (camera.Fx[N]) x = RELPHOT_GRID_X - X - 1;
    504597        y = Y;
    505         if (camera.Fy[i]) y = RELPHOT_GRID_Y - Y - 1;
     598        if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1;
    506599             
    507600        /* coordinates in the grid */
    508         ix = x + camera.Ox[i]*RELPHOT_GRID_X;
    509         iy = y + camera.Oy[i]*RELPHOT_GRID_Y;
     601        ix = x + camera.Ox[N]*RELPHOT_GRID_X;
     602        iy = y + camera.Oy[N]*RELPHOT_GRID_Y;
    510603             
    511604        bin = ix + iy*gridX;
    … …  
    519612# endif
    520613
     614  free (filename);
    521615}
    522616
Note: See TracChangeset for help on using the changeset viewer.