IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.