IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41603


Ignore:
Timestamp:
May 21, 2021, 10:03:06 AM (5 years ago)
Author:
eugene
Message:

major revision of the photometric flat-field correction / grid analysis; strip out code which uses astrometric information to define the grid

Location:
branches/eam_branches/relphot.20210521
Files:
21 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/Makefile

    r41556 r41603  
    2929$(SRC)/StatDataSetOps.$(ARCH).o  \
    3030$(SRC)/ConfigInit.$(ARCH).o      \
    31 $(SRC)/GridOps.$(ARCH).o         \
     31$(SRC)/GridOpsNew.$(ARCH).o      \
    3232$(SRC)/ImageOps.$(ARCH).o        \
    3333$(SRC)/ImageSubset.$(ARCH).o     \
    … …  
    9292$(SRC)/StatDataSetOps.$(ARCH).o  \
    9393$(SRC)/ConfigInit.$(ARCH).o      \
    94 $(SRC)/GridOps.$(ARCH).o         \
     94$(SRC)/GridOpsNew.$(ARCH).o      \
    9595$(SRC)/ImageOps.$(ARCH).o        \
    9696$(SRC)/ImageSubset.$(ARCH).o     \

  • branches/eam_branches/relphot.20210521/include/relphot.h

    r41557 r41603  
    44# include <signal.h>
    55# include <pthread.h>
    6 
    7 /* # define GRID_V1 */
    8 # define GRID_V2
    9 # define NO_IMAGE -100
    106
    117# define ID_SECF_STACK_PRIMARY 0x00004000
    … …  
    9995  char skipCal;               // if TRUE, this mosaic is incomplete and should not be calibrated
    10096  char inTGroup;
    101   Coords coords;
     97  Coords coords;              // this is only used to set the mosaic center for assign_images used by region hosts
    10298} Mosaic;
    10399
    … …  
    413409int    MOSAIC_ZEROPT;
    414410int    TGROUP_ZEROPT;
     411int    GRID_ZEROPT;
    415412
    416413int    FREEZE_IMAGES;
    … …  
    424421int    MANUAL_ITERATION;
    425422int    NLOOP;
    426 int    NGRID;
    427423int    RESET;
    428424int    RESET_ZEROPTS;
    … …  
    444440int    IMAGE_BAD;
    445441int    CALIBRATE_STACKS_AND_WARPS;
    446 int    USE_GRID;
     442
    447443int    KEEP_UBERCAL;
    448444char  *OUTROOT;
    … …  
    457453char  *PhotcodeList;
    458454
    459 int    RELPHOT_GRID_X;
    460 int    RELPHOT_GRID_Y;
    461 int    RELPHOT_GRID_BINNING;
    462 
    463455int      *photseclist;
    464456int      Nphotcodes;
    … …  
    489481SkyRegion UserPatch;
    490482char *UserCatalog;
     483
     484enum {GRID_ZPT_MODE_NONE, GRID_ZPT_MODE_ALL, };
     485int GRID_ZPT_MODE;
    491486
    492487enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, };
    … …  
    549544int SetSignals (void);
    550545
    551 void          freeGridBins        PROTO((int Ncatalog));
     546void          freeGridBins        PROTO((void));
    552547void          freeImageBins       PROTO((int Ncatalog, int doImageList));
    553548void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
    … …  
    561556float         getMcal             PROTO((off_t meas, int cat, dvoMagClassType class));
    562557float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
    563 float         getMgrid            PROTO((off_t meas, int cat));
     558float         getMgridTiny        PROTO((MeasureTiny *measure));
     559float         getMgrid            PROTO((Measure *measure));
    564560float         getMmos             PROTO((off_t meas, int cat));
    565561float         getMgrp             PROTO((off_t meas, int cat, float airmass, float *dZpt));
    … …  
    571567ImageSubset  *getimages_subset    PROTO((off_t *N));
    572568void          global_stats        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop));
    573 void          initGrid            PROTO((int dX, int dY));
    574 void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
     569void          initGrid            PROTO(());
     570void          initGridBins        PROTO(());
    575571void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
    576572void          initImagesSubset    PROTO((ImageSubset *input, off_t *line_number, off_t N));
    577573void          initImages          PROTO((Image *input, off_t *LineNumber, off_t N));
    578574void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    579 void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
     575void          initMosaicMcal      PROTO((Image *image, off_t Nimage));
    580576void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
    581577void          initTGroups         PROTO((Image *subset, off_t Nsubset));
    … …  
    621617void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
    622618void          plot_defaults       PROTO((Graphdata *graphdata));
    623 void          plot_grid           PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    624619void          plot_images         PROTO((void));
    625620void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) );
    … …  
    639634int           reload_images       PROTO((FITS_DB *db));
    640635int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
    641 void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
     636void          setMgrid            PROTO((Catalog *catalog, int Ncatalog));
     637void          resetMgrid          PROTO((void));
     638
    642639void          setMcalFromMosaics  PROTO((void));
    643640void          setMcalFromTGroups  PROTO((void));

  • branches/eam_branches/relphot.20210521/src/ConfigInit.c

    r41453 r41603  
    6363
    6464  GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
    65 
    66   GetConfig (config, "RELPHOT_GRID_X",         "%d",  0, &RELPHOT_GRID_X);
    67   GetConfig (config, "RELPHOT_GRID_Y",         "%d",  0, &RELPHOT_GRID_Y);
    68   GetConfig (config, "RELPHOT_GRID_BINNING",   "%d",  0, &RELPHOT_GRID_BINNING);
    6965  GetConfig (config, "CAMERA_CONFIG",          "%s",  0, CameraConfig);
    7066  GetConfig (config, "MOSAICNAME",             "%s",  0, MOSAICNAME);

  • branches/eam_branches/relphot.20210521/src/GridOps.c

    r41462 r41603  
    11# include "relphot.h"
    22
    3 enum {
    4     GRID_FITTED,
    5     GRID_FROZEN,
    6     GRID_REFERENCE,
    7 };
    8 
    9 static int     Ngrid;  // number of grid elements (gridX * gridY)
    10 static float   *gridM; // magnitude offset for this grid cell
    11 static float   *gridS; // stdev of the magnitude offset for this grid cell
    12 static int     *gridN; // number of stars used to measure the magnitude offset for this grid cell
    13 static int     *gridV; // data mode for this cell: fitted, frozen, reference
    14 static int      gridX; // number of grid elements in X direction
    15 static int      gridY; // number of grid elements in Y direction
    16 
    17 static int    **bin;   // link from catalog, measure to grid element
    18 static int    **Xmeas; // grid x-coordinate for a measurement
    19 static int    **Ymeas; // grid y-coordinate for a measurement
    20 
    21 static int    **clist; // link from measurement on a cell to catalog containing measurement
    22 static off_t  **mlist; // link from measurement on a cell to measurement in a catalog
    23 static off_t   *Nlist; // number of measurements for each grid cell
    24 static off_t   *NLIST; // allocated number of measurements for each grid cell
    25 
    26 static struct {
    27   int Nchip;
    28   int Mx, My;  /* mosaic size in chips */
    29   int Nx, Ny;  /* chip size in pixels */
    30   int *Fx, *Fy;  /* chip flip */
    31   int *Ox, *Oy;  /* chip offset */
    32   int *valid;
    33   char **ccdname;
    34 } camera;   
    35 
    36 static int *ccdnum;
    37 static char *config;
    38 
    39 void initGrid (int dX, int dY) {
    40   OHANA_UNUSED_PARAM(dX);
    41   OHANA_UNUSED_PARAM(dY);
    42 
    43   int i, N, ccdnum_max, refX, refY, refBin;
    44   char *p, field[64], line[256];
    45   int *Fx, *Fy;  /* chip flip */
    46   int *Ox, *Oy;  /* chip offset */
    47   char **ccdname;
    48 
    49   /* load camera config file */
    50   config = LoadConfigFile (CameraConfig);
    51   if (config == (char *) NULL) {
    52     fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig);
    53     exit (1);
    54   }
    55 
    56   /* load basic mosaic parameters */
    57   ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip);
    58   ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx);
    59   ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My);
    60   ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
    61   ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
    62 
    63   ScanConfig (config, "REFCELL.X", "%d", 1, &refX);
    64   ScanConfig (config, "REFCELL.Y", "%d", 1, &refY);
    65 
    66   // temporary storage
    67   ALLOCATE (Fx,      int,    camera.Nchip);
    68   ALLOCATE (Fy,      int,    camera.Nchip);
    69   ALLOCATE (Ox,      int,    camera.Nchip);
    70   ALLOCATE (Oy,      int,    camera.Nchip);
    71   ALLOCATE (ccdname, char *, camera.Nchip);
    72   ALLOCATE (ccdnum,  int,    camera.Nchip);
    73 
    74   /* load per-chip parameters */
    75   for (i = 0; i < camera.Nchip; i++) {
    76     ALLOCATE (ccdname[i], char, 256);
    77     sprintf (field, "CCD.%d", i);
    78     ScanConfig (config, field, "%s", 1, line);
    79     // XXX get error status!
    80 
    81     sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]);
     3/*
     4   We define a 'grid correction', essentially the flat-field correction, as a correction per
     5   photcode subdivided into an NxN array.  The dimensions of the chips corresponding to a
     6   photcode will need to be added to the photcode table.  This means a schema update,
     7   which I detest.  For now (2021.05.16), I will hard-wire the GPC1 / GPC2 chip size and
     8   worry about HSC & Megacam in the future.
     9
     10   We have a collection of photcodes, with photcodeID limited by design to 64k (unsigned short).
     11   
     12   Thus we can generate an array of pointers to the grid correction structures and access them
     13   by photcode.
     14
     15
     16 */
     17
     18typedef struct {
     19  unsigned short photcode;
     20  float **Mgrid; // grid of average corrections
     21  float **dMgrid; // grid of stdev of corrections
     22    int **nMgrid; // grid of number of stars for corrections
     23  int Nx; // bin = int(Xchip * (Nx / NxChip))
     24  int Ny; // bin = int(Ychip * (Ny / NyChip))
     25  float dX; // bin = int(Xchip * dX), dX = Nx / NxChip
     26  float dY; // bin = int(Ychip * dY), dY = Ny / NyChip
     27  // NxChip, NyChip = 4900,4900 for now
     28} GridCorrectionType;
     29
     30static GridCorrectionType **GridCorr = NULL;
     31static int                 NGridCorr = 0;
     32
     33// PS1: # define NX_CHIP 4900
     34// PS1: # define NY_CHIP 4900
     35// PS1: # define NX_BIN 16
     36// PS1: # define NY_BIN 16
     37
     38# define NX_CHIP 1000
     39# define NY_CHIP 1000
     40# define NX_BIN 2
     41# define NY_BIN 2
     42
     43void initGridBins () {
     44
     45  // allocate the full possible range of GridCorrectionType pointers, Nphotcode
     46  // loop over all images to find actual existing photcodes
     47  // generate initial grid values for each existing photcode
     48
     49  if (!GRID_ZEROPT) return; // skip if we are ignoring the grid correction
     50
     51  if (GridCorr) return;
     52
     53  PhotCodeData *photcodes = GetPhotcodeTable();
     54  if (!photcodes) return;
     55
     56  // we have photcodes->Ncodes actually loaded.  loop over them and allocate
     57  // array only as large as the max photcodes->code[i].code
     58
     59  int maxCode = 0;
     60  for (int i = 0; i < photcodes->Ncode; i++) {
     61    maxCode = MAX(maxCode, photcodes->code[i].code);
     62    myAssert (maxCode < 0x10000, "oops");
     63  }
     64
     65  NGridCorr = maxCode + 1;
     66  ALLOCATE (GridCorr, GridCorrectionType *, NGridCorr);
     67  for (int i = 0; i < NGridCorr; i++) {
     68    GridCorr[i] = NULL;
     69  }
     70  return;
     71}
     72
     73void freeGridBins() {
     74
     75  if (!GridCorr) return;
     76
     77  for (int code = 0; code < NGridCorr; code++) {
     78    if (!GridCorr[code]) continue;
     79
     80    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
     81      FREE (GridCorr[code]-> Mgrid[ix]);
     82      FREE (GridCorr[code]->dMgrid[ix]);
     83      FREE (GridCorr[code]->nMgrid[ix]);
     84    }     
     85
     86    FREE (GridCorr[code]-> Mgrid);
     87    FREE (GridCorr[code]->dMgrid);
     88    FREE (GridCorr[code]->nMgrid);
     89
     90    FREE(GridCorr[code]);
     91  }
     92
     93  FREE (GridCorr);
     94  GridCorr = NULL;
     95
     96  return;
     97}
     98
     99/* for GPC1, we have 60 chips, 5 filters, 16x16 grid cells = 2MB of memory for this stuff
     100   if we go to 64x64 grid cells (~75 pixels), then it is still only 29MB */
     101void initGrid () {
     102
     103  if (!GRID_ZEROPT) return;
     104
     105  off_t Nimages = 0;
     106  Image *images = getimages (&Nimages, NULL);
     107
     108  for (int i = 0; i < Nimages; i++) {
     109    int code = images[0].photcode;
     110    myAssert (code >= 0, "oops");
     111    myAssert (code < NGridCorr, "oops");
    82112   
    83     p = ccdname[i];
    84     while (!isdigit(*p) && *p) p++;
    85     if (*p == 0) continue;
    86     ccdnum[i] = atoi (p);
    87   }
    88 
    89   /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */
    90   ccdnum_max = 0;
    91   for (i = 0; i < camera.Nchip; i++) {
    92     ccdnum_max = MAX(ccdnum_max, ccdnum[i]);
    93   }
    94   ccdnum_max ++;
    95 
    96   if (ccdnum_max < camera.Nchip) {
    97     fprintf (stderr, "problem with camera config: duplicate ccd IDs\n");
    98     exit (1);
    99   }
    100  
    101   if (ccdnum_max > 0x1000) {
    102     fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max);
    103     exit (1);
    104   }
    105 
    106   ALLOCATE (camera.Fx,      int,    ccdnum_max);
    107   ALLOCATE (camera.Fy,      int,    ccdnum_max);
    108   ALLOCATE (camera.Ox,      int,    ccdnum_max);
    109   ALLOCATE (camera.Oy,      int,    ccdnum_max);
    110   ALLOCATE (camera.valid,   int,    ccdnum_max);
    111   ALLOCATE (camera.ccdname, char *, ccdnum_max);
    112 
    113   for (i = 0; i < ccdnum_max; i++) {
    114     camera.valid[i] = FALSE;
    115   }
    116 
    117   for (i = 0; i < camera.Nchip; i++) {
    118     N = ccdnum[i];
    119     camera.Fx[N] = Fx[i];
    120     camera.Fy[N] = Fy[i];
    121     camera.Ox[N] = Ox[i];
    122     camera.Oy[N] = Oy[i];
    123     camera.ccdname[N] = ccdname[i];
    124     camera.valid[N] = TRUE;
    125   }
    126 
    127   /* define mosaic 2d correction grid:
    128    * GRID_X is the number of grid pixels per chip in the X direction */
    129   gridX = RELPHOT_GRID_X * camera.Mx;
    130   gridY = RELPHOT_GRID_Y * camera.My;
    131   Ngrid = gridX * gridY;
    132 
    133   ALLOCATE (gridM, float, Ngrid);
    134   ALLOCATE (gridS, float, Ngrid);
    135   ALLOCATE (gridN, int,   Ngrid);
    136   ALLOCATE (gridV, int,   Ngrid);
    137 
    138   // the grid bins may have one of three possible states: fitted, frozen, reference
    139   // set the initial values to indicate that the bins are frozen
    140   for (i = 0; i < Ngrid; i++) {
    141     gridM[i] = 0.0;
    142     gridS[i] = 0.0;
    143     gridN[i] = 0;
    144     gridV[i] = GRID_FROZEN;
    145   }
    146 
    147   // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this)
    148   refBin = refX + refY*gridX;
    149   gridV[refBin] = GRID_REFERENCE;
    150 }
    151 
    152 void initGridBins (Catalog *catalog, int Ncatalog) {
    153 
    154   int i, j;
    155 
    156   if (!USE_GRID) return;
    157 
    158   /* define cat,meas -> grid pointers */
    159   ALLOCATE (bin, int *, Ncatalog);
    160   ALLOCATE (Xmeas, int *, Ncatalog);
    161   ALLOCATE (Ymeas, int *, Ncatalog);
    162   for (i = 0; i < Ncatalog; i++) {
    163     ALLOCATE (bin[i],   int, MAX (catalog[i].Nmeasure, 1));
    164     ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
    165     ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
    166     for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
    167   }
    168 
    169   /* define grid -> cat,meas pointers */
    170   ALLOCATE (Nlist, off_t,   Ngrid);
    171   ALLOCATE (NLIST, off_t,   Ngrid);
    172   ALLOCATE (clist, int *,   Ngrid);
    173   ALLOCATE (mlist, off_t *, Ngrid);
    174 
    175   for (i = 0; i < Ngrid; i++) {
    176     Nlist[i] = 0;
    177     NLIST[i] = 100;
    178     ALLOCATE (clist[i], int,   NLIST[i]);
    179     ALLOCATE (mlist[i], off_t, NLIST[i]);
    180   }
    181 }
    182 
    183 void freeGridBins (int Ncatalog) {
    184 
    185   int i;
    186 
    187   if (!USE_GRID) return;
    188 
    189   /* define cat,meas -> grid pointers */
    190   for (i = 0; i < Ncatalog; i++) {
    191     free (bin[i]);
    192     free (Xmeas[i]);
    193     free (Ymeas[i]);
    194   }
    195   free (bin);
    196   free (Xmeas);
    197   free (Ymeas);
    198 
    199   /* define grid -> cat,meas pointers */
    200   for (i = 0; i < Ngrid; i++) {
    201     free (clist[i]);
    202     free (mlist[i]);
    203   }
    204   free (Nlist);
    205   free (NLIST);
    206   free (clist);
    207   free (mlist);
    208 }
    209 
    210 int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum) {
    211 
    212   int ix, iy, Cx, Cy, i;
    213   double x, y;
    214 
    215   /* X, Y are chip coords on chip ccdnum */
    216 
    217   /* normalize X & Y */
    218   x = X;
    219   if (camera.Fx[ccdnum]) x = camera.Nx - X;
    220   y = Y;
    221   if (camera.Fy[ccdnum]) y = camera.Ny - Y;
    222 
    223   /* grid coords on the chip */
    224   Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
    225   Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
    226  
    227   /* coordinates in the grid */
    228   ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
    229   iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
    230 
    231   i = ix + iy*gridX;
    232 
    233   bin[cat][meas] = i;
    234   Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
    235   Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
    236   clist[i][Nlist[i]] = cat;
    237   mlist[i][Nlist[i]] = meas;
    238 
    239   // for the moment, add up the total grid count
    240   gridN[i]++;
    241 
    242   Nlist[i] ++;
    243   if (Nlist[i] == NLIST[i]) {
    244     NLIST[i] += 100;
    245     REALLOCATE (clist[i], int,   NLIST[i]);
    246     REALLOCATE (mlist[i], off_t, NLIST[i]);
    247   }     
    248   return (TRUE);
    249 
    250   fprintf (stderr, "error: star out of grid\n");
    251   exit (1);
    252 }
    253 
    254 int showGridCount() {
    255 
    256     int ix, iy, i;
    257 
    258     for (iy = 0; iy < gridY; iy++) {
    259         for (ix = 0; ix < gridX; ix++) {
    260             i = ix + iy*gridX;
    261             fprintf (stderr, "%3d  ", gridN[i]);
    262         }
    263         fprintf (stderr, "\n");
     113    // valid photcodes values (code) are in range 1 <= code < 0x10000
     114    // photcode == 0 are e.g., PHU (mosaic) images, and should be ignored here
     115    if (!code) continue;
     116
     117    if (GridCorr[code]) continue; // already created this one
     118
     119    ALLOCATE(GridCorr[code], GridCorrectionType, 1);
     120    GridCorr[code]->photcode = code;
     121    GridCorr[code]->Nx = NX_BIN;
     122    GridCorr[code]->Ny = NY_BIN;
     123    GridCorr[code]->dX = NX_BIN / (float) NX_CHIP;
     124    GridCorr[code]->dY = NY_BIN / (float) NX_CHIP;
     125
     126    // we are normally accessing this array randomly, so there is no advantage to
     127    // doing Mgrid[y][x] vs Mgrid[x][y]
     128    ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN);
     129    ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN);
     130    ALLOCATE (GridCorr[code]->nMgrid,   int *, NX_BIN);
     131    for (int ix = 0; ix < NX_BIN; ix++) {
     132      ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NX_BIN);
     133      ALLOCATE (GridCorr[code]->dMgrid[ix], float, NX_BIN);
     134      ALLOCATE (GridCorr[code]->nMgrid[ix],   int, NX_BIN);
     135    }     
     136  }
     137  resetMgrid(); // start with values of 0
     138}
     139
     140// reset the values in the arrays to 0
     141void resetMgrid () {
     142
     143  if (!GRID_ZEROPT) return;
     144  if (!GridCorr) return;
     145
     146  for (int code = 0; code < NGridCorr; code++) {
     147    if (!GridCorr[code]) continue;
     148
     149    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
     150      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
     151        GridCorr[code]-> Mgrid[ix][iy] = 0.0;
     152        GridCorr[code]->dMgrid[ix][iy] = 0.0;
     153        GridCorr[code]->nMgrid[ix][iy] =   0;
     154      }
     155    }     
     156  }
     157}
     158
     159void setMgrid (Catalog *catalog, int Ncatalog) {
     160
     161  // check if we are actually doing this step
     162  if (!GRID_ZEROPT) return;
     163  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return;
     164
     165  resetMgrid(); // start with values of 0
     166
     167  // loop over all measurements, accumulate Sum (in Mgrid), Sum2 (in dMgrid), and Npts
     168
     169  int Nsecfilt = GetPhotcodeNsecfilt ();
     170
     171  for (int nc = 0; nc < Ncatalog; nc++) {
     172    for (int na = 0; na < catalog[nc].Naverage; na++) {
     173
     174      int nm = catalog[nc].averageT[na].measureOffset;
     175      for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) {
     176
     177        // XXX I need to skip bad measurements, but I am not certain this
     178        // flag is set correctly.  see note below
     179        if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue;
     180
     181        float Mcal = getMcal  (nm, nc, MAG_CLASS_PSF);
     182        if (isnan(Mcal)) continue;
     183
     184        float Mgrp = getMgrp  (nm, nc, catalog[nc].measureT[nm].airmass, NULL);
     185        if (isnan(Mgrp)) continue;
     186
     187        float Mmos  = getMmos  (nm, nc);
     188        if (isnan(Mmos)) continue;
     189
     190        // Mrel* is the average magnitude for this star.  For PS1 stacks, we have too much
     191        // PSF variability.  We need to calibrate the PSF magnitudes separately from the
     192        // Aperture-like magnitues.  (We have an option to use the kron magnitudes or the
     193        // other apertures here).  I basically need to do this analysis separately for each
     194        // magnitude type
     195   
     196        float MrelPSF = getMrel  (catalog, nm, nc, MAG_CLASS_PSF, MAG_SRC_CHP);
     197        if (isnan(MrelPSF)) continue;
     198     
     199        float MsysPSF = PhotSysTiny (&catalog[nc].measureT[nm], &catalog[nc].averageT[na], &catalog[nc].secfilt[na*Nsecfilt], MAG_CLASS_PSF);
     200        if (isnan(MsysPSF)) continue;
     201
     202        // what about Mflat?
     203        float Moff =  Mcal + Mgrp + Mmos;
     204
     205        // Msys = Mrel + Moff + Mgrid
     206        // thus Mgrid = Msys - Mrel - Moff
     207
     208        int code = catalog[nc].measureT[nm].photcode;
     209        if (code <= 0) continue;
     210        if (code >= NGridCorr) continue; // does not match one of our image, skip
     211       
     212        GridCorrectionType *grid = GridCorr[code];
     213        if (!grid) continue; // does not match one of our images, skip
     214
     215        // edge effects could cause some positions to be slightly out of range
     216        // probably should trap extreme outliers
     217        int ix = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Xccd * grid->dX)), grid->Nx - 1);
     218        int iy = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Yccd * grid->dY)), grid->Ny - 1);
     219
     220        float dM = MsysPSF - MrelPSF - Moff;
     221       
     222        // XXX by the time we get here, we should have already mostly fixed up the zero
     223        // points.  reject measurements which are way off
     224        if (fabs(dM) > 0.5) continue;
     225
     226        grid-> Mgrid[ix][iy] += dM;
     227        grid->dMgrid[ix][iy] += dM*dM;
     228        grid->nMgrid[ix][iy] ++;
     229      }
    264230    }
    265     return (TRUE);
    266 }
    267 
    268 float getMgrid (off_t meas, int cat) {
    269 
    270   int i;
    271   float value;
    272 
    273   if (!USE_GRID) return (0);
    274   i = bin[cat][meas];
    275   if (i == -1) return (NAN);
    276 
    277   // during the grid annealing process, we skip over grid cells until they have enough
    278   // valid stars to be fitted
    279   if (gridV[i] == GRID_FROZEN) {
    280     return (NAN);
    281   }
    282 
    283   value = gridM[i];
    284   return (value);
    285 }
    286 
    287 /* direct (non-iterative) solution for Mgrid values for all grid bins */
    288 void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    289 
    290   int **gotstar, **gridmeas;
    291   int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp;
    292   double **A, **B, *Mjx, *Wjx;
    293   float Msys, Mcal, Mmos, Mgrp, Merr, Wsys;
    294   double Mj, Wj;
    295  
    296   if (!USE_GRID) return;
    297 
    298   ALLOCATE (A, double *, Ngrid);
    299   ALLOCATE (B, double *, Ngrid);
    300   for (i = 0; i < Ngrid; i++) {
    301     ALLOCATE (A[i], double, Ngrid);
    302     ALLOCATE (B[i], double, 1);
    303     memset (A[i], 0, Ngrid*sizeof(double));
    304     memset (B[i], 0, sizeof(double));
    305   }
    306 
    307   Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0;
    308 
    309   ALLOCATE (gotstar, int *, Ncatalog);
    310   for (i = 0; i < Ncatalog; i++) {
    311     ALLOCATE (gotstar[i], int, catalog[i].Naverage);
    312   } 
    313 
    314   // set up gridmeas table : grid index for each measurement
    315   ALLOCATE (gridmeas, int *, Ncatalog);
    316   for (i = 0; i < Ncatalog; i++) {
    317     ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure);
    318     for (j = 0; j < catalog[i].Nmeasure; j++) {
    319       gridmeas[i][j] = -1;
    320     }
    321   } 
    322   for (i = 0; i < Ngrid; i++) {
    323     for (j = 0; j < Nlist[i]; j++) {
    324       int m, c;
    325       m = mlist[i][j];
    326       c = clist[i][j];
    327       gridmeas[c][m] = i;
    328     }
    329   }
    330 
    331   // as we loop over the grid cells, we need to accumulate the values of Wjx for each star
    332   ALLOCATE (Wjx, double, Ngrid);
    333   ALLOCATE (Mjx, double, Ngrid);
    334 
    335   // accumulate the elements of the matrix equation.  We have an equation of the form: Ax = B
    336   // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix,
    337   // and B is an Ngrid vector.  For a cell A(x,y), we need the following elements from each
    338   // star which touches grid cell (x):
    339   //
    340   // Mj  : sum over all measurements of Msys / dMsys^2
    341   // Wj  : sum over all measurements of 1.0 / dMsys^2
    342   // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2
    343   // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2
    344 
    345   // Mjx and Wjx can be calculated by summing over all measurements which touch the cell
    346   // Mj requires looping over stars which touch (x)
    347 
    348   // this is tricky because we need to know both the measurements which touch a cell
    349   // and the stars for which any measurement touches a cell.  we need to accumulate
    350   // sums for each star which touches as cell on both bases.
    351 
    352   int Nsecfilt = GetPhotcodeNsecfilt ();
    353   int thisCode = photcodes[0][0].code;
    354   int Nsec = GetPhotcodeNsec(thisCode);
    355 
    356   for (i = 0; i < Ngrid; i++) {
    357    
    358     for (j = 0; j < Ncatalog; j++) {
    359       memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int));
    360     } 
    361 
    362     // we are looping over the stars, but doing so by looping over the set of measurements:
    363     // every star which touches this grid cell has a measurement in Nlist[i]
    364     for (j = 0; j < Nlist[i]; j++) {
    365      
    366       int mx, c, n, m0, Npts;
    367 
    368       mx = mlist[i][j];
    369       c  = clist[i][j];
    370       n  = catalog[c].measureT[mx].averef;
    371      
    372       // if we have already visited this star, skip the stuff below
    373       if (gotstar[c][n]) continue;
    374       gotstar[c][n] = TRUE;
    375 
    376       // skip stars marked as BAD
    377       if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) {
    378         Nrel ++;
    379         continue;
    380       }
    381 
    382       m0 = catalog[c].average[n].measureOffset;
    383 
    384       // we accumuate an entry for each cell
    385       memset (Wjx, 0, Ngrid*sizeof(double));
    386       memset (Mjx, 0, Ngrid*sizeof(double));
    387       Npts = Mj = Wj = 0.0;
    388 
    389       // if we have not yet visited this star, accumulate the Mj, Wj entries
    390       for (k = 0; k < catalog[c].average[n].Nmeasure; k++) {
    391 
    392         int m, Ng;
    393 
    394         m = m0 + k;
    395 
    396         // skip measurements marked as BAD
    397         if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    398           Nbad ++;
     231  }
     232
     233  // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
     234  for (int code = 0; code < NGridCorr; code++) {
     235    if (!GridCorr[code]) continue;
     236
     237    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
     238      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
     239
     240        // cells without sufficient coverage stay at 0.0
     241        if (GridCorr[code]->nMgrid[ix][iy] < 5) {
     242          GridCorr[code]-> Mgrid[ix][iy] = 0.0;
     243          GridCorr[code]->dMgrid[ix][iy] = 0.0;
    399244          continue;
    400245        }
    401246
    402         // skip images marked as BAD
    403         Mcal = getMcal  (m, c, MAG_CLASS_PSF);
    404         if (isnan(Mcal)) {
    405           Ncal ++;
    406           continue;
    407         }
    408 
    409         // skip mosaics marked as BAD
    410         Mmos = getMmos  (m, c);
    411         if (isnan(Mmos)) {
    412           Nmos ++;
    413           continue;
    414         }
    415 
    416         // skip mosaics marked as BAD
    417         Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    418         if (isnan(Mgrp)) {
    419           Ngrp ++;
    420           continue;
    421         }
    422 
    423         // select the color- and airmass-corrected observed magnitude for this star
    424         // XXX need to be able to turn off the color-correction until initial average mags are found
    425         Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF);
    426         if (isnan(Msys)) {
    427           Nsys++;
    428           continue;
    429         }
    430 
    431         // mag-error for this measurement
    432         Merr =  MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    433 
    434         // disable Wsys for now
    435         Wsys = TRUE ? 1.0 : 1.0 / SQ(Merr);
    436 
    437         Ng = gridmeas[c][m];
    438         if (Ng == -1) continue;  // skip measurements which do not touch any cell
    439 
    440         Mj += Msys * Wsys;  // we are only including measurements touching this cell
    441         Wj += Wsys;  // we are only including measurements touching this cell
    442         Npts ++;
    443         Ngood ++;
    444 
    445         Mjx[Ng] += Msys * Wsys;  // we are only including measurements touching cell (x)
    446         Wjx[Ng] += Wsys;  // we are only including measurements touching cell (x)
    447       }
    448 
    449       // some stars will not have any valid measurements, skip these
    450       if (Npts == 0) continue;
    451 
    452       B[i][0] += Mj*Wjx[i]/Wj - Mjx[i];
    453       A[i][i] -= Wjx[i];
    454       for (k = 0; k < Ngrid; k++) {
    455         A[i][k] += Wjx[i]*Wjx[k]/Wj;
    456         // fprintf (stderr, "%3.0f ", Wjx[k]);
     247        float Mgrid  = GridCorr[code]-> Mgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid
     248        float Mgrid2 = GridCorr[code]-> dMgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid^2
     249
     250        float r = GridCorr[code]->nMgrid[ix][iy] / (float) (GridCorr[code]->nMgrid[ix][iy] - 1.0); // pop -> sample stdev
     251
     252        GridCorr[code]-> Mgrid[ix][iy] = Mgrid;
     253        GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
     254        fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy,
     255                 GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
    457256      }
    458257    }
    459258  }
    460 
    461   if (1) {
    462 
    463     FILE *f;
    464     Header theader;
    465     Matrix matrix;
    466 
    467     /* we are writing to this file */
    468     f = fopen ("matrix.fits", "w");
    469     if (f == (FILE *) NULL) {
    470       fprintf (stderr, "cannot open matrix.fits for output\n");
    471       return;
    472     }
    473 
    474     /* save grid mag values */
    475     gfits_init_header (&theader);
    476     theader.Naxes = 2;
    477     theader.Naxis[0] = Ngrid;
    478     theader.Naxis[1] = Ngrid;
    479     theader.bitpix   = -32;
    480     gfits_create_Theader (&theader, "IMAGE");
    481     gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX");
    482     gfits_create_matrix  (&theader, &matrix);
    483     for (i = 0; i < Ngrid; i++) {
    484       for (j = 0; j < Ngrid; j++) {
    485         gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]);
    486       }
    487     }
    488     gfits_fwrite_header (f, &theader);
    489     gfits_fwrite_matrix (f, &matrix);
    490     gfits_free_matrix (&matrix);
    491 
    492     gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE");
    493     gfits_create_matrix  (&theader, &matrix);
    494     for (i = 0; i < Ngrid; i++) {
    495       for (j = 0; j < Ngrid; j++) {
    496         gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]);
    497       }
    498     }
    499     gfits_fwrite_header (f, &theader);
    500     gfits_fwrite_matrix (f, &matrix);
    501     gfits_free_matrix (&matrix);
    502     fclose (f);
    503   }
    504 
    505   dgaussjordan (A, B, Ngrid, 1);
    506 
    507   fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys);
    508 
    509   for (i = 0; i < Ngrid; i++) {
    510     gridM[i] = B[i][0];
    511     gridS[i] = sqrt(A[i][i]);
    512     gridN[i] = Ngood;
    513   }
    514 
    515   free (Wjx);
    516   free (Mjx);
    517 
    518   for (i = 0; i < Ngrid; i++) {
    519     free (A[i]);
    520   }
    521   free (A);
    522   free (B);
    523 
    524   for (i = 0; i < Ncatalog; i++) {
    525     free (gotstar[i]);
    526     free (gridmeas[i]);
    527   } 
    528   free (gotstar);
    529   free (gridmeas);
    530 }
    531 
    532 /* determine Mgrid values for all grid bins */
    533 void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    534 
    535   int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit;
    536   double *list, *dlist;
    537   float Msys, Mrel, Mcal, Mmos, Mgrp;
    538 
    539   StatType stats;
    540   liststats_setmode (&stats, "INNER_WTMEAN");
     259  return;
     260}
     261
     262float getMgrid (Measure *measure) {
     263
     264  if (!GRID_ZEROPT) return 0.0;
     265  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
     266
     267  int code = measure->photcode;
     268  if (code <= 0) return 0.0;
     269  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
     270       
     271  GridCorrectionType *grid = GridCorr[code];
     272  if (!grid) return 0.0; // does not match one of our images, skip
    541273 
    542   if (!USE_GRID) return;
    543 
    544   int Nsecfilt = GetPhotcodeNsecfilt ();
    545 
    546   Nmax = Nlist[0];
    547   for (i = 0; i < Ngrid; i++) {
    548     Nmax = MAX (Nmax, Nlist[i]);
    549   }
    550   ALLOCATE (list, double, Nmax);
    551   ALLOCATE (dlist, double, Nmax);
    552 
    553   Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0;
    554 
    555   for (i = 0; i < Ngrid; i++) {
    556 
    557     N = 0;
    558     for (j = 0; j < Nlist[i]; j++) {
    559      
    560       m = mlist[i][j];
    561       c = clist[i][j];
    562      
    563       if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    564         Nbad ++;
    565         continue;
    566       }
    567       Mcal = getMcal  (m, c, MAG_CLASS_PSF);
    568       if (isnan(Mcal)) {
    569         Ncal ++;
    570         continue;
    571       }
    572       Mmos = getMmos  (m, c);
    573       if (isnan(Mmos)) {
    574         Nmos ++;
    575         continue;
    576       }
    577       Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    578       if (isnan(Mgrp)) {
    579         Ngrp ++;
    580         continue;
    581       }
    582       Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    583       if (isnan(Mrel)) {
    584         Nrel ++;
    585         continue;
    586       }
    587      
    588       n = catalog[c].measureT[m].averef;
    589       Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
    590       if (isnan(Msys)) {
    591         Nsys++;
    592         continue;
    593       }
    594       list[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
    595       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    596       N++;
    597     }
    598 
    599     // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED
    600     if (gridV[i] == GRID_REFERENCE) {
    601       gridM[i] = 0.0;
    602       gridS[i] = 0.0;
    603       gridN[i] = N;
    604       continue;
    605     }
    606 
    607     // until we have enough valid measurements on this grid cell, skip it
    608     if (N < GRID_TOOFEW) {
    609       gridV[i] = GRID_FROZEN;
    610       continue;
    611     }
    612 
    613     liststats (list, dlist, NULL, N, &stats);
    614     gridM[i] = stats.mean;
    615     gridS[i] = stats.sigma;
    616     gridN[i] = N;
    617     gridV[i] = GRID_FITTED;
    618     Nfit++;
    619   }
    620 
    621   fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys);
    622 
    623   free (list);
    624   free (dlist);
    625 }
    626 
    627 void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    628 
    629   int i, j, m, c, n, N, Narea;
    630   float Msys, Mrel, Mcal, Mmos, Mgrp;
    631   double *xlist, *Mlist, *dlist, *ylist;
    632   Graphdata graphdata;
    633 
    634   if (!USE_GRID) return;
    635 
    636   int Nsecfilt = GetPhotcodeNsecfilt ();
    637 
    638   N = 0;
    639   for (i = 0; i < Ngrid; i++)
    640     N += Nlist[i];
    641 
    642   ALLOCATE (xlist, double, N);
    643   ALLOCATE (ylist, double, N);
    644   ALLOCATE (Mlist, double, N);
    645   ALLOCATE (dlist, double, N);
    646 
    647   Narea = 0;
    648   N = 0;
    649   for (i = 0; i < Ngrid; i++) {
    650     for (j = 0; j < Nlist[i]; j++) {
    651      
    652       m = mlist[i][j];
    653       c = clist[i][j];
    654      
    655       if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    656         Narea ++;
    657         continue;
    658       }
    659       Mcal  = getMcal  (m, c, MAG_CLASS_PSF);
    660       if (isnan(Mcal)) continue;
    661       Mmos  = getMmos  (m, c);
    662       if (isnan(Mmos)) continue;
    663       Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    664       if (isnan(Mgrp)) continue;
    665       Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    666       if (isnan(Mrel)) continue;
    667 
    668       n = catalog[c].measureT[m].averef;
    669       Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
    670 
    671       xlist[N] = Xmeas[c][m];
    672       ylist[N] = Ymeas[c][m];
    673       Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
    674       dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i];
    675       N++;
    676     }
    677   }
    678 
    679   fprintf (stderr, "skipped %d meas for area\n", Narea);
    680 
    681   plot_defaults (&graphdata);
    682   graphdata.ymin = PlotdMmin;
    683   graphdata.ymax = PlotdMmax;
    684   plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "%s.XdM.png", OUTROOT);
    685   plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT);
    686   plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT);
    687 
    688   plot_defaults (&graphdata);
    689   plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT);
    690 
    691   free (ylist);
    692   free (xlist);
    693   free (Mlist);
    694   free (dlist);
    695 
    696 }
    697 
    698 void dump_grid () {
    699 
    700   off_t i, Nimage;
    701   int j, Nbytes, Nformat;
    702   FILE *f;
    703   Header header, theader;
    704   Matrix matrix;
    705   Mosaic *refmosaic;
    706   char *filename;
    707   char formatline[32], key[32], value[64];
    708    
    709   Nbytes = strlen (OUTROOT) + 6;
    710   ALLOCATE (filename, char, Nbytes);
    711   snprintf (filename, Nbytes, "%s.fits", OUTROOT);
    712 
    713   /* select reference mosaic image */
    714   // off_t *imlist = SelectRefMosaic (&refmosaic, &Nimage); return value ignored
    715   SelectRefMosaic (&refmosaic, &Nimage);
    716 
    717   /* we are writing to this file */
    718   f = fopen (filename, "w");
    719   if (f == (FILE *) NULL) {
    720     fprintf (stderr, "cannot open %s for output\n", filename);
    721     free (filename);
    722     return;
    723   }
    724 
    725   /* create empty phu */
    726   gfits_init_header (&header);
    727   header.extend = TRUE;
    728   gfits_create_header (&header);
    729   gfits_create_matrix (&header, &matrix);
    730   gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1,  Nimage + 3);
    731   gfits_modify (&header, "FILTER", "%s", 1, photcodes[0][0].name);  // XXXX note that this expects a single photcode, enforced in initialize.d
    732   gfits_modify_alt (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
    733 
    734   // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system
    735   // Note that config must have been loaded (and not freed) above.
    736   ScanConfig (config, "NFORMAT", "%d", 1, &Nformat);
    737   for (i = 1; i <= Nformat; i++) {
    738       ScanConfig (config, "FORMAT", "%s", i, formatline);
    739       sscanf (formatline, "%s %s", key, value);
    740       gfits_modify (&header, key, "%s", 1, value);
    741   }
    742 
    743   gfits_fwrite_header (f, &header);
    744   gfits_fwrite_matrix (f, &matrix);
    745   gfits_free_matrix (&matrix);
    746 
    747   /* save grid mag values */
    748   gfits_init_header (&theader);
    749   theader.Naxes = 2;
    750   theader.Naxis[0] = gridX;
    751   theader.Naxis[1] = gridY;
    752   theader.bitpix   = -32;
    753   gfits_create_Theader (&theader, "IMAGE");
    754   gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
    755   gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
    756   gfits_create_matrix  (&theader, &matrix);
    757   for (i = 0; i < gridX; i++) {
    758     for (j = 0; j < gridY; j++) {
    759       gfits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]);
    760     }
    761   }
    762   write_coords (&theader, &refmosaic[0].coords);
    763   gfits_fwrite_header (f, &theader);
    764   gfits_fwrite_matrix (f, &matrix);
    765   gfits_free_matrix (&matrix);
    766 
    767   /* save grid Nmeas values */
    768   gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
    769   gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
    770   gfits_create_matrix  (&theader, &matrix);
    771   for (i = 0; i < gridX; i++) {
    772     for (j = 0; j < gridY; j++) {
    773       gfits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]);
    774     }
    775   }
    776   write_coords (&theader, &refmosaic[0].coords);
    777   gfits_fwrite_header (f, &theader);
    778   gfits_fwrite_matrix (f, &matrix);
    779   gfits_free_matrix (&matrix);
    780 
    781   /* save grid sigma values */
    782   gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
    783   gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
    784   gfits_create_matrix  (&theader, &matrix);
    785   for (i = 0; i < gridX; i++) {
    786     for (j = 0; j < gridY; j++) {
    787       gfits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]);
    788     }
    789   }
    790   write_coords (&theader, &refmosaic[0].coords);
    791   gfits_fwrite_header (f, &theader);
    792   gfits_fwrite_matrix (f, &matrix);
    793   gfits_free_matrix (&matrix);
    794 
    795   /* calculate value for each CCD pixel, write out CCD images */
    796   /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
    797   for (i = 0; i < camera.Nchip; i++) {
    798     int N, ix, iy, x, y, X, Y, bin;
    799 
    800     N = ccdnum[i];
    801 
    802     gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
    803     gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
    804     gfits_modify (&theader, "NX", "%d", 1, camera.Nx);
    805     gfits_modify (&theader, "NY", "%d", 1, camera.Ny);
    806      
    807     theader.Naxis[0] = RELPHOT_GRID_X;
    808     theader.Naxis[1] = RELPHOT_GRID_Y;
    809     gfits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X);
    810     gfits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y);
    811     gfits_create_matrix  (&theader, &matrix);
    812 
    813     for (Y = 0; Y < RELPHOT_GRID_Y; Y++) {
    814       for (X = 0; X < RELPHOT_GRID_X; X++) {
    815              
    816         /* normalize X & Y */
    817         x = X;
    818         if (camera.Fx[N]) x = RELPHOT_GRID_X - X - 1;
    819         y = Y;
    820         if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1;
    821              
    822         /* coordinates in the grid */
    823         ix = x + camera.Ox[N]*RELPHOT_GRID_X;
    824         iy = y + camera.Oy[N]*RELPHOT_GRID_Y;
    825              
    826         bin = ix + iy*gridX;
    827         gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]);
    828       }
    829     }
    830     gfits_fwrite_header (f, &theader);
    831     gfits_fwrite_matrix (f, &matrix);
    832     gfits_free_matrix (&matrix);
    833   }
    834 
    835   free (filename);
    836 }
    837 
    838 void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
    839 
    840   int i, j;
    841   double x, y, r, d, X, Y, dx, dy;
    842   double V00, V01, V10, V11;
    843   double wV00, wV01, wV10, wV11;
    844   double dV00, dV01, dV10, dV11;
    845   double v1, v2, value;
    846   int ix, iy, N;
    847 
    848   for (i = 0; i < Nx; i++) {
    849     for (j = 0; j < Ny; j++) {
    850       x = i * RELPHOT_GRID_BINNING / 2;
    851       y = j * RELPHOT_GRID_BINNING / 2;
    852       XY_to_RD (&r, &d, x, y, ccd);
    853       RD_to_XY (&X, &Y, r, d, gcoords);
    854 
    855       X = X / RELPHOT_GRID_BINNING;
    856       Y = Y / RELPHOT_GRID_BINNING;
    857 
    858       ix = (int) X;
    859       dx = X - ix;
    860       iy = (int) Y;
    861       dy = Y - iy;
    862 
    863       if (ix < 0) continue;
    864       if (iy < 0) continue;
    865       if (ix >= gridX) continue;
    866       if (iy >= gridY) continue;
    867 
    868       N = ix + iy*gridX;
    869       V00 = gridM[N];
    870       V10 = gridM[N + 1];
    871       V01 = gridM[N + gridX];
    872       V11 = gridM[N + gridX + 1];
    873 
    874       dV00 = gridS[N];
    875       dV10 = gridS[N + 1];
    876       dV01 = gridS[N + gridX];
    877       dV11 = gridS[N + gridX + 1];
    878 
    879       wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
    880       wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
    881       wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
    882       wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);
    883 
    884       v1 = wV00*V00*(1 + dx*dy - dx - dy) +
    885         wV10*V10*(dx - dx*dy) +
    886         wV01*V01*(dy - dx*dy) +
    887         wV11*V11*(dx*dy);
    888 
    889       v2 = wV00*(1 + dx*dy - dx - dy) +
    890         wV10*(dx - dx*dy) +
    891         wV01*(dy - dx*dy) +
    892         wV11*(dx*dy);
    893 
    894       value = v1 / v2;
    895       buffer[j*Nx + i] = value;
    896     }
    897   }
    898 }
     274  // edge effects could cause some positions to be slightly out of range
     275  // probably should trap extreme outliers
     276  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
     277  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
     278
     279  float Mgrid = grid-> Mgrid[ix][iy];
     280  return Mgrid;
     281}
     282
     283float getMgridTiny (MeasureTiny *measure) {
     284
     285  if (!GRID_ZEROPT) return 0.0;
     286  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
     287
     288  int code = measure->photcode;
     289  if (code <= 0) return 0.0;
     290  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
     291       
     292  GridCorrectionType *grid = GridCorr[code];
     293  if (!grid) return 0.0; // does not match one of our images, skip
     294 
     295  // edge effects could cause some positions to be slightly out of range
     296  // probably should trap extreme outliers
     297  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
     298  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
     299
     300  float Mgrid = grid-> Mgrid[ix][iy];
     301  return Mgrid;
     302}

  • branches/eam_branches/relphot.20210521/src/ImageOps.c

    r41557 r41603  
    319319}
    320320
    321 int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) {
    322 
    323   int ccdnum;
    324   double X, Y;
    325   char *pname, *filter, *p, base[256];
    326 
    327   /* identify the ccd on the basis of the photcode name */
    328   pname = GetPhotcodeNamebyCode (image[idx].photcode);
    329 
    330   // skip measurements which do not match one of the requested photcodes (
    331   // (do we not already exclude in bcatalog -- maybe needed for reload_objects?
    332   int Ns = GetActivePhotcodeIndex (measure[0].photcode);
    333   if (Ns < 0) return FALSE;
    334 
    335   filter = photcodes[Ns][0].name;
    336   sprintf (base, "%s.%s.", MOSAICNAME, filter);
    337   if (strncmp (pname, base, strlen (base))) return (FALSE);
    338   p = pname + strlen(base);
    339 
    340   /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR
    341      for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
    342      and nn is an integer -> ccdnum */
    343      
    344   /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
    345 
    346   while (!isdigit(*p) && *p) p++;
    347   if (*p == 0) return (FALSE);
    348   ccdnum = atoi (p);
    349 
    350   /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
    351      constructed, there will be null values for undefined ccdnums */
    352 
    353   // old code to add this measurement to the grid cell for this chip
    354   // ave = measureT[0].averef;
    355   // ra  = catalog[cat].averageT[ave].R - measureT[0].dR / 3600.0;
    356   // dec = catalog[cat].averageT[ave].D - measureT[0].dD / 3600.0;
    357   // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    358 
    359   // XXX we can now use these values (but need to be careful about old formats)
    360   X = measure[0].Xccd;
    361   Y = measure[0].Yccd;
    362   setGridMeasure (meas, cat, X, Y, ccdnum);
    363 
    364   return (TRUE);
    365 }
    366 
    367321void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) {
    368322
    369323  off_t idx, ID;
    370   int status;
    371324  MeasureTiny *measure;
    372325 
    … …  
    380333  }
    381334  catalog[cat].measureT[meas].myDet = TRUE;
    382 
    383   if (USE_GRID) {
    384     status = findCCD (idx, meas, cat, measure);
    385     if (!status) {
    386       if (VERBOSE2) fprintf (stderr, "failed to determine CCD for "OFF_T_FMT", %d\n",  meas, cat);
    387       return;
    388     }
    389   }
    390335
    391336  // index for (catalog, measure) -> image
    … …  
    690635          continue;
    691636      }
    692       float Mgrid = getMgrid (m, c);
     637      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    693638      if (isnan(Mgrid)) {
    694639          Ngrid++;
    … …  
    889834          continue;
    890835      }
    891       float Mgrid = getMgrid (m, c);
     836      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    892837      if (isnan(Mgrid)) {
    893838          Ngrid++;
    … …  
    11791124      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    11801125      if (isnan(Mgrp)) continue;
    1181       Mgrid = getMgrid (m, c);
     1126      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    11821127      if (isnan(Mgrid)) continue;
    11831128      N++;

  • branches/eam_branches/relphot.20210521/src/MagResidSave.c

    r41561 r41603  
    192192  CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename);
    193193
     194  FitDataSetFree (&psfStars);
     195
    194196  return TRUE;
    195197}
    … …  
    225227    if (isnan(Mgrp)) continue;
    226228
    227     float Mgrid = getMgrid (m, c);
     229    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    228230    if (isnan(Mgrid)) continue;
    229231

  • branches/eam_branches/relphot.20210521/src/MosaicOps.c

    r41559 r41603  
    5050# undef COMPARE
    5151
    52 }
    53 
    54 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
    55 void initMosaics_old (Image *image, off_t Nimage) {
    56 
    57   off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
    58   unsigned int start, stop;
    59   char *pname;
    60 
    61   if (!MOSAIC_ZEROPT) return;
    62 
    63   Nmosaic = 0;
    64   NMOSAIC = 10;
    65   ALLOCATE (mosaic, Mosaic, NMOSAIC);
    66 
    67   ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
    68   ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    69   ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    70 
    71   ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
    72 
    73   /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
    74    * entry in the image table that represents this mosaic.  Instead, it is an
    75    * internal construct that defines a group of related images
    76    */
    77 
    78   /* generate list of unique mosaics */
    79   for (i = 0; i < Nimage; i++) {
    80     ImageToMosaic[i] = -1;
    81 
    82     /* select valid mosaic images by photcode */
    83     pname = GetPhotcodeNamebyCode (image[i].photcode);
    84     status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    85     if (status) continue;
    86 
    87     /* set image time range */
    88     start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    89     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    90 
    91     /* find existing mosaic with this time range */
    92     found = FALSE;
    93     for (j = 0; !found && (j < Nmosaic); j++) {
    94       if (stop  < mosaic[j].start) continue;
    95       if (start > mosaic[j].stop)  continue;
    96       found = TRUE;
    97 
    98       // add reference from image to mosaic
    99       ImageToMosaic[i] = j;
    100 
    101       /* add image to mosaic image list */
    102       MosaicToImage[j][MosaicN_Image[j]] = i;
    103       MosaicN_Image[j] ++;
    104       if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
    105         MosaicN_IMAGE[j] += 10;
    106         REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
    107       }
    108 
    109     }
    110     if (found) continue;
    111    
    112     /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
    113     mosaic[Nmosaic].start     = start;
    114     mosaic[Nmosaic].stop      = stop;
    115     mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    116     mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
    117     mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    118     mosaic[Nmosaic].dMsys     = 0.0;
    119     mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
    120     mosaic[Nmosaic].flags     = image[i].flags;
    121     mosaic[Nmosaic].secz      = image[i].secz;
    122     mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
    123 
    124     // XXX do we need to do something about flag consistency across a mosaic?
    125 
    126     /* add image to mosaic image list */
    127     MosaicN_IMAGE[Nmosaic] = 10;
    128     MosaicN_Image[Nmosaic] = 1;
    129     ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
    130     MosaicToImage[Nmosaic][0] = i;
    131 
    132     // add reference from image to mosaic
    133     ImageToMosaic[i] = Nmosaic;
    134    
    135     Nmosaic ++;
    136     if (Nmosaic == NMOSAIC) {
    137       NMOSAIC += 10;
    138       REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
    139       REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
    140       REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    141       REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    142     }
    143   }
    144 
    145   // free this or not?
    146   free (MosaicN_IMAGE);
    147 
    148   initMosaicGrid (image, Nimage);
    149 
    150   fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
    151   return;
    15252}
    15353
    … …  
    318218  free (startTimesMosaic);
    319219
    320   initMosaicGrid (subset, Nsubset);
     220  initMosaicMcal (subset, Nsubset);
    321221
    322222  fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
    … …  
    460360
    461361  if (mergeMcal) {
    462     initMosaicGrid (image, Nimage);
     362    initMosaicMcal (image, Nimage);
    463363  }
    464364
    … …  
    511411}
    512412
     413// this function sets mosaic->coords to the median of the individual chips.  This
     414// coordinate frame is used by the parallel region analysis to assign exposures (mosaics)
     415// to specific machines by a single center (rather than individual chips)
    513416void setMosaicCenters (Image *image, off_t Nimage) {
    514417  OHANA_UNUSED_PARAM(Nimage);
    … …  
    599502    mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
    600503    mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
    601 
    602     mosaic[i].McalPSF   = 0.0;
    603     mosaic[i].McalAPER  = 0.0;
    604     mosaic[i].dMcal     = 0.0;
    605     mosaic[i].dMsys     = 0.0;
    606     mosaic[i].McalChiSq = 0.0;
    607504  }
    608505  return;
    609506}
    610507
    611 void initMosaicGrid (Image *image, off_t Nimage) {
     508void initMosaicMcal (Image *image, off_t Nimage) {
    612509  OHANA_UNUSED_PARAM(Nimage);
    613510
    614   /* find max dR, dD range for all mosaics */
    615   /* define mosaic.coords to cover dR, dD */
    616   /* send results to initGridBins */
    617 
    618   off_t i, j, m, NX, NY;
    619   int dXmax, dYmax;
    620   double dS, dX, dY;
    621   double R, D, Rmin, Rmax, Dmin, Dmax;
    622511  double McalPSF, McalAPER, dMcal, McalChiSq;
    623512
    624513  fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
    625514
    626   dXmax = dYmax = 0.0;
    627   for (i = 0; i < Nmosaic; i++) {
    628     Dmin = Rmin =  360.0;
    629     Dmax = Rmax = -360.0;
    630     dS = 0.0;
     515  for (off_t i = 0; i < Nmosaic; i++) {
    631516    McalPSF = McalAPER = dMcal = McalChiSq = 0;
    632     for (j = 0; j < MosaicN_Image[i]; j++) {
    633       m = MosaicToImage[i][j];
    634       NX = image[m].NX;
    635       NY = image[m].NY;
    636       dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
    637 
    638       OhanaProjection proj = GetProjection (image[m].coords.ctype);
    639       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    640         XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
    641       } else {
    642         XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
    643       }
    644       Rmin = MIN (Rmin, R);
    645       Rmax = MAX (Rmax, R);
    646       Dmin = MIN (Dmin, D);
    647       Dmax = MAX (Dmax, D);
    648 
    649       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    650         XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
    651       } else {
    652         XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
    653       }
    654       Rmin = MIN (Rmin, R);
    655       Rmax = MAX (Rmax, R);
    656       Dmin = MIN (Dmin, D);
    657       Dmax = MAX (Dmax, D);
    658 
    659       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    660         XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
    661       } else {
    662         XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
    663       }
    664       Rmin = MIN (Rmin, R);
    665       Rmax = MAX (Rmax, R);
    666       Dmin = MIN (Dmin, D);
    667       Dmax = MAX (Dmax, D);
    668 
    669       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    670         XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
    671       } else {
    672         XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
    673       }
    674       Rmin = MIN (Rmin, R);
    675       Rmax = MAX (Rmax, R);
    676       Dmin = MIN (Dmin, D);
    677       Dmax = MAX (Dmax, D);
    678 
    679       // XXX : this probably does not handle mosaics at RA = 0,360 well
    680 
    681       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    682 
    683       // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but
    684       // solve for a single offset for each exposure (Mosaic.Mcal). 
    685       // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
    686       // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
     517    for (off_t j = 0; j < MosaicN_Image[i]; j++) {
     518      off_t m = MosaicToImage[i][j];
    687519
    688520      if (!isfinite(image[m].McalPSF)) {
    689         image[m].McalPSF = 0.0;
    690         image[m].McalAPER = 0.0;
    691         image[m].dMcal = 0.0;
     521        image[m].McalPSF   = 0.0;
     522        image[m].McalAPER  = 0.0;
     523        image[m].dMcal     = 0.0;
    692524        image[m].McalChiSq = 0.0;
    693525        fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
    … …  
    703535      image[m].dMcal     = NAN;
    704536      image[m].McalChiSq = NAN;
    705 
    706     }
    707     dS /= MosaicN_Image[i];
    708 
    709     InitCoords (&mosaic[i].coords, "DEC--TAN");
    710     mosaic[i].coords.crval1 = Rmin;
    711     mosaic[i].coords.crval2 = Dmin;
    712     mosaic[i].coords.cdelt1 = dS;
    713     mosaic[i].coords.cdelt2 = dS;
    714 
    715     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
     537    }
    716538
    717539    mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
    … …  
    720542    mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
    721543  }
    722   if (!USE_GRID) return;
    723 
    724   dXmax = MAX (dXmax, dX);
    725   dYmax = MAX (dYmax, dY);
    726   initGrid (dXmax, dYmax);
    727544  return;
    728545}
    … …  
    12311048      continue;
    12321049    }
    1233     float Mgrid = getMgrid (m, c);
     1050    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    12341051    if (isnan(Mgrid)) {
    12351052      info->Ngrid ++;
    … …  
    13821199    float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
    13831200    // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    1384     float Mgrid    = getMgrid (m, c);
     1201    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);
    13851202    float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    13861203    float Mflat    = getMflat (m, c, flatcorr, catalog);
    … …  
    14811298           summary.Nskip);
    14821299  free (threadinfo);
     1300  // XXX SetMmosInfoFree (&summary);
    14831301
    14841302  npass_output ++;
    … …  
    18111629      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    18121630      if (isnan(Mgrp)) continue;
    1813       Mgrid = getMgrid (m, c);
     1631      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    18141632      if (isnan(Mgrid)) continue;
    18151633      Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);

  • branches/eam_branches/relphot.20210521/src/StarOps.c

    r41556 r41603  
    495495        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    496496        if (isnan(Mgrp)) continue;
    497         float Mgrid = getMgrid (m, i);
     497        float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    498498        if (isnan(Mgrid)) continue;
    499499
    … …  
    737737          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    738738          if (isnan(Mgrp)) { Ngrp ++; continue; }
    739           Mgrid = getMgrid (m, i);
     739          Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    740740          if (isnan(Mgrid)) { Ngrid ++; continue; }
    741741
    … …  
    788788          Mmos  = getMmos  (m, i);
    789789          if (isnan(Mmos)) continue;
    790           Mgrid = getMgrid (m, i);
     790          Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    791791          if (isnan(Mgrid)) continue;
    792792
    … …  
    869869        Mmos = getMmos  (m, i);
    870870        if (isnan(Mmos)) { continue; }
    871         Mgrid = getMgrid (m, i);
     871        Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    872872        if (isnan(Mgrid)) { continue;}
    873873        N++;
    … …  
    11051105
    11061106      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
    1107       // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
     1107      // float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
    11081108      // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
    11091109

  • branches/eam_branches/relphot.20210521/src/TGroupOps.c

    r41562 r41603  
    853853      continue;
    854854    }
    855     float Mgrid = getMgrid (m, c);
     855    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    856856    if (isnan(Mgrid)) {
    857857      info->Ngrid ++;
    … …  
    10131013    float Mmos     = getMmos  (m, c);                    // mosaic zero point
    10141014
    1015     float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
     1015    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
    10161016    float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
    10171017    float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
    … …  
    11241124  free (threadinfo);
    11251125
     1126  // XXX not allocated SetMgrpInfoFree (&summary);
     1127
    11261128  npass_output ++;
    11271129
    … …  
    12681270
    12691271static float MinMaxChiSq = NAN;
     1272static float MaxMaxChiSq = 6.0;
    12701273static float MinMaxScatter = NAN;
     1274static float MaxMaxScatter = 0.1;
    12711275
    12721276void clean_tgroups () {
    … …  
    13011305  float ChiSqUpper90 = stats.Upper90;
    13021306  if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
    1303   float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90);
     1307  float MaxChiSq = MIN(MaxMaxChiSq, MAX (MinMaxChiSq, ChiSqUpper90)); // the upper limit should be between MinMaxChiSq and MaxMaxChiSq
    13041308
    13051309  liststats (slist, NULL, NULL, N, &stats);
    13061310  float ScatterUpper90 = stats.Upper90;
    13071311  if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
    1308   float MaxScatter = MAX (MinMaxScatter, ScatterUpper90);
     1312  float MaxScatter = MIN(MaxMaxScatter, MAX (MinMaxScatter, ScatterUpper90));
    13091313
    13101314  fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
    … …  
    15011505
    15021506        float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
    1503         // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
     1507        // float Mgrid    = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?)
    15041508        // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
    15051509

  • branches/eam_branches/relphot.20210521/src/ZeroPointModes.c

    r41598 r41603  
    4949      clean_mosaics(); // do this on every pass
    5050    }
    51     if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) {
     51    if ((CurrentLoop >= 18) && (CurrentLoop < 21)) {
    5252      // after iterating a few times on the TGroup & Mosaic zero points:
    5353      // * identify the good / bad Mosaics
    … …  
    6060      clean_images(); // do this on every pass
    6161    }
     62    if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) {
     63      // after iterating a few times on the TGroup & Mosaic zero points:
     64      // * identify the good / bad Mosaics
     65      // * fit the chips from the bad Mosaics
     66      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
     67      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights
     68      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights
     69      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
     70      clean_tgroups(); // do this on every pass
     71      clean_mosaics(); // do this on every pass
     72      clean_images(); // do this on every pass
     73    }
    6274  }
    6375
    … …  
    6880      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
    6981    }
    70     if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
     82    if ((CurrentLoop > 4) && (CurrentLoop <= 8)) {
    7183      // after iterating a few times on the TGroup zero points:
    7284      // * identify the photometric nights
    … …  
    7587      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
    7688      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT;
     89      clean_tgroups(); // do this on every pass or just sometimes?
     90    }
     91    if ((CurrentLoop > 8) && (CurrentLoop <= 999)) {
     92      // after iterating a few times on the TGroup zero points:
     93      // * identify the photometric nights
     94      // * fit the images from the non-photometric nights
     95      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
     96      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
     97      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT;
     98      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
    7799      clean_tgroups(); // do this on every pass or just sometimes?
    78100    }

  • branches/eam_branches/relphot.20210521/src/args.c

    r41562 r41603  
    198198  }
    199199
    200   NGRID = 8;
    201   if ((N = get_argument (argc, argv, "-ngrid"))) {
    202     remove_argument (N, &argc, argv);
    203     NGRID = atof (argv[N]);
    204     remove_argument (N, &argc, argv);
    205   }
    206 
    207200  TEST_IMAGE1 = -1;
    208201  if ((N = get_argument (argc, argv, "-test-image1"))) {
    … …  
    387380  }
    388381
    389   // USE_GRID is not valid for all cases, probably should be its own mode...
    390   USE_GRID = FALSE;
     382  // GRID_ZEROPT is not valid for all cases, probably should be its own mode...
     383  GRID_ZEROPT = FALSE;
    391384  if ((N = get_argument (argc, argv, "-grid"))) {
    392385    remove_argument (N, &argc, argv);
    393     USE_GRID = TRUE;
     386    GRID_ZEROPT = TRUE;
     387    GRID_ZPT_MODE  = GRID_ZPT_MODE_NONE; // start with grid off
    394388  }
    395389

  • branches/eam_branches/relphot.20210521/src/assign_images.c

    r41485 r41603  
    3636  if (MOSAIC_ZEROPT) {
    3737    makeMosaics (image, Nimage, FALSE);
    38 
    39     // center coords and Mcal, dMcal, Mchisq for the mosaics
     38    MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
     39
     40    // center coords and zero Mcal, dMcal, Mchisq for the mosaics
    4041    setMosaicCenters (image, Nimage);
    4142    MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
    … …  
    124125      Mosaic *mosaic = getMosaicForImage (j);
    125126      if (mosaic) {
    126         Rc = mosaic->coords.crval1;
    127         Dc = mosaic->coords.crval2;
    128         // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
    129         // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
    130         Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
     127        Rc = mosaic->coords.crval1;
     128        Dc = mosaic->coords.crval2;
     129        // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
     130        // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
     131        Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
    131132      }
    132133    }

  • branches/eam_branches/relphot.20210521/src/bcatalog.c

    r41561 r41603  
    112112
    113113      // require 0x01 in photFlags (fitted with a PSF)
    114       if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
     114// XXX EAM : test is inconsistent -- fix      if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
    115115
    116116      // very loose cut on PSF - Kron

  • branches/eam_branches/relphot.20210521/src/extra.c

    r41467 r41603  
    11# include "relphot.h"
     2# define DISABLE_CORNER_GPC1 0
    23
    34int isMosaicChip (int photcode) {
    … …  
    1415int whichGPC1filter (int photcode) {
    1516
    16 # if (1)
     17# if (DISABLE_CORNER_GPC1)
    1718  // disable the corner chips:
    1819  if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return PS1_g; // g-band
    … …  
    3738int isGPC1chip (int photcode) {
    3839
    39 # if (1)
     40# if (DISABLE_CORNER_GPC1)
    4041  if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return TRUE; // g-band
    4142  if (((photcode >  10100) && (photcode <  10177)) || (photcode == 4200)) return TRUE; // r-band

  • branches/eam_branches/relphot.20210521/src/initialize.c

    r41453 r41603  
    2727  }   
    2828   
    29   if (USE_GRID && (Nphotcodes > 1)) {
    30     fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
    31     exit (1);
    32   }
    33 
    3429  IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    3530  STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    … …  
    6762
    6863    fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
    69     fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n",
    70              RELPHOT_GRID_X,
    71              RELPHOT_GRID_Y,
    72              RELPHOT_GRID_BINNING,
    73              (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
    7464
    7565    fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);

  • branches/eam_branches/relphot.20210521/src/plot_scatter.c

    r40291 r41603  
    4444                Mmos = getMmos  (m, i);
    4545                if (isnan(Mmos)) continue;
    46                 Mgrid = getMgrid (m, i);
     46                Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    4747                if (isnan(Mgrid)) continue;
    4848

  • branches/eam_branches/relphot.20210521/src/reload_catalogs.c

    r41390 r41603  
    126126    freeMosaicBins (1, FALSE);
    127127    freeTGroupBins (1);
    128     freeGridBins (1);
     128    freeGridBins ();
    129129    TIMESTAMP(time7);
    130130

  • branches/eam_branches/relphot.20210521/src/relphot_images.c

    r41556 r41603  
    5454  FlatCorrectionTable *flatcorr = NULL;
    5555
     56  initGridBins (); // allocates the empty array of corrections (elements are build below)
     57
    5658  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
    5759    /* load catalog data from region files (hostID is 0 since we are not a client */
    … …  
    6870    initTGroupBins (catalog, Ncatalog);
    6971
    70     initGridBins   (catalog, Ncatalog);
    7172    initMrel (catalog, Ncatalog);
    7273
    … …  
    7475    MARKTIME("-- set up image indexes: %f sec\n", dtime);
    7576
    76     findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
     77    findMosaics (catalog, Ncatalog, TRUE);
    7778    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    7879
    7980    findTGroups (catalog, Ncatalog);
    8081    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
     82
     83    initGrid (); // allocate grid correction entries for existing photcodes
    8184
    8285    SAVEPLOT = FALSE;
    … …  
    9194    }
    9295
    93     // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
    94     if (USE_GRID) {
    95       int star_toofew;
    96 
    97       // until we finish the grid analysis, do not reject stars out-of-hand based on ID_OBJ_FEW
    98       // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
    99       star_toofew = STAR_TOOFEW;
    100       STAR_TOOFEW = 0;
    101       for (i = 0; i < NGRID; i++) {
    102         STAR_BAD = ID_OBJ_POOR;
    103         setMrel  (catalog, Ncatalog, flatcorr);
    104         STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    105         setMgrid (catalog, flatcorr);
    106       }
    107       STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    108       STAR_TOOFEW = star_toofew;
    109     }
    110 
    11196    /* determine fit values */
    11297    for (i = 0; i < NLOOP; i++) {
    … …  
    119104      setMgrp  (catalog, flatcorr);
    120105
    121       setMgrid (catalog, flatcorr);
     106      setMgrid (catalog, Ncatalog);
    122107      MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    123108   
    124109      if (PLOTSTUFF) {
    125110        plot_scatter (catalog, Ncatalog, flatcorr);
    126         plot_grid (catalog, flatcorr);
    127111        plot_mosaics ();
    128112        plot_images ();
    … …  
    140124    if (PLOTSTUFF) {
    141125      plot_scatter (catalog, Ncatalog, flatcorr);
    142       plot_grid (catalog, flatcorr);
    143126      plot_mosaics ();
    144127      plot_images ();
    … …  
    147130    }
    148131 
    149     if (USE_GRID) dump_grid ();
     132    // if (GRID_ZEROPT) dump_grid ();
    150133
    151134    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    … …  
    169152    freeMosaicBins (Ncatalog, TRUE);
    170153    freeTGroupBins (Ncatalog);
    171     freeGridBins (Ncatalog);
    172 
    173154    // end of if (NLOOP > 0) block : this loop determines the offsets per chip
    174155  } else {
    … …  
    205186  freeImages (db.ftable.buffer);
    206187  freeMosaics ();
     188  freeGridBins ();
     189
    207190  gfits_db_free (&db);
    208191

  • branches/eam_branches/relphot.20210521/src/relphot_objects.c

    r41507 r41603  
    147147      freeImageBins (1, FALSE);
    148148      freeMosaicBins (1, FALSE);
    149       freeGridBins (1);
     149      freeGridBins ();
    150150      continue;
    151151    }
    … …  
    183183    freeImageBins (1, FALSE);
    184184    freeMosaicBins (1, FALSE);
    185     freeGridBins (1);
     185    freeGridBins ();
    186186  }
    187187

  • branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c

    r41556 r41603  
    161161  freeImageBins(Ncatalog, TRUE);
    162162  freeMosaicBins (Ncatalog, TRUE);
    163   freeGridBins (Ncatalog);
     163  freeGridBins ();
    164164  freeImages((char *)image);
    165165  free (image);

  • branches/eam_branches/relphot.20210521/src/setMrelCatalog.c

    r41561 r41603  
    243243      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
    244244      if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
    245       Mgrid = getMgrid (meas, cat);
     245      Mgrid = getMgridTiny (&measureT[k]);
    246246      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
    247247      Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
    … …  
    732732      McalAPER  = USE_MCAL_PSF_FOR_STACK_APER ? getMcal  (measSeq, cat, MAG_CLASS_PSF) : getMcal  (measSeq, cat, MAG_CLASS_KRON);
    733733      Mmos      = getMmos  (measSeq, cat);
    734       Mgrid     = getMgrid (measSeq, cat);
     734      Mgrid     = getMgrid (&measure[meas]);
    735735      // XXX can Mmos and Mgrid exist for stacks?
    736736    }
Note: See TracChangeset for help on using the changeset viewer.