Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/Makefile =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/Makefile (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/Makefile (revision 37913) @@ -112,4 +112,5 @@ $(SRC)/setastrom.$(ARCH).o \ $(SRC)/initialize_setastrom.$(ARCH).o \ +$(SRC)/cam_correction.$(ARCH).o \ $(SRC)/dcr_correction.$(ARCH).o \ $(SRC)/kh_correction.$(ARCH).o \ @@ -127,4 +128,5 @@ $(SRC)/update_dvo_setastrom.$(ARCH).o \ $(SRC)/update_catalog_setastrom.$(ARCH).o \ +$(SRC)/cam_correction.$(ARCH).o \ $(SRC)/dcr_correction.$(ARCH).o \ $(SRC)/kh_correction.$(ARCH).o \ Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/include/setastrom.h =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/include/setastrom.h (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/include/setastrom.h (revision 37913) @@ -13,11 +13,21 @@ // we have one correction (an image) for each filter and chip typedef struct { - int Nvalues; - int Nfilters; // number of filters - int NchipX; // number of chips in x - int NchipY; // number of chips in y + int Nx; // number of chips in x + int Ny; // number of chips in y + int Nfilter; // number of filters + int Ndir; // number of correction dimensions + + int Nchips; // chip offset (Nx*Ny) + int Ngroup; // direction offset (Nx*Ny*Nfilters) + int Nvalues; // Nx*Ny*Nfilters*Ndir + + int dX; // superpixel size + int dY; // superpixel size + int NxCCD; // number of pixels + int NyCCD; // number of pixels + Matrix **matrix; // allocate an array of pointers - // index = ix + iy*Nx + filter*Nx*Ny + dir*Nx*Ny*Nfilter + // index = ix + iy*Nx + filter*Nchips + dir*Ngroup } CamCorrection; @@ -38,7 +48,9 @@ char *KH_FILE; char *DCR_FILE; +char *CAM_FILE; int KH_RESET; int DCR_RESET; +int CAM_RESET; SkyRegion UserPatch; @@ -76,2 +88,5 @@ int load_dcr_correction (char *filename); int get_dcr_correction (int filter, double *dX, double *dY, float Color); + +int load_cam_correction (char *filename); +int get_cam_correction (int chipID, int filter, float Xccd, float Yccd, double *dX, double *dY); Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/cam_correction.c =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/cam_correction.c (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/cam_correction.c (revision 37913) @@ -6,15 +6,10 @@ // note that we are allocating pointers for XY00 - XY77, but only the octal elements are set (and not the 00,07,70,77 ones) -CamCorrection *load_cam_correction (char *filename) { +static CamCorrection *cam = NULL; - int Ncol; - off_t Nrow; +int load_cam_correction (char *filename) { - char type[16]; - char extname[80]; - + int i, ix, iy, dir, filter; Header header; - Header theader; - FTable ftable; /* open file for input */ @@ -29,17 +24,19 @@ // read the PHU header - if (!gfits_load_header (f, &header)) return (FALSE); + if (!gfits_load_header (f, &header)) return FALSE; int Nbytes = gfits_data_size (&header); fseeko (f, Nbytes, SEEK_CUR); - ftable.header = &theader; + int Ndir, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD; + if (!gfits_scan (&header, "NDIR", "%d", 1, &Ndir)) return FALSE; + if (!gfits_scan (&header, "NFILTER", "%d", 1, &Nfilter)) return FALSE; + if (!gfits_scan (&header, "NX", "%d", 1, &Nx)) return FALSE; + if (!gfits_scan (&header, "NY", "%d", 1, &Ny)) return FALSE; + if (!gfits_scan (&header, "DX", "%d", 1, &dX)) return FALSE; + if (!gfits_scan (&header, "DY", "%d", 1, &dY)) return FALSE; - if (!gfits_scan (&theader, "NDIR", "%d", 1, &Ndir)) return (NULL); - if (!gfits_scan (&theader, "NFILTER", "%d", 1, &Nfilter)) return (NULL); - if (!gfits_scan (&theader, "NX_CHIP", "%d", 1, &Nx)) return (NULL); - if (!gfits_scan (&theader, "NY_CHIP", "%d", 1, &Ny)) return (NULL); - if (!gfits_scan (&theader, "DX_CHIP", "%d", 1, &dX)) return (NULL); - if (!gfits_scan (&theader, "DY_CHIP", "%d", 1, &dY)) return (NULL); + if (!gfits_scan (&header, "NX_CHIP", "%d", 1, &NxCCD)) return FALSE; + if (!gfits_scan (&header, "NY_CHIP", "%d", 1, &NyCCD)) return FALSE; ALLOCATE (cam->matrix, Matrix *, Ndir*Nfilter*Nx*Ny); @@ -47,55 +44,100 @@ cam->Nx = Nx; // for gpc1, should be 8 cam->Ny = Ny; // for gpc1, should be 8 + cam->Nfilter = Nfilter; // for gpc1, should be 5 + cam->Ndir = Ndir; // for gpc1, should be 2 + cam->dX = dX; // superpixel sampling, x-dir cam->dY = dY; // superpixel sampling, y-dir - cam->Ndir = Ndir; // for gpc1, should be 2 - cam->Nfilter = Nfilter; // for gpc1, should be 5 + cam->NxCCD = NxCCD; // pixels per chip + cam->NyCCD = NyCCD; // pixels per chip + + cam->Nchips = Nx * Ny; + cam->Ngroup = Nx * Ny * Nfilter; + cam->Nvalues = Nx * Ny * Nfilter * Ndir; + + ALLOCATE (cam->matrix, Matrix *, cam->Nvalues); + for (i = 0; i < cam->Nvalues; i++) { + cam->matrix[i] = NULL; + } + + int found = FALSE; while (TRUE) { + Header theader; + // load data for this header : if not found, assume we hit the end of the file if (!gfits_load_header (f, &theader)) break; - if (!gfits_scan (&theader, "FILTER", "%d", 1, &filter)) return (FALSE); - if (!gfits_scan (&theader, "DIR", "%d", 1, &dir)) return (FALSE); - if (!gfits_scan (&theader, "X_CHIP", "%d", 1, &ix)) return (FALSE); - if (!gfits_scan (&theader, "Y_CHIP", "%d", 1, &iy)) return (FALSE); + if (!gfits_scan (&theader, "FILTER", "%d", 1, &filter)) return FALSE; + if (!gfits_scan (&theader, "DIR", "%d", 1, &dir)) return FALSE; + if (!gfits_scan (&theader, "X_CHIP", "%d", 1, &ix)) return FALSE; + if (!gfits_scan (&theader, "Y_CHIP", "%d", 1, &iy)) return FALSE; Matrix *matrix = NULL; ALLOCATE (matrix, Matrix, 1); - if (!gfits_load_header (f, matrix, &theader)) break; + if (!gfits_load_matrix (f, matrix, &theader)) break; - int index = ix + iy*Nx + filter*Nx*Ny + dir*Nx*Ny*Nfilter; + int index = ix + iy*cam->Nx + filter*cam->Nchips + dir*cam->Ngroup; myAssert (index >= 0, "index too small"); - myAssert (index < Ndir*Nfilter*Nx*Ny, "index too big"); + myAssert (index < cam->Nvalues, "index too big"); + + myAssert (!cam->matrix[i], "entry already assigned?"); // assert that cam->matrix[index] is NULL? cam->matrix[index] = matrix; + + if (!found) { + + } } - return (TRUE); + return TRUE; } -int get_cam_correction (CamCorrection *cam, int chip, int filter, double *dX, double *dY, float Xccd, float Yccd) { +// input is: +// chipID == octal chip ID (eg 40 for XY40) +// filter == (01234 = grizy) +int get_cam_correction (int chipID, int filter, float Xccd, float Yccd, double *dX, double *dY) { *dX = 0.0; *dY = 0.0; - // split out the chip number (eg 43 for XY43) into X and Y elements - int ix = (int) (chip / 10); - int iy = (int) (chip % 10); + int index, jx, jy, NxModel; + Matrix *matrix; + float *buffer, value; - int dir = 0; // dX - int index = ix + iy*cam->Nx + filter*cam->Nx*cam->Ny + dir*cam->Nx*cam->Ny*cam->Nfilter; + // split out the chipID number (eg 43 for XY43) into X and Y elements + int ix = (int) (chipID / 10); + int iy = (int) (chipID % 10); - Matrix *matrix = cam->matrix[index]; + // dX (dir == 0) + index = ix + iy*cam->Nx + filter*cam->Nchips; - float *buffer = (float *) matrix->buffer; + matrix = cam->matrix[index]; + NxModel = cam->matrix[index]->Naxis[0]; - int jx = MAX(MIN(Nx, Xccd / cam->dX), 0); - int jy = MAX(MIN(Nx, Yccd / cam->dY), 0); + buffer = (float *) matrix->buffer; - float value = buffer[jx + jy*Nx]; + jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0); + jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0); + + value = buffer[jx + jy*NxModel]; *dX = isnan(value) ? 0.0 : value; - return (TRUE); + + // dY (dir == 1) uses the next block + index += cam->Ngroup; + + matrix = cam->matrix[index]; + NxModel = cam->matrix[index]->Naxis[0]; + + buffer = (float *) matrix->buffer; + + jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0); + jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0); + + value = buffer[jx + jy*NxModel]; + + *dY = isnan(value) ? 0.0 : value; + + return TRUE; } - Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/initialize_setastrom.c =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/initialize_setastrom.c (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/initialize_setastrom.c (revision 37913) @@ -4,6 +4,7 @@ fprintf (stderr, "USAGE: setastrom [options]\n"); fprintf (stderr, " options:\n"); + fprintf (stderr, " -KH (file) : supply Koppenhoefer correction splines\n"); fprintf (stderr, " -DCR (file) : supply DCR correction splines\n"); - fprintf (stderr, " -KH (file) : supply Koppenhoefer correction splines\n"); + fprintf (stderr, " -CAM (file) : supply camera-static correction file\n"); fprintf (stderr, " -v : verbose mode\n"); fprintf (stderr, " -region Rmin Rmax Dmin Dmax\n"); @@ -77,4 +78,17 @@ } + CAM_FILE = NULL; + if ((N = get_argument (argc, argv, "-CAM"))) { + remove_argument (N, &argc, argv); + char *tmpfile = strcreate (argv[N]); + CAM_FILE = abspath (tmpfile, DVO_MAX_PATH); + remove_argument (N, &argc, argv); + } + CAM_RESET = FALSE; + if ((N = get_argument (argc, argv, "-CAM-reset"))) { + remove_argument (N, &argc, argv); + CAM_RESET = TRUE; + } + SINGLE_CPT = NULL; if ((N = get_argument (argc, argv, "-cpt"))) { @@ -146,7 +160,6 @@ } - if (!KH_FILE && !DCR_FILE) { - fprintf (stderr, "at least one of -KH and -DCR must be supplied\n"); - exit (1); + if (!KH_FILE && !DCR_FILE && !CAM_FILE) { + fprintf (stderr, "WARNING: none of -CAM, -KH, -DCR supplied\n"); } @@ -159,6 +172,7 @@ fprintf (stderr, "USAGE: setastrom_client -hostID (hostID) -catdir (catdir) -hostdir (hostdir) [options]\n"); fprintf (stderr, " options:\n"); + fprintf (stderr, " -KH (file) : supply Koppenhoefer correction splines\n"); fprintf (stderr, " -DCR (file) : supply DCR correction splines\n"); - fprintf (stderr, " -KH (file) : supply Koppenhoefer correction splines\n"); + fprintf (stderr, " -CAM (file) : supply camera-static correction file\n"); fprintf (stderr, " -region Rmin Rmax Dmin Dmax\n"); fprintf (stderr, " -update-catformat (format) : change database schema on output\n"); @@ -221,4 +235,16 @@ remove_argument (N, &argc, argv); DCR_RESET = TRUE; + } + + CAM_FILE = NULL; + if ((N = get_argument (argc, argv, "-CAM"))) { + remove_argument (N, &argc, argv); + CAM_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + CAM_RESET = FALSE; + if ((N = get_argument (argc, argv, "-CAM-reset"))) { + remove_argument (N, &argc, argv); + CAM_RESET = TRUE; } @@ -282,7 +308,6 @@ if (!CATDIR) usage_setastrom_client(); - if (!KH_FILE && !DCR_FILE) { - fprintf (stderr, "at least one of -KH and -DCR must be supplied\n"); - exit (1); + if (!KH_FILE && !DCR_FILE && !CAM_FILE) { + fprintf (stderr, "WARNING: none of -CAM, -KH, -DCR supplied\n"); } Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_catalog_setastrom.c =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_catalog_setastrom.c (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_catalog_setastrom.c (revision 37913) @@ -121,13 +121,14 @@ double dX_CAM = 0.0; double dY_CAM = 0.0; - if (CAM_FILE && doCAM) { + if (CAM_FILE) { // camera systematic correction ("astroflat") depends on only the X,Y coordinate, // the filter, and the chip - int filtCode = (int)(measure->photcode / 100); // eg, 101 = r + // 10134 : r, XY34 -> filtCode = 1 + int filtCode = (int)((measure->photcode % 100) / 100); // correction is in arcseconds - get_cam_correction (0, chipID, filtCode, measure->Xccd, measure->Yccd, &dX_CAM, &dY_CAM); + get_cam_correction (chipID, filtCode, measure->Xccd, measure->Yccd, &dX_CAM, &dY_CAM); } @@ -260,6 +261,14 @@ } if (DCR_RESET) { - measure[0].XoffDCR = -dX_DCR / pltScale; - measure[0].YoffDCR = -dY_DCR / pltScale; + measure[0].XoffDCR = 0.0; + measure[0].YoffDCR = 0.0; + } + if (CAM_FILE) { + measure[0].XoffCAM = -dX_DCR / pltScale; + measure[0].YoffCAM = -dY_DCR / pltScale; + } + if (CAM_RESET) { + measure[0].XoffDCR = 0.0; + measure[0].YoffDCR = 0.0; } @@ -274,4 +283,8 @@ measure[0].Yfix += measure[0].YoffDCR; } + if (isfinite(measure[0].XoffCAM) && isfinite(measure[0].YoffCAM)) { + measure[0].Xfix += measure[0].XoffCAM; + measure[0].Yfix += measure[0].YoffCAM; + } found ++; } Index: /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_dvo_setastrom.c =================================================================== --- /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_dvo_setastrom.c (revision 37912) +++ /branches/eam_branches/ipp-20150112/Ohana/src/uniphot/src/update_dvo_setastrom.c (revision 37913) @@ -17,6 +17,7 @@ } - if (KH_FILE) load_kh_correction (KH_FILE); - if (DCR_FILE) load_dcr_correction (DCR_FILE); + if (KH_FILE) if (!load_kh_correction (KH_FILE)) { fprintf (stderr, "failed to load KH correction %s\n", KH_FILE ); exit (1); } + if (DCR_FILE) if (!load_dcr_correction (DCR_FILE)) { fprintf (stderr, "failed to load DCR correction %s\n", DCR_FILE); exit (1); } + if (CAM_FILE) if (!load_cam_correction (CAM_FILE)) { fprintf (stderr, "failed to load CAM correction %s\n", CAM_FILE); exit (1); } // determine the populated SkyRegions overlapping the requested area (default depth) @@ -104,8 +105,10 @@ if (KH_FILE) { snprintf (tmpline, 1024, "%s -KH %s", command, KH_FILE); strcpy (command, tmpline); } if (DCR_FILE) { snprintf (tmpline, 1024, "%s -DCR %s", command, DCR_FILE); strcpy (command, tmpline); } + if (CAM_FILE) { snprintf (tmpline, 1024, "%s -CAM %s", command, CAM_FILE); strcpy (command, tmpline); } if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT); strcpy (command, tmpline); } if (KH_RESET) { snprintf (tmpline, 1024, "%s -KH-reset", command); strcpy (command, tmpline); } if (DCR_RESET) { snprintf (tmpline, 1024, "%s -DCR-reset", command); strcpy (command, tmpline); } + if (CAM_RESET) { snprintf (tmpline, 1024, "%s -CAM-reset", command); strcpy (command, tmpline); } fprintf (stderr, "command: %s\n", command);