iterate to find Mcal, image.flags:
+
+
+
+
+
+
+
+ relphot overview
+
+ relphot has two primary purposes:
+
+- calculate Mcal for images / measures
+
- calculate Mrel for stars
+
+
+ relphot can also be used to determine the mosaic grid used to generate
+photometrically corrected flats (-grid option).
+
+
data exclusion
+
+ relphot uses only a subset of the photometry data to calculate
+Mcal and Mrel. In the first stage, calculation of the Mcal values,
+relphot loads the photometry data from each relevant catalog and
+creates an internal subset catalog with the function
+bcatalog, excluding some of the irrelevant data. In
+addition, it uses flags to mark some of the data as invalid for the
+processing.
+
+ bcatalog exclusions
+
+- measure.photcode not equivalent to requested photcode
+
- measure.dophot != 1
+
- measure.Mcat > MAG_LIM
+
- measure.dM > SIGMA_LIM
+
- measure.Minst out of range (ImagMin - ImagMax) [optional]
+
- measure.t out of range (TSTART, TSTOP)
+
+ flagged data
+
+ flagged image data
+
+images can be flagged by setting bits of image.code
+
+stars can be flagged by setting bits of average.code
+
+measures can be flagged by setting bits of measure.flag
+
+image.code
+
+ID_IMAGE_NOCAL : ignore, irrelevant
+ID_IMAGE_POOR : image measured bad
+ID_IMAGE_SKIP : externally known bad
+
+dMcal > VALUE FLAG_IMAGE_SCATTER clean_images
+fabs(Mcal) > VALUE FLAG_IMAGE_ZEROPT clean_images
+
+dMcal > VALUE FLAG_IMAGE_SCATTER clean_mosaics
+fabs(Mcal) > VALUE FLAG_IMAGE_ZEROPT clean_mosaics
+
+ mark_images does not seem to do anything useful?
+
+average.code
+
+Ngood < MEAS_TOOFEW setMrel
+Ngood < MEAS_TOOFEW clean_measures
+ChiSq > STAR_CHISQ clean_stars
+dM > STAR_SCATTER clean_stars
+
+average.code (STAR_BAD) is not saved by relphot: it is set by
+clean_stars, clean_measures, and setMrel, but not setMrelOutput.
+STAR_BAD should only be internal since it depends on the photcode, but
+is not associated with a specific photcode in the data. Just in case,
+it is reset to 0 in setMrelFinal.
+
+measure.flag
+
+X,Y out of range setExclusions
+3 sigma clipping clean_measures
+
+
setting Mrel final value
+
+setMrelFinal is used to set the final average.Mrel values. We do this
+in 4 stages. In each stage, we set the Mrel values for stars which
+have not already been set, based on the current exclusion settings.
+At successive stages, we relax the exclusions, allowing the more
+spurious objects to have a valid Mrel value to be set. In this loop,
+we actually run setMrelOutput twice: once to get the approximate Mrel
+value, then we flag the outlier measurements with clean_measure,
+then we redetermine the Mrel values on this basis, and mark the stars
+for exclusion from the next iteration.
+
+exclude on
+ photcode 0 1 2 3
+ time range 0 1 2 3
+ MEAS_POOR 0 1 2 3
+ MEAS_TOOFEW 0 1 2 3
+ dophot == 10 0 1 2
+ inst mag 0 1 2
+ dophot != 1,2 0 1
+ ID_IMAGE_POOR 0 1
+ ID_IMAGE_SKIP 0 1
+ dophot != 1 0
+ measure.dM 0
+
+for all relphot runs, Mrel is re-calculated, and measures are marked
+at least if they are outliers in mag or ccd area.
+
+setMrel.output needs to do a few things differently from setMrel:
+
+- set measure.Mcal (skipped in setMrel.basic)
+
- set average.Mrel if N < TOO_FEW (not STAR_BAD) (optional!)
+
- use MAX (stats.error, stats.sigma) (optionally)
+
- allow STAR_BAD?
+
+
+ imphotset
+
+imphotset allows you to set certain phot.image table entries. here
+are the options:
+
+imphotset [-photcode code] [-name foo] [-trange (start) (stop)] -flag and value
Index: /branches/elixir/Ohana/src/relphot/doc/notes.txt
===================================================================
--- /branches/elixir/Ohana/src/relphot/doc/notes.txt (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/doc/notes.txt (revision 2485)
@@ -0,0 +1,88 @@
+
+- select different calibration options
+- plot various terms
+
+- rationalize Mmos & Mcal (sum(Mcal) -> Mmos)
+
+- options:
+
+ -update: update phot database
+ -mosaic: apply per-mosaic corrections
+ -grid: apply X,Y-based offset
+ -chips: apply per-chip corrections
+ -trange: select images in time range
+
+
+find_images: selects images in region, filtered by photcode, time range
+
+bcatalog: selects stars & measurments, filtered by mag limit, photcode, time range
+
+image: a single ccd
+mosaic: a collection of associate ccds (mosaiclink[image])
+grid: x,y magnitude correction grid across mosaic
+fit: 2-d fit of mag offsets across field
+
+
+
+m = Mrel + Mcal + Mmos + Mxy + Mgrid
+
+Mrel (average)
+Mcal (image)
+Mmos (mosaic)
+Mxy (Xi, Yi) (Xi,Yi per image ??)
+Mgrid (Xm, Ym)
+
+
+ measure -> average catalog[i].measure[j].averef
+ average -> measure catalog[i].average[j].offset, catalog[i].average[j].Nm
+
+ measure -> image image: image[cat, meas]
+ image -> measure image: catlist[image], measlist[image]
+
+ measure -> mosaic mosaic: mosaic[cat, meas]
+ mosaic -> measure mosaic: catlist[mosaic], measlist[mosaic]
+
+ measure -> gridbin grid: bin[cat, meas]
+ gridbin -> measure grid: catlist[bin], measlist[bin]
+
+ measure -> Xi, Yi -> Mxy
+
+
+getMmos:
+
+ for (i = 0; i < Nmosaic; i++) {
+
+ for (j = 0; j < Nmeasure[i]; j++) {
+ cat = catlist[i][j]
+ meas = measlist[i][j]
+
+ M = catalog[cat].measure[meas].M
+ Mcal = applyMcal (cat, meas)
+ Mrel = applyMrel (cat, meas)
+ Mxy = applyMxy (cat, meas)
+ Mgrid = applyMgrid (cat, meas)
+
+ list[n] = M - Mcal - Mrel - Mxy - Mgrid
+ }
+ }
+
+
+
+excluding data points:
+
+star: N < 2
+image:
+
+
+--
+
+MAG_LIM : keep stars brighter than MAG_LIM (Mcat)
+SIGMA_LIM : keep stars with error < SIGMA_LIM
+STAR_SCATTER : exclude stars with dM > STAR_SCATTER
+IMAGE_SCATTER, NIMAGE_SCATTER : define image scatter limits
+
+
+
+load_images - exclude images
+
+
Index: /branches/elixir/Ohana/src/relphot/include/relphot.h
===================================================================
--- /branches/elixir/Ohana/src/relphot/include/relphot.h (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/include/relphot.h (revision 2485)
@@ -0,0 +1,158 @@
+# include
+# include
+
+/* # define GRID_V1 */
+# define GRID_V2
+
+static union { unsigned char c[4]; float f; } f_undef = { 0xff, 0xff, 0xff, 0xfe };
+# define fUNDEF (f_undef.f)
+
+static union { unsigned char c[8]; float d; } d_undef = { 0xff, 0xff, 0xff, 0xfe, 0x00, 0x00, 0x00, 0x00 };
+# define dUNDEF (d_undef.d)
+
+typedef struct {
+ double xmin, xmax, ymin, ymax;
+ int style, ptype, ltype, etype, color;
+ double lweight, size;
+} Graphdata;
+
+typedef struct {
+ unsigned int start;
+ unsigned int stop;
+ short int Mcal;
+ short int dMcal;
+ short int Xm;
+ short int secz;
+ char code;
+ Coords coords;
+} Mosaic;
+
+typedef struct {
+ double median;
+ double mean;
+ double sigma;
+ double error;
+ double chisq;
+ double min;
+ double max;
+ double total;
+ int Nmeas;
+} StatType;
+
+/* global variables set in parameter file */
+char ImageCat[256];
+char ImageTemplate[256];
+char CatTemplate[256];
+char GSCFILE[256];
+char CATDIR[256];
+char CameraConfig[256];
+
+double MAG_LIM;
+double SIGMA_LIM;
+double IMAGE_SCATTER;
+double IMAGE_OFFSET;
+double STAR_SCATTER;
+double STAR_CHISQ;
+double MIN_ERROR;
+
+int VERBOSE;
+
+int NLOOP;
+int RESET;
+int UPDATE;
+int PLOTSTUFF;
+int SAVEPLOT;
+int SHOW_PARAMS;
+char MOSAICNAME[256];
+char STATMODE[32];
+int STAR_BAD;
+int MEAS_BAD;
+int STAR_TOOFEW;
+int IMAGE_TOOFEW;
+double IMAGE_GOOD_FRACTION;
+int IMAGE_BAD;
+int FREEZE_IMAGES;
+int USE_GRID;
+int PLOTDELAY;
+
+int RELPHOT_GRID_X;
+int RELPHOT_GRID_Y;
+int RELPHOT_GRID_BINNING;
+
+PhotCode *photcode;
+int PhotPrimary;
+int PhotNsec;
+int PhotSec;
+
+int AreaSelect;
+double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
+
+int ImagSelect, ImagMin, ImagMax;
+
+double PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
+enum {black, white, red, orange, yellow, green, blue, indigo, violet};
+
+int TimeSelect;
+unsigned int TSTART, TSTOP;
+
+/***** prototypes ****/
+int main (int argc, char **argv);
+Catalog *load_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion);
+GSCRegion *load_images (char *seed, int *nregion, GSCRegion *fullregion);
+GSCRegion *name_region (char *name, int *Nregions);
+GSCRegion *find_regions (Image *image, int Nimage, int *Nregions, GSCRegion *fullregion);
+GSCRegion *get_regions (double minRa, double maxRa, double minDec, double maxDec, int *Nregions);
+Image *find_images (GSCRegion *region, int Nregion, int *Nimages, int **LineNum);
+double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3);
+int edge_check (double *x1, double *y1, double *x2, double *y2);
+int corner_check (double *x1, double *y1, double *x2, double *y2);
+int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string);
+int bcatalog (Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion);
+int findMosaics (Catalog *catalog, int Ncatalog);
+int setMmos (Catalog *catalog, int Poor);
+int setMrel (Catalog *catalog, int Ncatalog);
+int setrefcode (Image *image, int Nimage);
+int args (int argc, char **argv);
+short getMgrid (int meas, int cat);
+short getMcal (int meas, int cat);
+short getMmos (int meas, int cat);
+short getMrel (Catalog *catalog, int meas, int cat);
+void initGrid (int dX, int dY);
+void initGridBins (Catalog *catalog, int Ncatalog);
+void initMosaics (Image *image, int Nimage);
+void initMosaicGrid (Image *image, int Nimage);
+void initMosaicBins (Catalog *catalog, int Ncatalog);
+void initMrel (Catalog *catalog, int Ncatalog);
+void initialize (int argc, char **argv);
+void ConfigInit (int *argc, char **argv);
+void mark_images (Image *image, int Nimage, Image *timage, int Ntimage);
+void setMgrid (Catalog *catalog);
+void dumpGrid ();
+void initImages (Image *input, int N);
+void initImageBins (Catalog *catalog, int Ncatalog);
+void findImages (Catalog *catalog, int Ncatalog);
+void matchImage (Catalog *catalog, int meas, int cat);
+void matchMosaics (Catalog *catalog, int meas, int cat);
+void setMcal (Catalog *catalog, int Poor);
+void clean_images ();
+void clean_stars (Catalog *catalog, int Ncatalog);
+void reload_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion);
+void free_catalogs (Catalog *catalog, int Ncatalog);
+
+# ifdef GRID_V1
+int setGridMeasure (int meas, int cat, double X, double Y);
+# endif
+
+# ifdef GRID_V2
+int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum);
+# endif
+
+void set_ZP (double ZERO);
+
+char *GetPhotnamebyCode (PhotCodeData *photcodes, int code);
+StatType statsStarS (Catalog *catalog, int Ncatalog);
+StatType statsStarX (Catalog *catalog, int Ncatalog);
+StatType statsStarN (Catalog *catalog, int Ncatalog);
+StatType statsImageN (Catalog *catalog);
+StatType statsImageX (Catalog *catalog);
+StatType statsImageM (Catalog *catalog);
Index: /branches/elixir/Ohana/src/relphot/src/ConfigInit.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/ConfigInit.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/ConfigInit.c (revision 2485)
@@ -0,0 +1,70 @@
+# include "relphot.h"
+void GetConfig (char *config, char *field, char *format, int N, void *ptr);
+
+void ConfigInit (int *argc, char **argv) {
+
+ double ZERO_POINT;
+ char *config, *file;
+ char PhotCodeFile[256];
+
+ /*** load configuration info ***/
+ file = SelectConfigFile (argc, argv, "ptolemy");
+ config = LoadConfigFile (file);
+ if (config == (char *) NULL) {
+ fprintf (stderr, "ERROR: can't find configuration file %s\n", file);
+ if (file != (char *) NULL) free (file);
+ exit (0);
+ }
+ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
+
+ GetConfig (config, "MAG_LIM", "%lf", 0, &MAG_LIM);
+ GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM);
+ GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER);
+ GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER);
+ GetConfig (config, "IMAGE_OFFSET", "%lf", 0, &IMAGE_OFFSET);
+ MAG_LIM *= 1000;
+ SIGMA_LIM *= 1000;
+ STAR_SCATTER *= 1000;
+ IMAGE_SCATTER *= 1000;
+ IMAGE_OFFSET *= 1000;
+
+ GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ);
+ GetConfig (config, "STAR_TOOFEW", "%d", 0, &STAR_TOOFEW);
+ GetConfig (config, "IMAGE_TOOFEW", "%d", 0, &IMAGE_TOOFEW);
+ GetConfig (config, "IMAGE_GOOD_FRACTION", "%lf", 0, &IMAGE_GOOD_FRACTION);
+
+ GetConfig (config, "IMAGE_CATALOG", "%s", 0, ImageCat);
+ GetConfig (config, "IMAGE_CATALOG_TEMPLATE", "%s", 0, ImageTemplate);
+ GetConfig (config, "CATALOG_TEMPLATE", "%s", 0, CatTemplate);
+ GetConfig (config, "GSCFILE", "%s", 0, GSCFILE);
+ GetConfig (config, "CATDIR", "%s", 0, CATDIR);
+ GetConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile);
+ GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT);
+
+ GetConfig (config, "RELPHOT_GRID_X", "%d", 0, &RELPHOT_GRID_X);
+ GetConfig (config, "RELPHOT_GRID_Y", "%d", 0, &RELPHOT_GRID_Y);
+ GetConfig (config, "RELPHOT_GRID_BINNING", "%d", 0, &RELPHOT_GRID_BINNING);
+ GetConfig (config, "CAMERA_CONFIG", "%s", 0, CameraConfig);
+
+ if (!LoadPhotcodes (PhotCodeFile)) {
+ fprintf (stderr, "error loading photcodes\n");
+ exit (1);
+ }
+ SetZeroPoint (ZERO_POINT);
+
+ free (config);
+ free (file);
+
+}
+
+void GetConfig (char *config, char *field, char *format, int N, void *ptr) {
+
+ int status;
+
+ status = ScanConfig (config, field, format, N, ptr);
+ if (!status) {
+ fprintf (stderr, "error in config, cannot find %s\n", field);
+ exit (1);
+ }
+ return;
+}
Index: /branches/elixir/Ohana/src/relphot/src/GridOps.v1.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/GridOps.v1.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/GridOps.v1.c (revision 2485)
@@ -0,0 +1,427 @@
+# include "relphot.h"
+
+static int Ngrid;
+static short *gridM;
+static short *gridS;
+static int *gridN;
+static int gridX;
+static int gridY;
+
+static int **bin;
+static int **Xmeas;
+static int **Ymeas;
+
+static int **clist;
+static int **mlist;
+static int *Nlist;
+static int *NLIST;
+
+void initGrid (int dX, int dY) {
+
+ int i;
+
+ /* define mosaic 2d correction grid */
+ gridX = dX / RELPHOT_GRID_BINNING + 1;
+ gridY = dY / RELPHOT_GRID_BINNING + 1;
+ Ngrid = gridX * gridY;
+
+ ALLOCATE (gridM, short, Ngrid);
+ ALLOCATE (gridS, short, Ngrid);
+ ALLOCATE (gridN, int, Ngrid);
+ bzero (gridM, Ngrid*sizeof(short));
+ bzero (gridS, Ngrid*sizeof(short));
+ bzero (gridN, Ngrid*sizeof(int));
+
+}
+
+void initGridBins (Catalog *catalog, int Ncatalog) {
+
+ int i, j;
+
+ if (!USE_GRID) return;
+
+ /* define cat,meas -> grid pointers */
+ ALLOCATE (bin, int *, Ncatalog);
+ ALLOCATE (Xmeas, int *, Ncatalog);
+ ALLOCATE (Ymeas, int *, Ncatalog);
+ for (i = 0; i < Ncatalog; i++) {
+ ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+ ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
+ ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
+ for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
+ }
+
+ /* define grid -> cat,meas pointers */
+ ALLOCATE (Nlist, int, Ngrid);
+ ALLOCATE (NLIST, int, Ngrid);
+ ALLOCATE (clist, int *, Ngrid);
+ ALLOCATE (mlist, int *, Ngrid);
+
+ for (i = 0; i < Ngrid; i++) {
+ Nlist[i] = 0;
+ NLIST[i] = 100;
+ ALLOCATE (clist[i], int, NLIST[i]);
+ ALLOCATE (mlist[i], int, NLIST[i]);
+ }
+}
+
+void freeGridBins (int Ncatalog) {
+
+ int i;
+
+ if (!USE_GRID) return;
+
+ /* define cat,meas -> grid pointers */
+ for (i = 0; i < Ncatalog; i++) {
+ free (bin[i]);
+ free (Xmeas[i]);
+ free (Ymeas[i]);
+ }
+ free (bin);
+ free (Xmeas);
+ free (Ymeas);
+
+ /* define grid -> cat,meas pointers */
+ for (i = 0; i < Ngrid; i++) {
+ free (clist[i]);
+ free (mlist[i]);
+ }
+ free (Nlist);
+ free (NLIST);
+ free (clist);
+ free (mlist);
+}
+
+int setGridMeasure (int meas, int cat, double X, double Y) {
+
+ int ix, iy, i;
+
+ ix = X / RELPHOT_GRID_BINNING;
+ iy = Y / RELPHOT_GRID_BINNING;
+ if (ix < 0) goto escape;
+ if (iy < 0) goto escape;
+ if (ix >= gridX) goto escape;
+ if (iy >= gridY) goto escape;
+
+ i = ix + iy*gridX;
+
+ bin[cat][meas] = i;
+ Xmeas[cat][meas] = X;
+ Ymeas[cat][meas] = Y;
+ clist[i][Nlist[i]] = cat;
+ mlist[i][Nlist[i]] = meas;
+
+ Nlist[i] ++;
+ if (Nlist[i] == NLIST[i]) {
+ NLIST[i] += 100;
+ REALLOCATE (clist[i], int, NLIST[i]);
+ REALLOCATE (mlist[i], int, NLIST[i]);
+ }
+ return (TRUE);
+
+ escape:
+ fprintf (stderr, "error: star out of grid\n");
+ exit (1);
+}
+
+short getMgrid (int meas, int cat) {
+
+ int i;
+ short value;
+
+ if (!USE_GRID) return (0);
+ i = bin[cat][meas];
+ if (i == -1) return (NO_MAG);
+
+ value = gridM[i];
+ return (value);
+}
+
+/* determine Mgrid values for all grid bins */
+void setMgrid (Catalog *catalog) {
+
+ int i, j, m, c, n, N, Nmax;
+ double *list, *dlist;
+ short Msys, Mrel, Mcal, Mmos;
+ StatType stats;
+
+ if (!USE_GRID) return;
+
+ Nmax = Nlist[0];
+ for (i = 0; i < Ngrid; i++) {
+ Nmax = MAX (Nmax, Nlist[i]);
+ }
+ ALLOCATE (list, double, Nmax);
+ ALLOCATE (dlist, double, Nmax);
+
+ for (i = 0; i < Ngrid; i++) {
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, c)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+ list[N] = Msys - Mrel - Mcal - Mmos;
+ dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+
+ liststats (list, dlist, N, &stats);
+ gridM[i] = stats.mean;
+ gridS[i] = stats.sigma;
+ gridN[i] = N;
+ }
+ free (list);
+ free (dlist);
+}
+
+void plot_grid (Catalog *catalog) {
+
+ int i, j, m, c, n, N;
+ short Msys, Mrel, Mcal, Mmos;
+ double *xlist, *Mlist, *dlist, *ylist;
+ Graphdata graphdata;
+
+ if (!USE_GRID) return;
+
+ N = 0;
+ for (i = 0; i < Ngrid; i++)
+ N += Nlist[i];
+
+ ALLOCATE (xlist, double, N);
+ ALLOCATE (ylist, double, N);
+ ALLOCATE (Mlist, double, N);
+ ALLOCATE (dlist, double, N);
+
+ N = 0;
+ for (i = 0; i < Ngrid; i++) {
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, c)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+
+ xlist[N] = Xmeas[c][m];
+ ylist[N] = Ymeas[c][m];
+ Mlist[N] = Msys - Mrel - Mcal - Mmos;
+ dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i];
+ Mlist[N] *= 0.001;
+ dlist[N] *= 0.001;
+ N++;
+ }
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "XdM.png");
+ plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
+ plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
+
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
+
+ free (ylist);
+ free (xlist);
+ free (Mlist);
+ free (dlist);
+
+}
+
+int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage);
+Image *getimage (int N);
+
+void dump_grid () {
+
+ int i, j, Nx, Ny, Nimage;
+ int *imlist;
+ char *p, *pname, ccdname[80];
+ FILE *f;
+ Header header, theader;
+ Matrix matrix;
+ Mosaic *refmosaic;
+ Image *image;
+
+ /* select reference mosaic image */
+ imlist = SelectRefMosaic (&refmosaic, &Nimage);
+
+ /* we are writing to this file */
+ f = fopen ("mosaic.fits", "w");
+ if (f == (FILE *) NULL) {
+ fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
+ return;
+ }
+
+ /* create empty phu */
+ fits_init_header (&header);
+ header.extend = TRUE;
+ fits_create_header (&header);
+ fits_create_matrix (&header, &matrix);
+ fits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
+ fits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
+ fits_fwrite_header (f, &header);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid mag values */
+ fits_init_header (&theader);
+ theader.Naxes = 2;
+ theader.Naxis[0] = gridX;
+ theader.Naxis[1] = gridY;
+ theader.bitpix = -32;
+ fits_create_Theader (&theader, "IMAGE");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid Nmeas values */
+ fits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid sigma values */
+ fits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* calculate value for each CCD pixel, write out CCD images */
+ for (i = 0; i < Nimage; i++) {
+ image = getimage (imlist[i]);
+ pname = GetPhotcodeNamebyCode (image[0].source);
+
+ /* this is kind of bogus... */
+ /* pname is CAMERA.FILTER.CCD, grab the CCD */
+ p = strrchr (pname, '.');
+ if (p == (char *) NULL) {
+ fprintf (stderr, "error parsing photcode %s\n", pname);
+ exit (2);
+ }
+ p++;
+ sprintf (ccdname, "ccd%s", p);
+
+ fits_modify (&theader, "EXTNAME", "%s", 1, ccdname);
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&theader, "PHOTCODE", "%s", 1, pname);
+ fits_modify (&theader, "NX", "%d", 1, image[i].NX);
+ fits_modify (&theader, "NY", "%d", 1, image[i].NY);
+ write_coords (&theader, &image[0].coords);
+
+ Nx = 2 * image[i].NX / RELPHOT_GRID_BINNING;
+ Ny = 2 * image[i].NY / RELPHOT_GRID_BINNING;
+ theader.Naxis[0] = Nx;
+ theader.Naxis[1] = Ny;
+ fits_modify (&theader, "NAXIS1", "%d", 1, Nx);
+ fits_modify (&theader, "NAXIS2", "%d", 1, Ny);
+ fits_create_matrix (&theader, &matrix);
+
+ InterpolateGrid ((float *)matrix.buffer, Nx, Ny, &image[0].coords, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+ }
+}
+
+InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
+
+ int i, j;
+ double x, y, r, d, X, Y, dx, dy, dVx, dVy, dVxy;
+ double V00, V01, V10, V11;
+ double wV00, wV01, wV10, wV11;
+ double dV00, dV01, dV10, dV11;
+ double v1, v2, value;
+ int ix, iy, N;
+
+ for (i = 0; i < Nx; i++) {
+ for (j = 0; j < Ny; j++) {
+ x = i * RELPHOT_GRID_BINNING / 2;
+ y = j * RELPHOT_GRID_BINNING / 2;
+ XY_to_RD (&r, &d, x, y, ccd);
+ RD_to_XY (&X, &Y, r, d, gcoords);
+
+ X = X / RELPHOT_GRID_BINNING;
+ Y = Y / RELPHOT_GRID_BINNING;
+
+ ix = (int) X;
+ dx = X - ix;
+ iy = (int) Y;
+ dy = Y - iy;
+
+ if (ix < 0) continue;
+ if (iy < 0) continue;
+ if (ix >= gridX) continue;
+ if (iy >= gridY) continue;
+
+ N = ix + iy*gridX;
+ V00 = gridM[N];
+ V10 = gridM[N + 1];
+ V01 = gridM[N + gridX];
+ V11 = gridM[N + gridX + 1];
+
+ dV00 = gridS[N];
+ dV10 = gridS[N + 1];
+ dV01 = gridS[N + gridX];
+ dV11 = gridS[N + gridX + 1];
+
+ wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
+ wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
+ wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
+ wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);
+
+ v1 = wV00*V00*(1 + dx*dy - dx - dy) +
+ wV10*V10*(dx - dx*dy) +
+ wV01*V01*(dy - dx*dy) +
+ wV11*V11*(dx*dy);
+
+ v2 = wV00*(1 + dx*dy - dx - dy) +
+ wV10*(dx - dx*dy) +
+ wV01*(dy - dx*dy) +
+ wV11*(dx*dy);
+
+ value = v1 / v2;
+ buffer[j*Nx + i] = value;
+ }
+ }
+}
Index: /branches/elixir/Ohana/src/relphot/src/GridOps.v2.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/GridOps.v2.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/GridOps.v2.c (revision 2485)
@@ -0,0 +1,579 @@
+# include "relphot.h"
+
+static int Ngrid;
+static short *gridM;
+static short *gridS;
+static int *gridN;
+static int gridX;
+static int gridY;
+
+static int **bin;
+static int **Xmeas;
+static int **Ymeas;
+
+static int **clist;
+static int **mlist;
+static int *Nlist;
+static int *NLIST;
+
+static struct {
+ int Nchip;
+ int Mx, My; /* mosaic size in chips */
+ int Nx, Ny; /* chip size in pixels */
+ int *Fx, *Fy; /* chip flip */
+ int *Ox, *Oy; /* chip offset */
+ char **ccdname;
+} camera;
+
+# ifdef GRID_V1
+void initGrid (int dX, int dY) {
+
+ int i;
+
+ /* define mosaic 2d correction grid */
+ gridX = dX / RELPHOT_GRID_BINNING + 1;
+ gridY = dY / RELPHOT_GRID_BINNING + 1;
+ Ngrid = gridX * gridY;
+
+ ALLOCATE (gridM, short, Ngrid);
+ ALLOCATE (gridS, short, Ngrid);
+ ALLOCATE (gridN, int, Ngrid);
+ bzero (gridM, Ngrid*sizeof(short));
+ bzero (gridS, Ngrid*sizeof(short));
+ bzero (gridN, Ngrid*sizeof(int));
+
+}
+# endif
+
+# ifdef GRID_V2
+void initGrid (int dX, int dY) {
+
+ int i;
+ char field[64], *config, line[256];
+
+ /* load camera config file */
+ config = LoadConfigFile (CameraConfig);
+ if (config == (char *) NULL) {
+ fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig);
+ exit (1);
+ }
+
+ /* load basic mosaic parameters */
+ ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip);
+ ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx);
+ ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My);
+ ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
+ ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
+
+ ALLOCATE (camera.Fx, int, camera.Nchip);
+ ALLOCATE (camera.Fy, int, camera.Nchip);
+ ALLOCATE (camera.Ox, int, camera.Nchip);
+ ALLOCATE (camera.Oy, int, camera.Nchip);
+ ALLOCATE (camera.ccdname, char *, camera.Nchip);
+
+ /* load per-chip parameters */
+ for (i = 0; i < camera.Nchip; i++) {
+ ALLOCATE (camera.ccdname[i], char, 256);
+ sprintf (field, "CCD.%d", i);
+ ScanConfig (config, field, "%s", 1, line);
+ sscanf (line, "%s %d %d %d %d", camera.ccdname[i], &camera.Ox[i], &camera.Oy[i], &camera.Fx[i], &camera.Fy[i]);
+ }
+ free (config);
+
+ /* define mosaic 2d correction grid:
+ * GRID_X is the number of grid pixels per chip in the X direction */
+ gridX = RELPHOT_GRID_X * camera.Mx;
+ gridY = RELPHOT_GRID_Y * camera.My;
+ Ngrid = gridX * gridY;
+
+ ALLOCATE (gridM, short, Ngrid);
+ ALLOCATE (gridS, short, Ngrid);
+ ALLOCATE (gridN, int, Ngrid);
+ bzero (gridM, Ngrid*sizeof(short));
+ bzero (gridS, Ngrid*sizeof(short));
+ bzero (gridN, Ngrid*sizeof(int));
+
+}
+# endif
+
+void initGridBins (Catalog *catalog, int Ncatalog) {
+
+ int i, j;
+
+ if (!USE_GRID) return;
+
+ /* define cat,meas -> grid pointers */
+ ALLOCATE (bin, int *, Ncatalog);
+ ALLOCATE (Xmeas, int *, Ncatalog);
+ ALLOCATE (Ymeas, int *, Ncatalog);
+ for (i = 0; i < Ncatalog; i++) {
+ ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+ ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
+ ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
+ for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
+ }
+
+ /* define grid -> cat,meas pointers */
+ ALLOCATE (Nlist, int, Ngrid);
+ ALLOCATE (NLIST, int, Ngrid);
+ ALLOCATE (clist, int *, Ngrid);
+ ALLOCATE (mlist, int *, Ngrid);
+
+ for (i = 0; i < Ngrid; i++) {
+ Nlist[i] = 0;
+ NLIST[i] = 100;
+ ALLOCATE (clist[i], int, NLIST[i]);
+ ALLOCATE (mlist[i], int, NLIST[i]);
+ }
+}
+
+void freeGridBins (int Ncatalog) {
+
+ int i;
+
+ if (!USE_GRID) return;
+
+ /* define cat,meas -> grid pointers */
+ for (i = 0; i < Ncatalog; i++) {
+ free (bin[i]);
+ free (Xmeas[i]);
+ free (Ymeas[i]);
+ }
+ free (bin);
+ free (Xmeas);
+ free (Ymeas);
+
+ /* define grid -> cat,meas pointers */
+ for (i = 0; i < Ngrid; i++) {
+ free (clist[i]);
+ free (mlist[i]);
+ }
+ free (Nlist);
+ free (NLIST);
+ free (clist);
+ free (mlist);
+}
+
+# ifdef GRID_V1
+int setGridMeasure (int meas, int cat, double X, double Y) {
+
+ int ix, iy, i;
+
+ ix = X / RELPHOT_GRID_BINNING;
+ iy = Y / RELPHOT_GRID_BINNING;
+ if (ix < 0) goto escape;
+ if (iy < 0) goto escape;
+ if (ix >= gridX) goto escape;
+ if (iy >= gridY) goto escape;
+
+ i = ix + iy*gridX;
+
+ bin[cat][meas] = i;
+ Xmeas[cat][meas] = X;
+ Ymeas[cat][meas] = Y;
+ clist[i][Nlist[i]] = cat;
+ mlist[i][Nlist[i]] = meas;
+
+ Nlist[i] ++;
+ if (Nlist[i] == NLIST[i]) {
+ NLIST[i] += 100;
+ REALLOCATE (clist[i], int, NLIST[i]);
+ REALLOCATE (mlist[i], int, NLIST[i]);
+ }
+ return (TRUE);
+
+escape:
+ fprintf (stderr, "error: star out of grid\n");
+ exit (1);
+}
+# endif
+
+# ifdef GRID_V2
+int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum) {
+
+ int ix, iy, Cx, Cy, i;
+ double x, y;
+
+ /* X, Y are chip coords on chip ccdnum */
+
+ /* normalize X & Y */
+ x = X;
+ if (camera.Fx[ccdnum]) x = camera.Nx - X;
+ y = Y;
+ if (camera.Fy[ccdnum]) y = camera.Ny - Y;
+
+ /* grid coords on the chip */
+ Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
+ Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
+
+ /* coordinates in the grid */
+ ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
+ iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
+
+ i = ix + iy*gridX;
+
+ bin[cat][meas] = i;
+ Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
+ Ymeas[cat][meas] = Y + camera.Oy[ccdnum]*camera.Ny;
+ clist[i][Nlist[i]] = cat;
+ mlist[i][Nlist[i]] = meas;
+
+ Nlist[i] ++;
+ if (Nlist[i] == NLIST[i]) {
+ NLIST[i] += 100;
+ REALLOCATE (clist[i], int, NLIST[i]);
+ REALLOCATE (mlist[i], int, NLIST[i]);
+ }
+ return (TRUE);
+
+escape:
+ fprintf (stderr, "error: star out of grid\n");
+ exit (1);
+}
+# endif
+
+short getMgrid (int meas, int cat) {
+
+ int i;
+ short value;
+
+ if (!USE_GRID) return (0);
+ i = bin[cat][meas];
+ if (i == -1) return (NO_MAG);
+
+ value = gridM[i];
+ return (value);
+}
+
+/* determine Mgrid values for all grid bins */
+void setMgrid (Catalog *catalog) {
+
+ int i, j, m, c, n, N, Nmax;
+ double *list, *dlist;
+ short Msys, Mrel, Mcal, Mmos;
+ StatType stats;
+
+ if (!USE_GRID) return;
+
+ Nmax = Nlist[0];
+ for (i = 0; i < Ngrid; i++) {
+ Nmax = MAX (Nmax, Nlist[i]);
+ }
+ ALLOCATE (list, double, Nmax);
+ ALLOCATE (dlist, double, Nmax);
+
+ for (i = 0; i < Ngrid; i++) {
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, c)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+ list[N] = Msys - Mrel - Mcal - Mmos;
+ dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+
+ liststats (list, dlist, N, &stats);
+ gridM[i] = stats.mean;
+ gridS[i] = stats.sigma;
+ gridN[i] = N;
+ }
+ free (list);
+ free (dlist);
+}
+
+void plot_grid (Catalog *catalog) {
+
+ int i, j, m, c, n, N;
+ short Msys, Mrel, Mcal, Mmos;
+ double *xlist, *Mlist, *dlist, *ylist;
+ Graphdata graphdata;
+
+ if (!USE_GRID) return;
+
+ N = 0;
+ for (i = 0; i < Ngrid; i++)
+ N += Nlist[i];
+
+ ALLOCATE (xlist, double, N);
+ ALLOCATE (ylist, double, N);
+ ALLOCATE (Mlist, double, N);
+ ALLOCATE (dlist, double, N);
+
+ N = 0;
+ for (i = 0; i < Ngrid; i++) {
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, c)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+
+ xlist[N] = Xmeas[c][m];
+ ylist[N] = Ymeas[c][m];
+ Mlist[N] = Msys - Mrel - Mcal - Mmos;
+ dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i];
+ Mlist[N] *= 0.001;
+ dlist[N] *= 0.001;
+ N++;
+ }
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "XdM.png");
+ plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
+ plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
+
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
+
+ free (ylist);
+ free (xlist);
+ free (Mlist);
+ free (dlist);
+
+}
+
+int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage);
+Image *getimage (int N);
+
+void dump_grid () {
+
+ int i, j, Nx, Ny, Nimage;
+ int *imlist;
+ char *p, *pname, ccdname[80];
+ FILE *f;
+ Header header, theader;
+ Matrix matrix;
+ Mosaic *refmosaic;
+ Image *image;
+
+ /* select reference mosaic image */
+ imlist = SelectRefMosaic (&refmosaic, &Nimage);
+
+ /* we are writing to this file */
+ f = fopen ("mosaic.fits", "w");
+ if (f == (FILE *) NULL) {
+ fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
+ return;
+ }
+
+ /* create empty phu */
+ fits_init_header (&header);
+ header.extend = TRUE;
+ fits_create_header (&header);
+ fits_create_matrix (&header, &matrix);
+ fits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
+ fits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
+ fits_fwrite_header (f, &header);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid mag values */
+ fits_init_header (&theader);
+ theader.Naxes = 2;
+ theader.Naxis[0] = gridX;
+ theader.Naxis[1] = gridY;
+ theader.bitpix = -32;
+ fits_create_Theader (&theader, "IMAGE");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid Nmeas values */
+ fits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+ /* save grid sigma values */
+ fits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_create_matrix (&theader, &matrix);
+ for (i = 0; i < gridX; i++) {
+ for (j = 0; j < gridY; j++) {
+ fits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]);
+ }
+ }
+ write_coords (&theader, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+
+# ifdef GRID_V1
+ /* calculate pixel values for each CCD pixel, write out CCD images */
+ /* grid pixels are in RA,DEC coords, transform to image and interpolate */
+ for (i = 0; i < Nimage; i++) {
+ image = getimage (imlist[i]);
+ pname = GetPhotcodeNamebyCode (image[0].source);
+
+ /* this is kind of bogus... */
+ /* pname is CAMERA.FILTER.CCD, grab the CCD */
+ p = strrchr (pname, '.');
+ if (p == (char *) NULL) {
+ fprintf (stderr, "error parsing photcode %s\n", pname);
+ exit (2);
+ }
+ p++;
+ sprintf (ccdname, "ccd%s", p);
+
+ fits_modify (&theader, "EXTNAME", "%s", 1, ccdname);
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&theader, "PHOTCODE", "%s", 1, pname);
+ fits_modify (&theader, "NX", "%d", 1, image[i].NX);
+ fits_modify (&theader, "NY", "%d", 1, image[i].NY);
+ write_coords (&theader, &image[0].coords);
+
+ Nx = 2 * image[i].NX / RELPHOT_GRID_BINNING;
+ Ny = 2 * image[i].NY / RELPHOT_GRID_BINNING;
+ theader.Naxis[0] = Nx;
+ theader.Naxis[1] = Ny;
+ fits_modify (&theader, "NAXIS1", "%d", 1, Nx);
+ fits_modify (&theader, "NAXIS2", "%d", 1, Ny);
+ fits_create_matrix (&theader, &matrix);
+
+ InterpolateGrid ((float *)matrix.buffer, Nx, Ny, &image[0].coords, &refmosaic[0].coords);
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+ }
+# endif
+
+# ifdef GRID_V2
+ /* calculate value for each CCD pixel, write out CCD images */
+ /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
+ for (i = 0; i < camera.Nchip; i++) {
+ int ix, iy, x, y, X, Y, bin;
+
+ fits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[i]);
+ fits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+ fits_modify (&theader, "NX", "%d", 1, camera.Nx);
+ fits_modify (&theader, "NY", "%d", 1, camera.Ny);
+
+ theader.Naxis[0] = RELPHOT_GRID_X;
+ theader.Naxis[1] = RELPHOT_GRID_Y;
+ fits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X);
+ fits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y);
+ fits_create_matrix (&theader, &matrix);
+
+ for (Y = 0; Y < RELPHOT_GRID_Y; Y++) {
+ for (X = 0; X < RELPHOT_GRID_X; X++) {
+
+ /* normalize X & Y */
+ x = X;
+ if (camera.Fx[i]) x = RELPHOT_GRID_X - X - 1;
+ y = Y;
+ if (camera.Fy[i]) y = RELPHOT_GRID_Y - Y - 1;
+
+ /* coordinates in the grid */
+ ix = x + camera.Ox[i]*RELPHOT_GRID_X;
+ iy = y + camera.Oy[i]*RELPHOT_GRID_Y;
+
+ bin = ix + iy*gridX;
+ fits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]);
+ }
+ }
+ fits_fwrite_header (f, &theader);
+ fits_fwrite_matrix (f, &matrix);
+ fits_free_matrix (&matrix);
+ }
+# endif
+
+}
+
+InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
+
+ int i, j;
+ double x, y, r, d, X, Y, dx, dy, dVx, dVy, dVxy;
+ double V00, V01, V10, V11;
+ double wV00, wV01, wV10, wV11;
+ double dV00, dV01, dV10, dV11;
+ double v1, v2, value;
+ int ix, iy, N;
+
+ for (i = 0; i < Nx; i++) {
+ for (j = 0; j < Ny; j++) {
+ x = i * RELPHOT_GRID_BINNING / 2;
+ y = j * RELPHOT_GRID_BINNING / 2;
+ XY_to_RD (&r, &d, x, y, ccd);
+ RD_to_XY (&X, &Y, r, d, gcoords);
+
+ X = X / RELPHOT_GRID_BINNING;
+ Y = Y / RELPHOT_GRID_BINNING;
+
+ ix = (int) X;
+ dx = X - ix;
+ iy = (int) Y;
+ dy = Y - iy;
+
+ if (ix < 0) continue;
+ if (iy < 0) continue;
+ if (ix >= gridX) continue;
+ if (iy >= gridY) continue;
+
+ N = ix + iy*gridX;
+ V00 = gridM[N];
+ V10 = gridM[N + 1];
+ V01 = gridM[N + gridX];
+ V11 = gridM[N + gridX + 1];
+
+ dV00 = gridS[N];
+ dV10 = gridS[N + 1];
+ dV01 = gridS[N + gridX];
+ dV11 = gridS[N + gridX + 1];
+
+ wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
+ wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
+ wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
+ wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);
+
+ v1 = wV00*V00*(1 + dx*dy - dx - dy) +
+ wV10*V10*(dx - dx*dy) +
+ wV01*V01*(dy - dx*dy) +
+ wV11*V11*(dx*dy);
+
+ v2 = wV00*(1 + dx*dy - dx - dy) +
+ wV10*(dx - dx*dy) +
+ wV01*(dy - dx*dy) +
+ wV11*(dx*dy);
+
+ value = v1 / v2;
+ buffer[j*Nx + i] = value;
+ }
+ }
+}
Index: /branches/elixir/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/ImageOps.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/ImageOps.c (revision 2485)
@@ -0,0 +1,465 @@
+# include "relphot.h"
+
+static unsigned int *start;
+static unsigned int *stop;
+static int **bin;
+
+static int **clist;
+static int **mlist;
+static int *Nlist;
+static int *NLIST;
+
+static Image *image;
+static int Nimage;
+
+void initImages (Image *input, int N) {
+
+ int i, j;
+
+ image = input;
+ Nimage = N;
+
+ ALLOCATE (start, unsigned, Nimage);
+ ALLOCATE (stop, unsigned, Nimage);
+
+ for (i = 0; i < Nimage; i++) {
+ start[i] = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
+ stop[i] = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
+ }
+}
+
+void initImageBins (Catalog *catalog, int Ncatalog) {
+
+ int i, j;
+
+ ALLOCATE (bin, int *, Ncatalog);
+ for (i = 0; i < Ncatalog; i++) {
+ ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+ for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
+ }
+
+ ALLOCATE (Nlist, int, Nimage);
+ ALLOCATE (NLIST, int, Nimage);
+ ALLOCATE (clist, int *, Nimage);
+ ALLOCATE (mlist, int *, Nimage);
+
+ for (i = 0; i < Nimage; i++) {
+ Nlist[i] = 0;
+ NLIST[i] = 100;
+ ALLOCATE (clist[i], int, NLIST[i]);
+ ALLOCATE (mlist[i], int, NLIST[i]);
+ }
+}
+
+void freeImageBins (int Ncatalog) {
+
+ int i;
+
+ for (i = 0; i < Ncatalog; i++) {
+ free (bin[i]);
+ }
+ free (bin);
+ for (i = 0; i < Nimage; i++) {
+ free (clist[i]);
+ free (mlist[i]);
+ }
+ free (clist);
+ free (mlist);
+}
+
+/* select all image equivalent to the current photcode */
+void findImages (Catalog *catalog, int Ncatalog) {
+
+ int i, j, ecode;
+
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Nmeasure; j++) {
+ ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
+ if (photcode[0].code != ecode) continue;
+ matchImage (catalog, j, i);
+ }
+ }
+}
+
+void matchImage (Catalog *catalog, int meas, int cat) {
+
+ int i, ave, ccdnum;
+ double ra, dec, X, Y;
+ char *pname, *filter, *p, base[256];
+ Measure *measure;
+
+ measure = &catalog[cat].measure[meas];
+ for (i = 0; i < Nimage; i++) {
+ if (image[0].source == -1) continue;
+ if (measure[0].source != image[i].source) continue;
+ if (measure[0].t < start[i]) continue;
+ if (measure[0].t > stop[i]) continue;
+
+# ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
+ if (USE_GRID) {
+
+ /* identify the ccd on the basis of the photcode name */
+ pname = GetPhotcodeNamebyCode (image[i].source);
+ filter = photcode[0].name;
+ sprintf (base, "%s.%s.", MOSAICNAME, filter);
+ if (strncmp (pname, base, strlen (base))) continue;
+ p = pname + strlen(base);
+ if (*p == 0) continue;
+ ccdnum = atoi (p);
+ /* some ambiguity here between seq number and id number */
+
+ /* add this measurement to the grid cell for this chip */
+ ave = measure[0].averef;
+ ra = catalog[cat].average[ave].R - measure[0].dR / 360000.0;
+ dec = catalog[cat].average[ave].D - measure[0].dD / 360000.0;
+
+ /* X,Y always positive-definite in range 0,0 - dX, dY */
+ RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
+ setGridMeasure (meas, cat, X, Y, ccdnum);
+ }
+# endif
+
+ bin[cat][meas] = i;
+
+ clist[i][Nlist[i]] = cat;
+ mlist[i][Nlist[i]] = meas;
+ Nlist[i] ++;
+
+ if (Nlist[i] == NLIST[i]) {
+ NLIST[i] += 100;
+ REALLOCATE (clist[i], int, NLIST[i]);
+ REALLOCATE (mlist[i], int, NLIST[i]);
+ }
+ return;
+ }
+ /* fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].source); */
+}
+
+short getMcal (int meas, int cat) {
+
+ int i;
+ short value;
+
+ i = bin[cat][meas];
+ if (i == -1) return (NO_MAG);
+
+ if (image[i].code & IMAGE_BAD) return (NO_MAG);
+ value = image[i].Mcal;
+ return (value);
+}
+
+Coords *getCoords (int meas, int cat) {
+
+ int i;
+ short value;
+
+ i = bin[cat][meas];
+ if (i == -1) return (NULL);
+ return (&image[i].coords);
+}
+
+/* determine Mcal values for all images */
+void setMcal (Catalog *catalog, int PoorImages) {
+
+ int i, j, m, c, n, N, Nmax, mark, bad;
+ short Msys, Mrel, Mmos, Mgrid;
+ double *list, *dlist;
+ StatType stats;
+
+ if (FREEZE_IMAGES) return;
+
+ if (PoorImages) {
+ IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
+ }
+
+ Nmax = 0;
+ for (i = 0; i < Nimage; i++) {
+ Nmax = MAX (Nmax, Nlist[i]);
+ }
+ ALLOCATE (list, double, Nmax);
+ ALLOCATE (dlist, double, Nmax);
+
+ for (i = 0; i < Nimage; i++) {
+
+ /* on PoorImages run, skip good images */
+ if (PoorImages) {
+ bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
+ if (!bad) continue;
+ }
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mmos = getMmos (m, c)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+ list[N] = Msys - Mrel - Mmos - Mgrid;
+ dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+ /* Nlist[i] is all measurements, N is good measurements */
+
+ /* too few good measurements or too many bad measurements */
+ if (!PoorImages) {
+ mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
+ if (mark) {
+ image[i].code |= ID_IMAGE_FEW;
+ } else {
+ image[i].code &= ~ID_IMAGE_FEW;
+ }
+ }
+
+ liststats (list, dlist, N, &stats);
+ image[i].Mcal = stats.mean;
+ image[i].dMcal = stats.sigma;
+ image[i].Xm = 100.0*log10(stats.chisq);
+ }
+ free (list);
+ free (dlist);
+ if (PoorImages) {
+ IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
+ STAR_BAD = ID_STAR_POOR | ID_STAR_FEW;
+ MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
+ }
+ return;
+}
+
+/* mark image if: abs(Mcal) too large, dMcal too large */
+void clean_images () {
+
+ int i, N, mark, Nmark;
+ double *mlist, *slist, *dlist;
+ double MaxOffset, MaxScatter, MedOffset;
+ StatType stats;
+
+ if (FREEZE_IMAGES) return;
+
+ if (VERBOSE) fprintf (stderr, "marking poor images\n");
+
+ ALLOCATE (mlist, double, Nimage);
+ ALLOCATE (slist, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+
+ for (i = N = 0; i < Nimage; i++) {
+ if (image[i].code & IMAGE_BAD) continue;
+ mlist[N] = fabs (image[i].Mcal);
+ slist[N] = image[i].dMcal;
+ dlist[N] = 1;
+ N++;
+ }
+ initstats ("MEAN");
+ liststats (mlist, dlist, N, &stats);
+ MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
+ MedOffset = stats.median;
+ liststats (slist, dlist, N, &stats);
+ MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
+ fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", 0.001*MedOffset, 0.001*stats.median, 0.001*MaxScatter, 0.001*MaxOffset);
+
+ Nmark = 0;
+ for (i = 0; i < Nimage; i++) {
+ mark = FALSE;
+ image[i].code &= ~ID_IMAGE_POOR;
+ mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
+ if (mark) {
+ Nmark ++;
+ image[i].code |= ID_IMAGE_POOR;
+ } else {
+ image[i].code &= ~ID_IMAGE_POOR;
+ }
+ }
+
+ fprintf (stderr, "%d images marked poor\n", Nmark);
+ initstats (STATMODE);
+ free (mlist);
+ free (slist);
+ free (dlist);
+}
+
+void plot_images () {
+
+ int i, j, m, c, N, Nmax, bin;
+ short Mrel, Mcal, Mmos;
+ double *xlist, *Mlist, *dlist;
+ Graphdata graphdata;
+
+ if (FREEZE_IMAGES) return;
+
+ ALLOCATE (xlist, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+ ALLOCATE (Mlist, double, Nimage);
+
+ /**** dMcal vs airmass ****/
+ for (i = 0; i < Nimage; i++) {
+ Mlist[i] = 0.001*image[i].Mcal;
+ dlist[i] = 0.001*image[i].dMcal;
+ xlist[i] = 0.001*image[i].secz;
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png");
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);
+
+# define NBIN 200
+ REALLOCATE (xlist, double, NBIN);
+ REALLOCATE (Mlist, double, NBIN);
+
+ /**** dMcal histgram ****/
+ for (i = 0; i < NBIN; i++) xlist[i] = 0.0025*i;
+ bzero (Mlist, NBIN*sizeof(double));
+ for (i = 0; i < Nimage; i++) {
+ bin = 0.001*image[i].dMcal / 0.0025;
+ bin = MAX (0, MIN (NBIN - 1, bin));
+ Mlist[bin] += 1.0;
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.style = 1;
+ plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png");
+
+ free (dlist);
+ free (xlist);
+ free (Mlist);
+}
+
+StatType statsImageN (Catalog *catalog) {
+
+ int i, j, k, m, c, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (FREEZE_IMAGES) return (stats);
+
+ ALLOCATE (list, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+
+ n = 0;
+ for (i = 0; i < Nimage; i++) {
+ if (image[i].code & IMAGE_BAD) continue;
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (getMcal (m, c) == NO_MAG) continue;
+ if (getMmos (m, c) == NO_MAG) continue;
+ if (getMgrid (m, c) == NO_MAG) continue;
+ N++;
+ }
+ list[n] = N;
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsImageX (Catalog *catalog) {
+
+ int i, j, k, m, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (FREEZE_IMAGES) return (stats);
+
+ ALLOCATE (list, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+
+ n = 0;
+ for (i = 0; i < Nimage; i++) {
+
+ if (image[i].code & IMAGE_BAD) continue;
+
+ list[n] = pow (10.0, 0.01*image[i].Xm);
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsImageM (Catalog *catalog) {
+
+ int i, j, k, m, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (FREEZE_IMAGES) return (stats);
+
+ ALLOCATE (list, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+
+ n = 0;
+ for (i = 0; i < Nimage; i++) {
+
+ if (image[i].code & IMAGE_BAD) continue;
+
+ list[n] = 0.001*image[i].Mcal;
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsImagedM (Catalog *catalog) {
+
+ int i, j, k, m, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (FREEZE_IMAGES) return (stats);
+
+ ALLOCATE (list, double, Nimage);
+ ALLOCATE (dlist, double, Nimage);
+
+ n = 0;
+ for (i = 0; i < Nimage; i++) {
+
+ if (image[i].code & IMAGE_BAD) continue;
+
+ list[n] = 0.001*image[i].dMcal;
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+Image *getimages (int *N) {
+
+ *N = Nimage;
+ return (image);
+}
+
+Image *getimage (int N) {
+ return (&image[N]);
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/MosaicOps.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/MosaicOps.c (revision 2485)
@@ -0,0 +1,671 @@
+# include "relphot.h"
+
+Image *getimages (int *N);
+
+static int Nmosaic;
+static Mosaic *mosaic;
+
+static int *Nimlist;
+static int **imlist; /* mosaic -> image[] */
+static int **bin; /* catalog, measure -> mosaic */
+
+static int **clist; /* mosaic -> catalog[] */
+static int **mlist; /* mosiac -> measure[] */
+static int *Nlist;
+static int *NLIST;
+
+/* find mosaic frames (unique time periods & photcode name matches mosaic) */
+void initMosaics (Image *image, int Nimage) {
+
+ int i, j, status, found, NMOSAIC, *NIMLIST;
+ unsigned int start, stop;
+ char *pname;
+
+ if (!MOSAICNAME[0]) return;
+
+ Nmosaic = 0;
+ NMOSAIC = 10;
+ ALLOCATE (mosaic, Mosaic, NMOSAIC);
+
+ ALLOCATE (imlist, int *, NMOSAIC);
+ ALLOCATE (Nimlist, int, NMOSAIC);
+ ALLOCATE (NIMLIST, int, NMOSAIC);
+
+ /* generate list of unique mosaics */
+ for (i = 0; i < Nimage; i++) {
+
+ /* select valid mosaic images by photcode */
+ pname = GetPhotcodeNamebyCode (image[i].source);
+ status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
+ if (status) continue;
+
+ /* set image time range */
+ start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
+ stop = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
+
+ /* find existing mosaic with this time range */
+ found = FALSE;
+ for (j = 0; !found && (j < Nmosaic); j++) {
+ if (stop < mosaic[j].start) continue;
+ if (start > mosaic[j].stop) continue;
+ found = TRUE;
+
+ /* add image to mosaic image list */
+ imlist[j][Nimlist[j]] = i;
+ Nimlist[j] ++;
+ if (Nimlist[j] == NIMLIST[j]) {
+ NIMLIST[j] += 10;
+ REALLOCATE (imlist[j], int, NIMLIST[j]);
+ }
+
+ }
+ if (found) continue;
+
+ /* a new mosaic, define ranges */
+ mosaic[Nmosaic].start = start;
+ mosaic[Nmosaic].stop = stop;
+ mosaic[Nmosaic].Mcal = 0.0;
+ mosaic[Nmosaic].dMcal = 0.0;
+ mosaic[Nmosaic].Xm = 0.0;
+ mosaic[Nmosaic].code = image[i].code;
+ mosaic[Nmosaic].secz = image[i].secz;
+
+ /* add image to mosaic image list */
+ NIMLIST[Nmosaic] = 10;
+ Nimlist[Nmosaic] = 1;
+ ALLOCATE (imlist[Nmosaic], int, NIMLIST[Nmosaic]);
+ imlist[Nmosaic][0] = i;
+
+ Nmosaic ++;
+ if (Nmosaic == NMOSAIC) {
+ NMOSAIC += 10;
+ REALLOCATE (mosaic, Mosaic, NMOSAIC);
+ REALLOCATE (imlist, int *, NMOSAIC);
+ REALLOCATE (Nimlist, int, NMOSAIC);
+ REALLOCATE (NIMLIST, int, NMOSAIC);
+ }
+ }
+
+ initMosaicGrid (image, Nimage);
+ return;
+}
+
+void initMosaicGrid (Image *image, int Nimage) {
+
+ /* find max dR, dD range for all mosaics */
+ /* define mosaic.coords to cover dR, dD */
+ /* send results to initGridBins */
+
+ int i, j, m, NX, NY;
+ int dXmax, dYmax;
+ double dR, dD, dS, dX, dY;
+ double R, D, Rmin, Rmax, Dmin, Dmax;
+ double Mcal, dMcal, Xm;
+
+ dXmax = dYmax = 0.0;
+ for (i = 0; i < Nmosaic; i++) {
+ Dmin = Rmin = 360.0;
+ Dmax = Rmax = -360.0;
+ dS = 0.0;
+ Mcal = dMcal = Xm = 0;
+ for (j = 0; j < Nimlist[i]; j++) {
+ m = imlist[i][j];
+ NX = image[m].NX;
+ NY = image[m].NY;
+ dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
+ XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
+ Rmin = MIN (Rmin, R);
+ Rmax = MAX (Rmax, R);
+ Dmin = MIN (Dmin, D);
+ Dmax = MAX (Dmax, D);
+ XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
+ Rmin = MIN (Rmin, R);
+ Rmax = MAX (Rmax, R);
+ Dmin = MIN (Dmin, D);
+ Dmax = MAX (Dmax, D);
+ XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
+ Rmin = MIN (Rmin, R);
+ Rmax = MAX (Rmax, R);
+ Dmin = MIN (Dmin, D);
+ Dmax = MAX (Dmax, D);
+ XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
+ Rmin = MIN (Rmin, R);
+ Rmax = MAX (Rmax, R);
+ Dmin = MIN (Dmin, D);
+ Dmax = MAX (Dmax, D);
+ Mcal += image[m].Mcal;
+ dMcal += image[m].dMcal;
+ Xm += image[m].Xm;
+ /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
+ image[m].Mcal = 0.0;
+ image[m].dMcal = image[m].Xm = NO_MAG;
+ }
+ dS /= Nimlist[i];
+ strcpy (mosaic[i].coords.ctype, "RA---TAN");
+ mosaic[i].coords.crval1 = Rmin;
+ mosaic[i].coords.crval2 = Dmin;
+ mosaic[i].coords.crpix1 = 0.0;
+ mosaic[i].coords.crpix2 = 0.0;
+ mosaic[i].coords.cdelt1 = dS;
+ mosaic[i].coords.cdelt2 = dS;
+ mosaic[i].coords.pc1_1 = 1.0;
+ mosaic[i].coords.pc2_2 = 1.0;
+ mosaic[i].coords.pc1_2 = 0.0;
+ mosaic[i].coords.pc2_1 = 0.0;
+ RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
+
+ mosaic[i].Mcal = Mcal / Nimlist[i];
+ mosaic[i].dMcal = dMcal / Nimlist[i];
+ mosaic[i].Xm = Xm / Nimlist[i];
+ }
+ if (!USE_GRID) return;
+
+ dXmax = MAX (dXmax, dX);
+ dYmax = MAX (dYmax, dY);
+ initGrid (dXmax, dYmax);
+ return;
+}
+
+void setMcalFinal (Image *image, int Nimage) {
+
+ int i, j, im;
+
+ if (!MOSAICNAME[0]) return;
+
+ for (i = 0; i < Nmosaic; i++) {
+ for (j = 0; j < Nimlist[i]; j++) {
+ im = imlist[i][j];
+ image[im].Mcal = mosaic[i].Mcal;
+ image[im].dMcal = mosaic[i].dMcal;
+ image[im].Xm = mosaic[i].Xm;
+ image[im].code |= (mosaic[i].code & ID_IMAGE_FEW);
+ image[im].code |= (mosaic[i].code & ID_IMAGE_POOR);
+ }
+ }
+}
+
+void initMosaicBins (Catalog *catalog, int Ncatalog) {
+
+ int i, j;
+
+ /* measure -> mosaic */
+ if (!MOSAICNAME[0]) return;
+ ALLOCATE (bin, int *, Ncatalog);
+ for (i = 0; i < Ncatalog; i++) {
+ ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+ for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
+ }
+
+ /* mosaic -> measure */
+ ALLOCATE (Nlist, int, Nmosaic);
+ ALLOCATE (NLIST, int, Nmosaic);
+ ALLOCATE (clist, int *, Nmosaic);
+ ALLOCATE (mlist, int *, Nmosaic);
+
+ for (i = 0; i < Nmosaic; i++) {
+ Nlist[i] = 0;
+ NLIST[i] = 100;
+ ALLOCATE (clist[i], int, NLIST[i]);
+ ALLOCATE (mlist[i], int, NLIST[i]);
+ }
+}
+
+void freeMosaicBins (Ncatalog) {
+
+ int i;
+
+ /* measure -> mosaic */
+ if (!MOSAICNAME[0]) return;
+ for (i = 0; i < Ncatalog; i++) {
+ free (bin[i]);
+ }
+ free (bin);
+
+ /* mosaic -> measure */
+ for (i = 0; i < Nmosaic; i++) {
+ free (clist[i]);
+ free (mlist[i]);
+ }
+ free (Nlist);
+ free (NLIST);
+ free (clist);
+ free (mlist);
+}
+
+int findMosaics (Catalog *catalog, int Ncatalog) {
+
+ int i, j, ecode;
+
+ if (!MOSAICNAME[0]) return (FALSE);
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Nmeasure; j++) {
+ if (TimeSelect) {
+ if (catalog[i].measure[j].t < TSTART) continue;
+ if (catalog[i].measure[j].t > TSTOP) continue;
+ }
+ ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
+ if (photcode[0].code != ecode) continue;
+ matchMosaics (catalog, j, i);
+ }
+ }
+
+}
+
+void matchMosaics (Catalog *catalog, int meas, int cat) {
+
+ int i, ave;
+ double ra, dec, X, Y;
+
+ for (i = 0; i < Nmosaic; i++) {
+ if (catalog[cat].measure[meas].t < mosaic[i].start) continue;
+ if (catalog[cat].measure[meas].t > mosaic[i].stop) continue;
+
+# ifdef GRID_V1
+ if (USE_GRID) {
+ ave = catalog[cat].measure[meas].averef;
+ ra = catalog[cat].average[ave].R - catalog[cat].measure[meas].dR / 360000.0;
+ dec = catalog[cat].average[ave].D - catalog[cat].measure[meas].dD / 360000.0;
+
+ /* X,Y always positive-definite in range 0,0 - dX, dY */
+ RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords);
+ setGridMeasure (meas, cat, X, Y);
+ }
+# endif
+
+ bin[cat][meas] = i;
+
+ clist[i][Nlist[i]] = cat;
+ mlist[i][Nlist[i]] = meas;
+ Nlist[i] ++;
+
+ if (Nlist[i] == NLIST[i]) {
+ NLIST[i] += 100;
+ REALLOCATE (clist[i], int, NLIST[i]);
+ REALLOCATE (mlist[i], int, NLIST[i]);
+ }
+ return;
+ }
+ fprintf (stderr, "missed measurement\n");
+ return;
+}
+
+short getMmos (int meas, int cat) {
+
+ int i;
+ short value;
+
+ if (!MOSAICNAME[0]) return (0);
+ i = bin[cat][meas];
+ if (i == -1) return (NO_MAG);
+
+ if (mosaic[i].code & IMAGE_BAD) return (NO_MAG);
+ value = mosaic[i].Mcal;
+ return (value);
+}
+
+int setMmos (Catalog *catalog, int PoorImages) {
+
+ int i, j, m, c, n, N, Nmax, mark, bad;
+ short Msys, Mrel, Mcal, Mgrid;
+ double *list, *dlist, *Mlist, *dMlist;
+ StatType stats;
+ Image *image;
+ Graphdata graphdata;
+
+ if (!MOSAICNAME[0]) return (FALSE);
+
+ image = getimages (&N);
+
+ if (PoorImages) {
+ IMAGE_BAD = 0;
+ }
+
+ Nmax = 0;
+ for (i = 0; i < Nmosaic; i++) {
+ Nmax = MAX (Nmax, Nlist[i]);
+ }
+ ALLOCATE (list, double, Nmax);
+ ALLOCATE (dlist, double, Nmax);
+ ALLOCATE (Mlist, double, Nmax);
+ ALLOCATE (dMlist, double, Nmax);
+
+ for (i = 0; i < Nmosaic; i++) {
+
+ /* on PoorImages run, skip good images */
+ if (PoorImages) {
+ bad = mosaic[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
+ if (!bad) continue;
+ }
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (catalog[c].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, c)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
+ if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue;
+
+ n = catalog[c].measure[m].averef;
+ Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+ list[N] = Msys - Mrel - Mcal - Mgrid;
+ dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+ Mlist[N] = 0.001*Msys;
+ dMlist[N] = 0.001*list[N];
+ N++;
+ }
+ /* Nlist[i] is all measurements, N is good measurements */
+
+ /* too few good measurements or too many bad measurements (skip in PoorImages run) */
+
+ if (!PoorImages) {
+ mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
+ if (mark) {
+ fprintf (stderr, "marked image %s (%d), (%d < %d) || (%d < %f*%d)\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]);
+ mosaic[i].code |= ID_IMAGE_FEW;
+ } else {
+ mosaic[i].code &= ~ID_IMAGE_FEW;
+ }
+ }
+ liststats (list, dlist, N, &stats);
+ if (PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
+ mosaic[i].Mcal = stats.mean;
+ mosaic[i].dMcal = stats.sigma;
+ mosaic[i].Xm = 100.0*log10(stats.chisq);
+ }
+ free (list);
+ free (dlist);
+ free (Mlist);
+ free (dMlist);
+
+ if (PoorImages) {
+ IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
+ STAR_BAD = ID_STAR_POOR | ID_STAR_FEW;
+ MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
+ }
+ return (TRUE);
+}
+
+StatType statsMosaicM (Catalog *catalog) {
+
+ int i, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (!MOSAICNAME[0]) return (stats);
+
+ ALLOCATE (list, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+
+ n = 0;
+ for (i = 0; i < Nmosaic; i++) {
+ if (mosaic[i].code & IMAGE_BAD) continue;
+ list[n] = 0.001*mosaic[i].Mcal;
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsMosaicdM (Catalog *catalog) {
+
+ int i, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (!MOSAICNAME[0]) return (stats);
+
+ ALLOCATE (list, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+
+ n = 0;
+ for (i = 0; i < Nmosaic; i++) {
+
+ if (mosaic[i].code & IMAGE_BAD) continue;
+ list[n] = 0.001*mosaic[i].dMcal;
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsMosaicN (Catalog *catalog) {
+
+ int i, j, k, m, c, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (!MOSAICNAME[0]) return (stats);
+
+ ALLOCATE (list, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+
+ n = 0;
+ for (i = 0; i < Nmosaic; i++) {
+ if (mosaic[i].code & IMAGE_BAD) continue;
+
+ N = 0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ if (getMcal (m, c) == NO_MAG) continue;
+ if (getMgrid (m, c) == NO_MAG) continue;
+ if (getMrel (catalog, m, c) == NO_MAG) continue;
+ N++;
+ }
+ list[n] = N;
+ dlist[n] = 1;
+ n++;
+ }
+ fprintf (stderr, "Nmosaic: %d, n: %d\n", Nmosaic, n);
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsMosaicX (Catalog *catalog) {
+
+ int i, n, N, Nmax;
+ double *list, *dlist;
+ StatType stats;
+
+ bzero (&stats, sizeof (StatType));
+ if (!MOSAICNAME[0]) return (stats);
+
+ ALLOCATE (list, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+
+ n = 0;
+ for (i = 0; i < Nmosaic; i++) {
+
+ if (mosaic[i].code & IMAGE_BAD) continue;
+ list[n] = pow(10.0, 0.01*mosaic[i].Xm);
+ dlist[n] = 1;
+ n++;
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+/* mark mosaic if: abs(Mcal - ) too large, dMcal too large */
+void clean_mosaics () {
+
+ int i, j, N, mark, Nmark;
+ double *mlist, *slist, *dlist;
+ double MaxOffset, MedOffset, MaxScatter;
+ StatType stats;
+ Image *image;
+
+ if (!MOSAICNAME[0]) return;
+
+ if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
+
+ ALLOCATE (mlist, double, Nmosaic);
+ ALLOCATE (slist, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+
+ for (i = N = 0; i < Nmosaic; i++) {
+ if (mosaic[i].code & IMAGE_BAD) continue;
+ mlist[N] = mosaic[i].Mcal;
+ slist[N] = mosaic[i].dMcal;
+ dlist[N] = 1;
+ N++;
+ }
+ initstats ("MEAN");
+ liststats (mlist, dlist, N, &stats);
+ MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
+ MedOffset = stats.median;
+ liststats (slist, dlist, N, &stats);
+ MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
+ fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", 0.001*MedOffset, 0.001*stats.median, 0.001*MaxScatter, 0.001*MaxOffset);
+
+ Nmark = 0;
+ for (i = 0; i < Nmosaic; i++) {
+ mark = FALSE;
+ mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
+ if (mark) {
+ Nmark ++;
+ mosaic[i].code |= ID_IMAGE_POOR;
+ } else {
+ mosaic[i].code &= ~ID_IMAGE_POOR;
+ }
+ }
+
+ fprintf (stderr, "%d mosaics marked poor\n", Nmark);
+ initstats (STATMODE);
+ free (mlist);
+ free (slist);
+ free (dlist);
+}
+
+void plot_mosaic_fields (Catalog *catalog) {
+
+ int i, j, m, c, N, ave, Nimage;
+ short Mrel, Mcal, Mmos;
+ double *xlist, *ylist;
+ double Xmin, Xmax, Ymin, Ymax;
+ char string[64];
+ Image *image;
+ Graphdata graphdata;
+
+ if (!MOSAICNAME[0]) return;
+
+ image = getimages (&Nimage);
+
+ N = 0;
+ for (i = 0; i < Nmosaic; i++)
+ N = MAX (N, Nlist[i]);
+
+ ALLOCATE (xlist, double, N);
+ ALLOCATE (ylist, double, N);
+
+ for (i = 0; i < Nmosaic; i++) {
+ N = 0;
+ Xmin = Ymin = +360.0;
+ Xmax = Ymax = -360.0;
+ for (j = 0; j < Nlist[i]; j++) {
+
+ m = mlist[i][j];
+ c = clist[i][j];
+
+ ave = catalog[c].measure[m].averef;
+ xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 360000.0;
+ ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 360000.0;
+ N++;
+ }
+
+ sprintf (string, "Mosaic %d", i);
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, xlist, ylist, N, string, NULL);
+ }
+
+ free (ylist);
+ free (xlist);
+}
+
+void plot_mosaics () {
+
+ int i, j, m, c, N, Nmax, bin;
+ short Mrel, Mcal, Mmos;
+ double *xlist, *Mlist, *dlist;
+ Graphdata graphdata;
+
+ if (!MOSAICNAME[0]) return;
+
+ ALLOCATE (xlist, double, Nmosaic);
+ ALLOCATE (dlist, double, Nmosaic);
+ ALLOCATE (Mlist, double, Nmosaic);
+
+ for (i = 0; i < Nmosaic; i++) {
+ Mlist[i] = 0.001*mosaic[i].Mcal;
+ dlist[i] = 0.001*mosaic[i].dMcal;
+ xlist[i] = 0.001*mosaic[i].secz;
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.xmin = 0.95;
+ graphdata.xmax = 2.50;
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, xlist, Mlist, Nmosaic, "airmass vs Mcal", "airmass.png");
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, Mlist, dlist, Nmosaic, "Mcal vs dMcal", NULL);
+
+# define NBIN 200
+ REALLOCATE (xlist, double, NBIN);
+ REALLOCATE (Mlist, double, NBIN);
+
+ /**** dMcal histgram ****/
+ for (i = 0; i < NBIN; i++) xlist[i] = 0.0025*i;
+ bzero (Mlist, NBIN*sizeof(double));
+ for (i = 0; i < Nmosaic; i++) {
+ bin = 0.001*mosaic[i].dMcal / 0.0025;
+ bin = MAX (0, MIN (NBIN - 1, bin));
+ Mlist[bin] += 1.0;
+ }
+ plot_defaults (&graphdata);
+ graphdata.style = 1;
+ plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png");
+
+ free (dlist);
+ free (xlist);
+ free (Mlist);
+}
+
+int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage) {
+
+ int i, Imax, Nmax;
+
+ Imax = 0;
+ Nmax = Nimlist[0];
+ for (i = 0; i < Nmosaic; i++) {
+ if (Nimlist[i] > Nmax) {
+ Imax = i;
+ Nmax = Nimlist[i];
+ }
+ }
+
+ *refmosaic = &mosaic[Imax];
+ *Nimage = Nmax;
+ return (imlist[Imax]);
+}
Index: /branches/elixir/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/StarOps.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/StarOps.c (revision 2485)
@@ -0,0 +1,528 @@
+# include "relphot.h"
+
+static int Nmax;
+static double *list;
+static double *dlist;
+
+void initMrel (Catalog *catalog, int Ncatalog) {
+
+ int i, j;
+
+ Nmax = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ Nmax = MAX (Nmax, catalog[i].average[j].Nm);
+ }
+ }
+
+ ALLOCATE (list, double, MAX (1, Nmax));
+ ALLOCATE (dlist, double, MAX (1, Nmax));
+}
+
+short getMrel (Catalog *catalog, int meas, int cat) {
+
+ int ave;
+ short value;
+
+ ave = catalog[cat].measure[meas].averef;
+ if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG);
+
+ value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
+ return (value);
+}
+
+int setMrel (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m, N;
+ short Msys, Mcal, Mmos, Mgrid;
+ StatType stats;
+
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* calculate the average value for a single star */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ m = catalog[i].average[j].offset;
+
+ N = 0;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+
+ Msys = iPhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+ list[N] = Msys - Mcal - Mmos - Mgrid;
+ dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+ if (N < STAR_TOOFEW) { /* too few measurements */
+ catalog[i].average[j].code |= ID_STAR_FEW;
+ } else {
+ catalog[i].average[j].code &= ~ID_STAR_FEW;
+ }
+
+ liststats (list, dlist, N, &stats);
+
+ if (PhotPrimary) {
+ catalog[i].average[j].M = stats.mean;
+ catalog[i].average[j].dM = stats.sigma;
+ catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
+ } else {
+ catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean;
+ catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
+ catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
+ }
+ }
+ }
+}
+
+int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) {
+
+ int i, j, k, m, N, Nmax;
+ short Msys, Mcal, Mmos, Mgrid;
+ double *list, *dlist;
+ StatType stats;
+
+ /* Nmeasure is now different, need to reallocate */
+ Nmax = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ Nmax = MAX (Nmax, catalog[i].average[j].Nm);
+ }
+ }
+ ALLOCATE (list, double, MAX (1, Nmax));
+ ALLOCATE (dlist, double, MAX (1, Nmax));
+
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* skip stars already calibrated */
+ if (catalog[i].found[j]) continue;
+
+ N = 0;
+ m = catalog[i].average[j].offset;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+
+ Msys = iPhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+ list[N] = Msys - Mcal - Mmos - Mgrid;
+ dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+ if (N < 1) continue;
+
+ liststats (list, dlist, N, &stats);
+ if (mark) catalog[i].found[j] = TRUE;
+
+ /* use sigma or error in dM for output? */
+ if (PhotPrimary) {
+ catalog[i].average[j].M = stats.mean;
+ catalog[i].average[j].dM = MAX (stats.sigma, stats.error);
+ catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
+ } else {
+ catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean;
+ catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
+ catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
+ }
+ }
+ }
+
+ free (list);
+ free (dlist);
+}
+
+/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_NOCAL */
+int setMcalOutput (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m;
+ short Mcal, Mmos, Mgrid;
+
+ MEAS_BAD = ID_MEAS_NOCAL;
+ IMAGE_BAD = ID_IMAGE_NOCAL;
+
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ m = catalog[i].average[j].offset;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+ catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
+ }
+ }
+ }
+ return (TRUE);
+}
+
+void clean_stars (Catalog *catalog, int Ncatalog) {
+
+ int i, j, Ndel, Nave, Ntot, mark;
+ short dM, Xm;
+ double Chisq, MaxScatter, MaxChisq;
+ double *xlist, *slist, *dlist;
+ StatType stats;
+
+ if (VERBOSE) fprintf (stderr, "marking poor stars\n");
+
+ /* find Xm median -> ChiSq lim must be > median */
+ for (i = Ntot = 0; i < Ncatalog; i++) {
+ Ntot += catalog[i].Naverage;
+ }
+ ALLOCATE (xlist, double, Ntot);
+ ALLOCATE (slist, double, Ntot);
+ ALLOCATE (dlist, double, Ntot);
+ for (i = Ntot = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+ if (Xm == -1) continue;
+ Chisq = pow (10.0, 0.01*Xm);
+ xlist[Ntot] = Chisq;
+ slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+ dlist[Ntot] = 1;
+ Ntot ++;
+ }
+ }
+
+ initstats ("MEAN");
+ liststats (xlist, dlist, Ntot, &stats);
+ MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
+ liststats (slist, dlist, Ntot, &stats);
+ MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
+ fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", 0.001*MaxScatter, MaxChisq);
+
+ Ndel = Nave = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+ Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+ Chisq = pow (10.0, 0.01*Xm);
+ mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq);
+ if (mark) {
+ catalog[i].average[j].code |= ID_STAR_POOR;
+ Ndel ++;
+ } else {
+ catalog[i].average[j].code &= ~ID_STAR_POOR;
+ }
+ Nave ++;
+ }
+ }
+ fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
+ initstats (STATMODE);
+ free (xlist);
+ free (slist);
+ free (dlist);
+}
+
+void clean_measures (Catalog *catalog, int Ncatalog, int final) {
+
+ int i, j, k, m, N, Ndel, Nave, Nmax, image_bad, TOOFEW;
+ int *ilist;
+ double *tlist, *list, *dlist, Ns;
+ short Msys, Mcal, Mmos, Mgrid;
+ StatType stats;
+
+ if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+ /* Nmeasure is now different, need to reallocate */
+ Nmax = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ Nmax = MAX (Nmax, catalog[i].average[j].Nm);
+ }
+ }
+ ALLOCATE (list, double, MAX (1, Nmax));
+ ALLOCATE (dlist, double, MAX (1, Nmax));
+ ALLOCATE (ilist, int, Nmax);
+ ALLOCATE (tlist, double, Nmax);
+
+ /* it makes no sense to mark 3-sigma outliers with <5 measurements */
+ TOOFEW = MAX (5, STAR_TOOFEW);
+
+ Ns = 3;
+ Ndel = Nave = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* skip bad stars to prevent them from becoming good (on inner sample) */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+
+ /* on final processing, skip stars already measured */
+ if (final && catalog[i].found[j]) continue;
+
+ /* accumulate list of valid measurements */
+ m = catalog[i].average[j].offset;
+ N = 0;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ /* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+
+ Msys = iPhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+ list[N] = Msys - Mcal - Mmos - Mgrid;
+ dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+ N++;
+ }
+ if (N < TOOFEW) continue;
+
+ /* 3-sigma clip based on stats of inner 50% */
+ initstats ("INNER_MEAN");
+ liststats (list, dlist, N, &stats);
+ stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
+ for (k = m = 0; k < N; k++) {
+ if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
+ list[m] = list[k];
+ m++;
+ }
+ }
+ initstats ("MEAN");
+ liststats (list, dlist, m, &stats);
+ stats.sigma = MAX (MIN_ERROR, stats.sigma);
+
+ /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
+ image_bad = IMAGE_BAD;
+ IMAGE_BAD = ID_IMAGE_NOCAL;
+ m = catalog[i].average[j].offset;
+ N = 0;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ /* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+
+ Msys = iPhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+ list[N] = Msys - Mcal - Mmos - Mgrid;
+ dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+ ilist[N] = m;
+ N++;
+ Nave ++;
+ }
+ if (N < TOOFEW) continue;
+
+ /* mark bad measures */
+ for (k = 0; k < N; k++) {
+ if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
+ catalog[i].measure[ilist[k]].flags |= ID_MEAS_POOR;
+ Ndel ++;
+ }
+ }
+ IMAGE_BAD = image_bad;
+ }
+ }
+ initstats (STATMODE);
+ fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
+ free (ilist);
+ free (tlist);
+}
+
+StatType statsStarN (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m, n, N, Ntot;
+ double *list, *dlist;
+ StatType stats;
+
+ Ntot = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ Ntot += catalog[i].Naverage;
+ }
+
+ ALLOCATE (list, double, Ntot);
+ ALLOCATE (dlist, double, Ntot);
+
+ n = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* calculate the average value for a single star */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ m = catalog[i].average[j].offset;
+
+ N = 0;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ if (getMcal (m, i) == NO_MAG) continue;
+ if (getMmos (m, i) == NO_MAG) continue;
+ if (getMgrid (m, i) == NO_MAG) continue;
+ N++;
+ }
+
+ list[n] = N;
+ dlist[n] = 1;
+ n++;
+ }
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsStarX (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m, n, N, Ntot, Xm;
+ double *list, *dlist;
+ StatType stats;
+
+ Ntot = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ Ntot += catalog[i].Naverage;
+ }
+
+ ALLOCATE (list, double, Ntot);
+ ALLOCATE (dlist, double, Ntot);
+
+ n = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* calculate the average value for a single star */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+
+ Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+ if (Xm == NO_MAG) continue;
+ list[n] = pow (10.0, 0.01*Xm);
+ dlist[n] = 1;
+ n++;
+ }
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+StatType statsStarS (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m, n, N, Ntot, dM;
+ double *list, *dlist;
+ StatType stats;
+
+ Ntot = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ Ntot += catalog[i].Naverage;
+ }
+
+ ALLOCATE (list, double, Ntot);
+ ALLOCATE (dlist, double, Ntot);
+
+ n = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* calculate the average value for a single star */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+
+ dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+ list[n] = 0.001*dM;
+ dlist[n] = 1;
+ n++;
+ }
+ }
+
+ liststats (list, dlist, n, &stats);
+ free (list);
+ free (dlist);
+ return (stats);
+}
+
+void plot_stars (Catalog *catalog, int Ncatalog) {
+
+ int i, j, m, c, N, bin;
+ short dMrel;
+ double *xlist, *Mlist;
+ Graphdata graphdata;
+
+# define NBIN 200
+ ALLOCATE (xlist, double, NBIN);
+ ALLOCATE (Mlist, double, NBIN);
+
+ for (i = 0; i < NBIN; i++) xlist[i] = 0.0025*i;
+ bzero (Mlist, NBIN*sizeof(double));
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+ bin = 0.001*dMrel / 0.0025;
+ bin = MAX (0, MIN (NBIN-1, bin));
+ Mlist[bin] += 1.0;
+ }
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.style = 1;
+ plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "dMhist.png");
+
+ free (xlist);
+ free (Mlist);
+}
+
+void plot_chisq (Catalog *catalog, int Ncatalog) {
+
+ int i, j, N, Ntotal, value;
+ double *xlist, *ylist;
+ Graphdata graphdata;
+
+ Ntotal = 0;
+ for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage;
+
+ ALLOCATE (xlist, double, Ntotal);
+ ALLOCATE (ylist, double, Ntotal);
+
+ N = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ xlist[N] = PhotPrimary ? 0.001*catalog[i].average[j].M : 0.001*catalog[i].secfilt[PhotNsec*j+PhotSec].M;
+ value = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+ if (value == NO_MAG) continue;
+ ylist[N] = 0.01*value;
+ N++;
+ }
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.ymin = -3.0;
+ plot_list (&graphdata, xlist, ylist, N, "chisq", NULL);
+ free (xlist);
+ free (ylist);
+}
+
+void plot_star_coords (Catalog *catalog, int Ncatalog) {
+
+ int i, j, m, c, N, bin;
+ short Mrel, Mcal, Mmos;
+ double *xlist, *ylist;
+ double Xmin, Ymin, Xmax, Ymax;
+ Graphdata graphdata;
+
+ N = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ N += catalog[i].Naverage;
+ }
+ ALLOCATE (xlist, double, N);
+ ALLOCATE (ylist, double, N);
+
+ N = 0;
+ Xmin = Ymin = +360.0;
+ Xmax = Ymax = -360.0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ xlist[N] = catalog[i].average[j].R;
+ ylist[N] = catalog[i].average[j].D;
+ N++;
+ }
+ }
+ plot_defaults (&graphdata);
+ plot_list (&graphdata, xlist, ylist, N, "coords", NULL);
+
+ free (xlist);
+ free (ylist);
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/args.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/args.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/args.c (revision 2485)
@@ -0,0 +1,175 @@
+# include "relphot.h"
+
+int args (int argc, char **argv) {
+
+ int N;
+
+ /* define time */
+ TimeSelect = FALSE;
+ if (N = get_argument (argc, argv, "-time")) {
+ remove_argument (N, &argc, argv);
+ if (!str_to_time (argv[N], &TSTART)) {
+ fprintf (stderr, "ERROR: syntax error\n");
+ return (FALSE);
+ }
+ remove_argument (N, &argc, argv);
+ if (!str_to_time (argv[N], &TSTOP)) {
+ fprintf (stderr, "ERROR: syntax error\n");
+ return (FALSE);
+ }
+ remove_argument (N, &argc, argv);
+ TimeSelect = TRUE;
+ }
+
+ VERBOSE = FALSE;
+ if ((N = get_argument (argc, argv, "-v"))) {
+ VERBOSE = TRUE;
+ remove_argument (N, &argc, argv);
+ }
+
+ PLOTSTUFF = FALSE;
+ if ((N = get_argument (argc, argv, "-plot"))) {
+ PLOTSTUFF = TRUE;
+ remove_argument (N, &argc, argv);
+ }
+
+ PLOTDELAY = 500000;
+ if ((N = get_argument (argc, argv, "-plotdelay"))) {
+ remove_argument (N, &argc, argv);
+ PLOTDELAY = 1e6*atof(argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ strcpy (STATMODE, "CHI_INNER_WTMEAN");
+ if ((N = get_argument (argc, argv, "-statmode"))) {
+ remove_argument (N, &argc, argv);
+ strcpy (STATMODE, argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ NLOOP = 8;
+ if ((N = get_argument (argc, argv, "-n"))) {
+ remove_argument (N, &argc, argv);
+ NLOOP = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ RESET = FALSE;
+ if ((N = get_argument (argc, argv, "-reset"))) {
+ remove_argument (N, &argc, argv);
+ RESET = TRUE;
+ }
+
+ UPDATE = FALSE;
+ if ((N = get_argument (argc, argv, "-update"))) {
+ remove_argument (N, &argc, argv);
+ UPDATE = TRUE;
+ }
+
+ SHOW_PARAMS = FALSE;
+ if ((N = get_argument (argc, argv, "-params"))) {
+ remove_argument (N, &argc, argv);
+ SHOW_PARAMS = TRUE;
+ }
+
+ PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0;
+ if ((N = get_argument (argc, argv, "-plrange"))) {
+ remove_argument (N, &argc, argv);
+ PlotMmin = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ PlotMmax = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ PlotdMmin = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ PlotdMmax = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ /* group images by mosaic, find Mmos */
+ MOSAICNAME[0] = 0;
+ if ((N = get_argument (argc, argv, "-mosaic"))) {
+ remove_argument (N, &argc, argv);
+ strcpy (MOSAICNAME, argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ FREEZE_IMAGES = FALSE;
+ if ((N = get_argument (argc, argv, "-imfreeze"))) {
+ remove_argument (N, &argc, argv);
+ FREEZE_IMAGES = TRUE;
+ }
+
+ USE_GRID = FALSE;
+ if ((N = get_argument (argc, argv, "-grid"))) {
+ remove_argument (N, &argc, argv);
+ USE_GRID = TRUE;
+ if (!MOSAICNAME[0]) {
+ fprintf (stderr, "-grid is only valid with -mosaic\n");
+ exit (2);
+ }
+ }
+
+ MIN_ERROR = 1;
+ if ((N = get_argument (argc, argv, "-minerror"))) {
+ remove_argument (N, &argc, argv);
+ MIN_ERROR = 1000.0*atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ /* require MIN_ERROR > 0 */
+ }
+
+ AreaSelect = FALSE;
+ if ((N = get_argument (argc, argv, "-area"))) {
+ remove_argument (N, &argc, argv);
+ AreaXmin = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ AreaXmax = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ AreaYmin = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ AreaYmax = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ AreaSelect = TRUE;
+ }
+
+ ImagSelect = FALSE;
+ if ((N = get_argument (argc, argv, "-instmag"))) {
+ remove_argument (N, &argc, argv);
+ ImagMin = 1000*atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ ImagMax = 1000*atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ ImagSelect = TRUE;
+ }
+
+ if (argc != 3) {
+ fprintf (stderr, "ERROR: USAGE: nrphot (region) (photcode)\n");
+ fprintf (stderr, " options: \n");
+ fprintf (stderr, " -time (start) (stop)\n");
+ fprintf (stderr, " -v\n");
+ fprintf (stderr, " -plot\n");
+ fprintf (stderr, " -plotdelay (seconds)\n");
+ fprintf (stderr, " -statmode (mode)\n");
+ fprintf (stderr, " -n (nloop)\n");
+ fprintf (stderr, " -reset\n");
+ fprintf (stderr, " -update\n");
+ fprintf (stderr, " -params\n");
+ fprintf (stderr, " -mosaic (mosaic)\n");
+ fprintf (stderr, " -imfreeze\n");
+ fprintf (stderr, " -grid\n");
+ fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n");
+ fprintf (stderr, " -instmag min max\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ fprintf (stderr, " -\n");
+ exit (2);
+ }
+
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/bcatalog.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/bcatalog.c (revision 2485)
@@ -0,0 +1,119 @@
+# include "relphot.h"
+
+int bcatalog (Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion) {
+
+ int i, j, offset, ecode;
+ int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
+ int mag;
+
+ /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
+ NAVERAGE = 50;
+ NMEASURE = 1000;
+ ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
+ ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
+ ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+ Nmeasure = Naverage = 0;
+
+ /* exclude stars not in range or with too few measurements */
+ for (i = 0; i < catalog[0].Naverage; i++) {
+ if (catalog[0].average[i].Nm < 2) continue;
+ if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
+ if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
+ if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
+ if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
+
+ /* start with all stars good */
+ subcatalog[0].average[Naverage] = catalog[0].average[i];
+ subcatalog[0].average[Naverage].offset = Nmeasure;
+ for (j = 0; j < PhotNsec; j++) {
+ subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j];
+ }
+
+ if (RESET) {
+ short *p;
+ p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M;
+ *p = NO_MAG;
+ p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM;
+ *p = NO_MAG;
+ subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW;
+ subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR;
+ }
+
+ Nm = 0;
+ for (j = 0; j < catalog[0].average[i].Nm; j++) {
+
+ offset = catalog[0].average[i].offset + j;
+
+ /* select measurements by photcode */
+ ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
+ if (ecode != photcode[0].code) continue;
+
+ /* select measurements by time */
+ if (TimeSelect) {
+ if (catalog[0].measure[offset].t < TSTART) continue;
+ if (catalog[0].measure[offset].t > TSTOP) continue;
+ }
+
+ /* select measurements by quality */
+ if (catalog[0].measure[offset].dophot != 1) continue;
+
+ /* select measurements by mag limit */
+ mag = iPhotCat (&catalog[0].measure[offset]);
+ if (mag > MAG_LIM) continue;
+
+ /* select measurements by measurement error */
+ if (catalog[0].measure[offset].dM > SIGMA_LIM) continue;
+
+ /* select measurements by mag limit */
+ if (ImagSelect) {
+ mag = iPhotInst (&catalog[0].measure[offset]);
+ if (mag < ImagMin) continue;
+ if (mag > ImagMax) continue;
+ }
+
+ subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_SKIP;
+ subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
+ subcatalog[0].measure[Nmeasure].averef = Naverage;
+ if (RESET) {
+ subcatalog[0].measure[Nmeasure].Mcal = 0;
+ subcatalog[0].measure[Nmeasure].flags &= 0xff00;
+ subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_POOR;
+ subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_AREA;
+ subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_NOCAL;
+ }
+ Nmeasure ++;
+ Nm ++;
+ if (Nmeasure == NMEASURE) {
+ NMEASURE += 1000;
+ REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+ }
+ }
+ if (Nm < 2) { /* enough measurements in band? */
+ Nmeasure -= Nm;
+ continue;
+ }
+ subcatalog[0].average[Naverage].Nm = Nm;
+ Naverage ++;
+ if (Naverage == NAVERAGE) {
+ NAVERAGE += 50;
+ REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
+ REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
+ }
+ }
+ REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
+ REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
+ REALLOCATE (subcatalog[0].secfilt, SecFilt, PhotNsec*MAX (Naverage, 1));
+ subcatalog[0].Naverage = Naverage;
+ subcatalog[0].Nmeasure = Nmeasure;
+ subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+
+ if (catalog[0].Naverage > 0) free (catalog[0].average);
+ if (catalog[0].Nmeasure > 0) free (catalog[0].measure);
+ if (catalog[0].Naverage > 0) free (catalog[0].secfilt);
+
+ if (VERBOSE) {
+ fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
+ subcatalog[0].Naverage, subcatalog[0].Nmeasure);
+ }
+
+}
Index: /branches/elixir/Ohana/src/relphot/src/find_images.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/find_images.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/find_images.c (revision 2485)
@@ -0,0 +1,237 @@
+# include "relphot.h"
+
+/* this function returns a list of all images which overlap the given
+ set of region files. All images in the image catalog are tested
+ once, so there is no check that an image already has been included.
+ LineNum stores the locations in the Image database of the list of
+ images */
+
+Image *find_images (GSCRegion *region, int Nregion, int *Nimages, int **LineNum) {
+
+ Header header;
+ Image *timage, *image;
+ int i, j, k, m, found, nimage, Nimage, NIMAGE, NTIMAGE, Nloop, Nlast;
+ int n, Nim, status, InRange, ecode;
+ FILE *f;
+ double Xc[5], Yc[5], Xi[5], Yi[5], r, d, dx, dy;
+ char line[256];
+ int *line_number;
+ Coords tcoords;
+
+ if (VERBOSE) fprintf (stderr, "finding images\n");
+
+ /* we make positional comparisons in the projection of catalog */
+ tcoords.crval1 = 0.5*(region[0].RA[0] + region[0].RA[1]);
+ tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]);
+ tcoords.crpix1 = tcoords.crpix2 = 0.0;
+ tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+ tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+ tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+ strcpy (tcoords.ctype, "RA---TAN");
+
+ /* check if image datafile exists, get header, number of images */
+ if (!fits_read_header (ImageCat, &header)) {
+ fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat);
+ exit (0);
+ }
+ Nimage = 0;
+ fits_scan (&header, "NIMAGES", "%d", 1, &Nimage);
+ if (Nimage == 0) {
+ fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat);
+ exit (0);
+ }
+
+ /* get ready to read data on images */
+ f = fopen (ImageCat, "r");
+ if (f == (FILE *) NULL) {
+ fprintf (stderr, "ERROR: trouble opening Image catalog: %s (2)\n", ImageCat);
+ exit (0);
+ }
+ fseek (f, header.size, SEEK_SET);
+
+ /* load all images */
+ ALLOCATE (timage, Image, Nimage);
+ status = Fread (timage, sizeof(Image), Nimage, f, "image");
+ if (status != Nimage) {
+ fprintf (stderr, "ERROR: couldn't read images from image catalog: %s\n", ImageCat);
+ exit (0);
+ }
+
+ nimage = 0;
+ NIMAGE = 100;
+ ALLOCATE (image, Image, NIMAGE);
+ ALLOCATE (line_number, int, NIMAGE);
+
+ for (i = 0; i < Nimage; i++) {
+
+ /* select images by photcode */
+ ecode = GetPhotcodeEquivCodebyCode (timage[i].source);
+ if (ecode != photcode[0].code) continue;
+
+ /* select images by time */
+ if (TimeSelect) {
+ if (timage[i].tzero < TSTART) continue;
+ if (timage[i].tzero > TSTOP) continue;
+ }
+
+ /* ignore image set to SKIP
+ if (timage[i].code & ID_IMAGE_SKIP) continue; */
+
+ /* define image corners */
+ Xi[0] = 0; Yi[0] = 0;
+ Xi[1] = timage[i].NX; Yi[1] = 0;
+ Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
+ Xi[3] = 0; Yi[3] = timage[i].NY;
+ Xi[4] = 0; Yi[4] = 0;
+ found = FALSE;
+ /* transform to tcoords */
+ for (j = 0; j < 5; j++) {
+ XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords);
+ InRange = RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords);
+ if (!InRange) goto imskip;
+ }
+ /* compare with each region file */
+ for (m = 0; (m < Nregion) && !found; m++) {
+ /* define catalog corners */
+ Xc[0] = region[m].RA[0]; Yc[0] = region[m].DEC[0];
+ Xc[1] = region[m].RA[1]; Yc[1] = region[m].DEC[0];
+ Xc[2] = region[m].RA[1]; Yc[2] = region[m].DEC[1];
+ Xc[3] = region[m].RA[0]; Yc[3] = region[m].DEC[1];
+ Xc[4] = region[m].RA[0]; Yc[4] = region[m].DEC[0];
+ for (j = 0; j < 5; j++) {
+ r = Xc[j]; d = Yc[j];
+ RD_to_XY (&Xc[j], &Yc[j], r, d, &tcoords);
+ }
+ dx = 0.02*(Xc[2] - Xc[0]);
+ dy = 0.02*(Yc[2] - Yc[0]);
+ Xc[0] -= dx; Yc[0] -= dy;
+ Xc[1] += dx; Yc[1] -= dy;
+ Xc[2] += dx; Yc[2] += dy;
+ Xc[3] -= dx; Yc[3] += dy;
+ Xc[4] -= dx; Yc[4] -= dy;
+
+ /* check if image corner inside catalog */
+ for (j = 0; (j < 4) && !found; j++) {
+ found |= corner_check (&Xi[j], &Yi[j], &Xc[0], &Yc[0]);
+ }
+ /* check if catalog corner inside image */
+ for (j = 0; (j < 4) && !found; j++) {
+ found |= corner_check (&Xc[j], &Yc[j], &Xi[0], &Yi[0]);
+ }
+ /* check if edges cross */
+ for (j = 0; (j < 4) && !found; j++) {
+ for (k = 0; (k < 4) && !found; k++) {
+ found |= edge_check (&Xi[j], &Yi[j], &Xc[k], &Yc[k]);
+ }
+ }
+ if (found) {
+ image[nimage] = timage[i];
+ /* always allow a 'few' image to succeed, if possible */
+ if (image[nimage].code & ID_IMAGE_FEW) {
+ image[nimage].code &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR);
+ }
+ if (RESET) {
+ assignMcal (&image[nimage], (double *) NULL, -1);
+ image[nimage].dMcal = NO_MAG;
+ image[nimage].code &= ~ID_IMAGE_POOR;
+ }
+ line_number[nimage] = i;
+ nimage ++;
+ if (nimage == NIMAGE) {
+ NIMAGE += 100;
+ REALLOCATE (image, Image, NIMAGE);
+ REALLOCATE (line_number, int, NIMAGE);
+ }
+ }
+ }
+ imskip:
+ }
+
+ if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);
+
+ REALLOCATE (image, Image, MAX (nimage, 1));
+ REALLOCATE (line_number, int, MAX (nimage, 1));
+
+ free (timage);
+ *Nimages = nimage;
+ *LineNum = line_number;
+ fclose (f);
+ return (image);
+}
+
+/* check if line between points 0 and 1 of x1
+ crosses line between points 0 and 1 of x2 */
+int edge_check (double *x1, double *y1, double *x2, double *y2) {
+
+ double theta1, theta2;
+ double Theta1, Theta2;
+
+ theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]);
+ theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]);
+
+ if (theta1*theta2 < 0.0) {
+ return (FALSE);
+ }
+
+ if (fabs(theta1) < fabs(theta2)) {
+ return (FALSE);
+ }
+
+ Theta1 = theta1;
+ Theta2 = theta2;
+ theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]);
+ theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]);
+
+
+ if (theta1*theta2 < 0.0) {
+ return (FALSE);
+ }
+
+ if (fabs(theta1) < fabs(theta2)) {
+ return (FALSE);
+ }
+
+ return (TRUE);
+
+}
+
+/* check if point x1,y1 is in box formed by x2[0-4] */
+int corner_check (double *x1, double *y1, double *x2, double *y2) {
+
+ int i;
+ double theta;
+
+ theta = 0;
+
+ for (i = 0; i < 4; i++) {
+ theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]);
+ }
+ if (fabs(theta) > 6) {
+ return (TRUE);
+ } else {
+ return (FALSE);
+ }
+}
+
+/* returns the opening angle between the three points (2 is in middle)
+ in range -pi to pi */
+
+double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) {
+
+ double dx1, dy1, dx2, dy2, ct, st, theta;
+
+ dx1 = x1 - x2;
+ dy1 = y1 - y2;
+
+ dx2 = x3 - x2;
+ dy2 = y3 - y2;
+
+ ct = (dx1*dx2 + dy1*dy2);
+ st = (dx1*dy2 - dx2*dy1);
+
+ theta = atan2 (st, ct);
+
+ return (theta);
+
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/find_regions.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/find_regions.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/find_regions.c (revision 2485)
@@ -0,0 +1,102 @@
+# include "relphot.h"
+
+/* given a list of images, find all region files which overlap these images.
+ if a region file does not exist, we should recognize this fact somehow */
+
+GSCRegion *find_regions (Image *image, int Nimage, int *Nregions, GSCRegion *fullregion) {
+
+ int i, nregions;
+ double x, y;
+ GSCRegion *regions;
+ double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec;
+
+ if (VERBOSE) fprintf (stderr, "finding regions: ");
+
+ /* find RA and DEC range of images */
+ MinRA = MinDEC = +361;
+ MaxRA = MaxDEC = -361;
+ for (i = 0; i < Nimage; i++) {
+ for (x = 0.0; x <= 1.0; x+=0.25) {
+ for (y = 0.0; y <= 1.0; y+=0.25) {
+ XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords);
+ MinRA = MIN (ra, MinRA);
+ MaxRA = MAX (ra, MaxRA);
+ MinDEC = MIN (dec, MinDEC);
+ MaxDEC = MAX (dec, MaxDEC);
+ }
+ }
+ }
+
+ /* if images overlap 0,360 boundary, then use -180 to +180 range */
+ if (MaxRA - MinRA > 180) {
+ MinRA = +361;
+ MaxRA = -361;
+ for (i = 0; i < Nimage; i++) {
+ for (x = 0.0; x <= 1.0; x+=0.25) {
+ for (y = 0.0; y <= 1.0; y+=0.25) {
+ XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords);
+ if (ra > 180) ra -= 360;
+ MinRA = MIN (ra, MinRA);
+ MaxRA = MAX (ra, MaxRA);
+ }
+ }
+ }
+ }
+
+ fullregion[0].RA[0] = MinRA;
+ fullregion[0].RA[1] = MaxRA;
+ fullregion[0].DEC[0] = MinDEC;
+ fullregion[0].DEC[1] = MaxDEC;
+
+ /* find all region files containing specified ra, dec range */
+ regions = get_regions (MinRA, MaxRA, MinDEC, MaxDEC, &nregions);
+
+ if (VERBOSE) fprintf (stderr, "%d found\n", nregions);
+
+ *Nregions = nregions;
+ return (regions);
+}
+
+
+getfullregion (Image *image, int Nimage, GSCRegion *fullregion) {
+
+ int i;
+ double x, y;
+ double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec;
+
+ /* find RA and DEC range of images */
+ MinRA = MinDEC = +361;
+ MaxRA = MaxDEC = -361;
+ for (i = 0; i < Nimage; i++) {
+ for (x = 0.0; x <= 1.0; x+=0.25) {
+ for (y = 0.0; y <= 1.0; y+=0.25) {
+ XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords);
+ MinRA = MIN (ra, MinRA);
+ MaxRA = MAX (ra, MaxRA);
+ MinDEC = MIN (dec, MinDEC);
+ MaxDEC = MAX (dec, MaxDEC);
+ }
+ }
+ }
+
+ /* if images overlap 0,360 boundary, then use -180 to +180 range */
+ if (MaxRA - MinRA > 180) {
+ MinRA = +361;
+ MaxRA = -361;
+ for (i = 0; i < Nimage; i++) {
+ for (x = 0.0; x <= 1.0; x+=0.25) {
+ for (y = 0.0; y <= 1.0; y+=0.25) {
+ XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords);
+ if (ra > 180) ra -= 360;
+ MinRA = MIN (ra, MinRA);
+ MaxRA = MAX (ra, MaxRA);
+ }
+ }
+ }
+ }
+
+ fullregion[0].RA[0] = MinRA;
+ fullregion[0].RA[1] = MaxRA;
+ fullregion[0].DEC[0] = MinDEC;
+ fullregion[0].DEC[1] = MaxDEC;
+}
Index: /branches/elixir/Ohana/src/relphot/src/free_catalogs.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/free_catalogs.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/free_catalogs.c (revision 2485)
@@ -0,0 +1,13 @@
+# include "relphot.h"
+
+void free_catalogs (Catalog *catalog, int Ncatalog) {
+
+ int i;
+
+ for (i = 0; i < Ncatalog; i++) {
+ free (catalog[i].secfilt);
+ free (catalog[i].average);
+ free (catalog[i].measure);
+ }
+
+}
Index: /branches/elixir/Ohana/src/relphot/src/gcatalog.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/gcatalog.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/gcatalog.c (revision 2485)
@@ -0,0 +1,26 @@
+# include "relphot.h"
+
+int gcatalog (Catalog *catalog, int FINAL) {
+
+ char mode, status;
+
+ if (FINAL) {
+ mode = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+ } else {
+ mode = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
+ }
+
+ if (!load_catalog (catalog, mode, VERBOSE)) {
+ fprintf (stderr, "ERROR: failure loading catalog %s\n", catalog[0].filename);
+ exit (1);
+ }
+
+ if (VERBOSE) {
+ fprintf (stderr, "read %d stars from catalog file %s (%d measurements)\n",
+ catalog[0].Naverage, catalog[0].filename, catalog[0].Nmeasure);
+ }
+ return (TRUE);
+}
+
+
+
Index: /branches/elixir/Ohana/src/relphot/src/get_regions.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/get_regions.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/get_regions.c (revision 2485)
@@ -0,0 +1,208 @@
+# include "relphot.h"
+double BigDecBounds[] = {0.0, 7.5, 15.0, 22.5, 30.0, 37.5, 45.0,
+ 52.5, 60.0, 67.5, 75.0, 82.5, 90.0,
+ 0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0,
+ -52.5, -60.0, -67.5, -75.0, -82.5, -90.0};
+
+double Areas[] = {2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0,
+ 2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0};
+
+char *DecSections[] = {"N0000", "N0730", "N1500", "N2230", "N3000", "N3730", "N4500",
+ "N5230", "N6000", "N6730", "N7500", "N8230", "weirdness",
+ "S0000", "S0730", "S1500", "S2230", "S3000", "S3730", "S4500",
+ "S5230", "S6000", "S6730", "S7500", "S8230", "weirdness"};
+
+char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500",
+ "n5230", "n6000", "n6730", "n7500", "n8230", "weirdness",
+ "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500",
+ "s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"};
+
+char *disk[] = {"disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1",
+ "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "weirdness",
+ "disk 1", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2",
+ "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "weirdness"};
+
+int NBigRASections [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3, 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3};
+
+int NDecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24,
+ 0, 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48};
+
+/* returns a list of region files within the desired RA, DEC region */
+GSCRegion *get_regions (double minRa, double maxRa, double minDec, double maxDec, int *Nregions) {
+
+ char buffer[28800], temp[50];
+ GSCRegion *regions;
+ FILE *f;
+ double RA0, RA1, DEC0, DEC1;
+ int i, j, NBigDec;
+ int NLINES, done, NREGIONS, nregions;
+
+ f = fopen (GSCFILE, "r");
+ if (f == NULL) {
+ fprintf (stderr, "ERROR: can't find regions file %s\n", GSCFILE);
+ exit (0);
+ }
+
+ NREGIONS = 50;
+ ALLOCATE (regions, GSCRegion, NREGIONS);
+ nregions = 0;
+
+ if ((minDec <= -90) || (maxDec >= 90)) {
+ minRa = 0;
+ maxRa = 360;
+ }
+
+ /* use the pole regions, if near pole */
+ if (maxDec > 86.25) {
+ sprintf (regions[nregions].filename, "n8230/pole.cpt\0");
+ regions[nregions].RA[0] = 0;
+ regions[nregions].RA[1] = 360;
+ regions[nregions].DEC[0] = 86.25;
+ regions[nregions].DEC[1] = 90.0;
+ nregions ++;
+ if (nregions == NREGIONS) {
+ NREGIONS += 50;
+ REALLOCATE (regions, GSCRegion, NREGIONS);
+ }
+ }
+
+ if (minDec > 86.25) {
+ *Nregions = nregions;
+ return (regions);
+ }
+
+ if ((minDec < 0) && (maxDec > 0)) {
+ /* Search Both Sides */
+ NBigDec = 0;
+ } else {
+ /* find large DEC region (directory) */
+ NBigDec = -1;
+ for (i = 0; i < 12; i++) {
+ if ((minDec >= BigDecBounds[i]) && (minDec < BigDecBounds[i+1])) {
+ NBigDec = i;
+ break;
+ }
+ }
+ if (NBigDec < 0) {
+ for (i = 13; i < 24; i++) {
+ if ((maxDec < BigDecBounds[i]) && (maxDec >= BigDecBounds[i+1])) {
+ NBigDec = i;
+ break;
+ }
+ }
+ }
+ }
+ if (NBigDec < 0) {
+ fprintf (stderr, "ERROR: Dec out of range: %f\n", minDec);
+ *Nregions = 0;
+ return ((GSCRegion *) NULL);
+ }
+
+ /* count lines before section */
+ NLINES = 0;
+ for (i = 0; i < NBigDec; i++) {
+ NLINES += NDecLines[i];
+ }
+ fseek (f, 5*2880 + 48*NLINES, SEEK_SET);
+
+ /* should be in this section. if not, there is a problem counting... */
+ /* careful with the 0,360.0 boundary **/
+ done = FALSE;
+ for (j = 0; !done && (NBigDec + j < 25); j++) {
+ fread (buffer, 48*NDecLines[NBigDec + j], 1, f);
+ for (i = 0; (i < NDecLines[NBigDec + j]); i++) {
+ strncpy (temp, &buffer[i*48], 48);
+ temp[49] = 0;
+ hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]);
+ if (RA1 < RA0) RA1 += 360.0;
+ if ((DEC1 > 0) && (minDec < DEC1) && (maxDec > DEC0) && (minRa < RA1) && (maxRa > RA0)) {
+ temp[5] = 0;
+ sprintf (regions[nregions].filename, "%s/%s.cpt\0", Dec2Sections[NBigDec + j], &temp[1]);
+ regions[nregions].RA[0] = RA0;
+ regions[nregions].RA[1] = RA1;
+ regions[nregions].DEC[0] = DEC0;
+ regions[nregions].DEC[1] = DEC1;
+ nregions ++;
+ if (nregions == NREGIONS) {
+ NREGIONS += 50;
+ REALLOCATE (regions, GSCRegion, NREGIONS);
+ }
+ }
+ if ((DEC1 < 0) && (minDec < DEC0) && (maxDec > DEC1) && (minRa < RA1) && (maxRa > RA0)) {
+ temp[5] = 0;
+ sprintf (regions[nregions].filename, "%s/%s.cpt\0", Dec2Sections[NBigDec + j], &temp[1]);
+ regions[nregions].RA[0] = RA0;
+ regions[nregions].RA[1] = RA1;
+ regions[nregions].DEC[0] = DEC0;
+ regions[nregions].DEC[1] = DEC1;
+ nregions ++;
+ if (nregions == NREGIONS) {
+ NREGIONS += 50;
+ REALLOCATE (regions, GSCRegion, NREGIONS);
+ }
+ }
+ if (((DEC1 > 0) && (maxDec <= DEC1)) || ((DEC1 < 0) && (minDec >= DEC1))) {
+ done = TRUE;
+ }
+ }
+ if (done && (minDec < 0) && (maxDec > 0) && (BigDecBounds[NBigDec + j + 1] > 0)) {
+ /* skip remaining north sections, try south sections */
+ /* count lines before section */
+ NLINES = 0;
+ for (i = 0; i < 13; i++) {
+ NLINES += NDecLines[i];
+ }
+ fseek (f, 5*2880 + 48*NLINES, SEEK_SET);
+ done = FALSE;
+ j = 12;
+ }
+ }
+
+ REALLOCATE (regions, GSCRegion, MAX (1,nregions));
+ *Nregions = nregions;
+
+ fclose (f);
+ return (regions);
+
+}
+
+/**********/
+int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) {
+
+ int flag_d0, flag_d1, flag_h0, flag_h1;
+ double tmp;
+
+ *d0 = *h0 = *d1 = *h1 = 0;
+
+ flag_h0 = dparse (h0, 1, string);
+ flag_h1 = dparse (h1, 4, string);
+ flag_d0 = dparse (d0, 7, string);
+ flag_d1 = dparse (d1, 9, string);
+ *h0 *= flag_h0;
+ *h1 *= flag_h1;
+ *d0 *= flag_d0;
+ *d1 *= flag_d1;
+
+ dparse (&tmp, 2, string);
+ *h0 += tmp/60.0;
+ dparse (&tmp, 3, string);
+ *h0 += tmp/3600.0;
+
+ dparse (&tmp, 5, string);
+ *h1 += tmp/60.0;
+ dparse (&tmp, 6, string);
+ *h1 += tmp/3600.0;
+
+ dparse (&tmp, 8, string);
+ *d0 += tmp/60.0;
+
+ dparse (&tmp, 10, string);
+ *d1 += tmp/60.0;
+
+ *h0 *= 15*flag_h0;
+ *h1 *= 15*flag_h1;
+ *d0 *= flag_d0;
+ *d1 *= flag_d1;
+
+ return (TRUE);
+}
Index: /branches/elixir/Ohana/src/relphot/src/global_stats.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/global_stats.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/global_stats.c (revision 2485)
@@ -0,0 +1,51 @@
+# include "relphot.h"
+
+StatType statsImagedM (Catalog *catalog);
+
+StatType statsMosaicdM (Catalog *catalog);
+StatType statsMosaicM (Catalog *catalog);
+StatType statsMosaicN (Catalog *catalog);
+StatType statsMosaicX (Catalog *catalog);
+
+global_stats (Catalog *catalog, int Ncatalog) {
+
+ StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD;
+
+ initstats ("MEAN");
+
+ stN = statsStarN (catalog, Ncatalog);
+ stX = statsStarX (catalog, Ncatalog);
+ stS = statsStarS (catalog, Ncatalog);
+
+ imN = statsImageN (catalog);
+ imX = statsImageX (catalog);
+ imM = statsImageM (catalog);
+ imD = statsImagedM (catalog);
+
+ msN = statsMosaicN (catalog);
+ msM = statsMosaicM (catalog);
+ msD = statsMosaicdM (catalog);
+ msX = statsMosaicX (catalog);
+
+ fprintf (stderr, "STATS median mean sigma min max Nmeas\n");
+ fprintf (stderr, "meas / image: %6.0f %6.1f %6.1f %6.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas);
+ fprintf (stderr, "meas / mosaic: %6.0f %6.1f %6.1f %6.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas);
+ fprintf (stderr, "meas / star: %6.0f %6.1f %6.1f %6.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas);
+
+ fprintf (stderr, "chisq image: %6.1f %6.1f %6.1f %6.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas);
+ fprintf (stderr, "chisq mosaic: %6.1f %6.1f %6.1f %6.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);
+ fprintf (stderr, "chisq star: %6.1f %6.1f %6.1f %6.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas);
+
+ fprintf (stderr, "Mcal image: %6.3f %6.3f %6.3f %6.3f %7.3f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas);
+ fprintf (stderr, "dMcal image: %6.3f %6.3f %6.3f %6.3f %7.3f %6d\n", imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas);
+
+ fprintf (stderr, "Mcal mosaic: %6.3f %6.3f %6.3f %6.3f %7.3f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas);
+ fprintf (stderr, "dMcal mosaic: %6.3f %6.3f %6.3f %6.3f %7.3f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
+
+ fprintf (stderr, "dMrel star: %6.3f %6.3f %6.3f %6.3f %7.3f %6d\n\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas);
+
+ initstats (STATMODE);
+
+}
+
+
Index: /branches/elixir/Ohana/src/relphot/src/initialize.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/initialize.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/initialize.c (revision 2485)
@@ -0,0 +1,56 @@
+# include "relphot.h"
+
+void initialize (int argc, char **argv) {
+
+ int Np;
+
+ /* are these set correctly? */
+ ConfigInit (&argc, argv);
+ args (argc, argv);
+
+ if ((photcode = GetPhotcodebyName (argv[2])) == NULL) {
+ fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[2]);
+ exit (1);
+ }
+ if ((photcode[0].type != PHOT_PRI) && (photcode[0].type != PHOT_SEC)) {
+ fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[2]);
+ exit (1);
+ }
+
+ PhotPrimary = (photcode[0].type == PHOT_PRI);
+ PhotNsec = GetPhotcodeNsecfilt ();
+ PhotSec = GetPhotcodeNsec (photcode[0].code);
+
+ initstats (STATMODE);
+
+ IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
+ STAR_BAD = ID_STAR_POOR | ID_STAR_FEW;
+ MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
+
+ if (SHOW_PARAMS) {
+ fprintf (stderr, "current parameter settings:\n");
+ if (TimeSelect) {
+ fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", sec_to_date (TSTART), sec_to_date (TSTOP));
+ } else {
+ fprintf (stderr, "TimeSelect: FALSE\n");
+ }
+ fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
+ fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n",
+ RELPHOT_GRID_X,
+ RELPHOT_GRID_Y,
+ RELPHOT_GRID_BINNING,
+ (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
+
+ fprintf (stderr, "MAG_LIM %lf\n", MAG_LIM);
+ fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER);
+ fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER);
+ fprintf (stderr, "IMAGE_OFFSET %lf\n", IMAGE_OFFSET);
+ fprintf (stderr, "IMAGE_CATALOG %s\n", ImageCat);
+ fprintf (stderr, "IMAGE_CATALOG_TEMPLATE %s\n", ImageTemplate);
+ fprintf (stderr, "CATALOG_TEMPLATE %s\n", CatTemplate);
+ fprintf (stderr, "GSCFILE %s\n", GSCFILE);
+ fprintf (stderr, "CATDIR %s\n", CATDIR);
+ exit (0);
+ }
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/liststats.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/liststats.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/liststats.c (revision 2485)
@@ -0,0 +1,117 @@
+# include "relphot.h"
+
+enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
+ M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
+
+static int statmode;
+
+initstats (char *mode) {
+
+ statmode = -1;
+ if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
+ if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
+ if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN;
+ if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
+ if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
+ if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
+ if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
+
+ if (statmode == -1) {
+ fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
+ exit (1);
+ }
+}
+
+int liststats (double *value, double *dvalue, int N, StatType *stats) {
+
+ int i, ks, ke, Nm;
+ double Mo, dMo, M, dM, X2, dS, *chi;
+
+ stats[0].Nmeas = N;
+ stats[0].mean = 0;
+ stats[0].sigma = 0;
+ stats[0].error = 0;
+ stats[0].chisq = 0;
+ if (N < 1) return (FALSE);
+
+ sortB (value, dvalue, N);
+ stats[0].median = value[(int)(0.5*N)];
+ stats[0].min = value[0];
+ stats[0].max = value[N-1];
+
+ switch (statmode) {
+ case M_MEDIAN:
+ ks = 0;
+ ke = N;
+ Mo = stats[0].median;
+ Nm = N;
+ goto chisq;
+ break;
+ case M_MEAN:
+ case M_WT_MEAN:
+ ks = 0;
+ ke = N;
+ break;
+ case M_INNER_MEAN:
+ case M_INNER_WTMEAN:
+ case M_CHI_INNER_MEAN:
+ case M_CHI_INNER_WTMEAN:
+ ks = 0.25*N + 0.50;
+ ke = 0.75*N + 0.25;
+ if (N <= 3) {
+ ks = 0;
+ ke = N;
+ }
+ break;
+ }
+
+ if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
+ ALLOCATE (chi, double, N);
+ for (i = 0; i < N; i++) {
+ chi[i] = (value[i] - stats[0].median) / dvalue[i];
+ }
+ sortD (chi, value, dvalue, N);
+ free (chi);
+ }
+
+ /* calculating the per-star offset based on the weighted average */
+ M = dM = Nm = 0;
+ if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
+ for (i = ks; i < ke; i++) {
+ M += value[i] / SQ (dvalue[i]);
+ dM += 1.0 / SQ (dvalue[i]);
+ Nm ++;
+ }
+ Mo = M / dM;
+ dMo = sqrt (1.0 / dM);
+ } else {
+ for (i = ks; i < ke; i++) {
+ M += value[i];
+ dM += SQ (dvalue[i]);
+ Nm ++;
+ }
+ Mo = M / (double) Nm;
+ dMo = sqrt (dM / (double) Nm);
+ }
+
+ chisq:
+ /* find sigma and chisq */
+ X2 = dS = 0;
+ for (i = ks; i < ke; i++) {
+ M = SQ (value[i] - Mo);
+ dM = SQ (dvalue[i]);
+ X2 += M / dM;
+ dS += M;
+ }
+ X2 = X2 / Nm;
+ dS = sqrt (dS / Nm);
+
+ stats[0].mean = Mo;
+ stats[0].Nmeas = Nm;
+ stats[0].chisq = X2;
+ stats[0].sigma = dS;
+ stats[0].error = dMo;
+
+ return (TRUE);
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/load_catalogs.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/load_catalogs.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/load_catalogs.c (revision 2485)
@@ -0,0 +1,45 @@
+# include "relphot.h"
+
+Catalog *load_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion) {
+
+ int i, Nstar;
+ char filename[256];
+ Catalog *catalog, tcatalog;
+
+ if (VERBOSE) fprintf (stderr, "loading catalog data\n");
+
+ ALLOCATE (catalog, Catalog, Nregion);
+
+ /* load data from each region file, only use bright stars */
+ for (i = 0; i < Nregion; i++) {
+ sprintf (filename, "%s/%s", CATDIR, region[i].filename);
+ tcatalog.filename = filename; /* don't free region before catalog! */
+ switch (lock_catalog (&tcatalog, LCK_SOFT)) {
+ case 0:
+ fprintf (stderr, "ERROR: can't lock file %s\n", region[i].filename);
+ exit (1);
+ case 1:
+ gcatalog (&tcatalog, FALSE); /* load from disk */
+ break;
+ case 2:
+ if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
+ tcatalog.Naverage = 0;
+ tcatalog.Nmeasure = 0;
+ break;
+ default:
+ fprintf (stderr, "weird lock_catalog exit state\n");
+ exit (1);
+ }
+ bcatalog (&catalog[i], &tcatalog, fullregion);
+ }
+
+ Nstar = 0;
+ for (i = 0; i < Nregion; i++) {
+ Nstar += catalog[i].Naverage;
+ }
+ if (Nstar < 2) {
+ fprintf (stderr, "insufficient stars %d\n", Nstar);
+ exit (0);
+ }
+ return (catalog);
+}
Index: /branches/elixir/Ohana/src/relphot/src/load_images.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/load_images.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/load_images.c (revision 2485)
@@ -0,0 +1,67 @@
+# include "relphot.h"
+
+static int Nimage;
+static int *Line;
+static Image *image;
+
+GSCRegion *load_images (char *seed, int *nregion, GSCRegion *fullregion) {
+
+ GSCRegion *region;
+ int Nregion, Ntimage;
+ Image *timage;
+
+ region = name_region (seed, &Nregion);
+
+ /* find images which supplied all measurements */
+ timage = find_images (region, Nregion, &Ntimage, &Line);
+ free (Line);
+
+ /* find all region files which touch images */
+ region = find_regions (timage, Ntimage, &Nregion, fullregion);
+ free (timage);
+
+ /* find images which supplied all measurements for all regions */
+ image = find_images (region, Nregion, &Nimage, &Line);
+ getfullregion (image, Nimage, fullregion);
+
+ initImages (image, Nimage);
+ initMosaics (image, Nimage);
+
+ *nregion = Nregion;
+ return (region);
+
+}
+
+void wimages () {
+
+ int i, status, mode;
+ FILE *f;
+ Header header;
+
+ setMcalFinal (image, Nimage);
+
+ if (!fits_read_header (ImageCat, &header)) {
+ fprintf (stderr, "ERROR: can't find %s\n", ImageCat);
+ exit (0);
+ }
+
+ f = fopen (ImageCat, "r+");
+ if (f == NULL) {
+ fprintf (stderr, "ERROR: can't open image catalog file to append: %s\n", ImageCat);
+ exit (0);
+ }
+ mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
+ chmod (ImageCat, mode);
+
+ if (VERBOSE) fprintf (stderr, "writing out %d images\n", Nimage);
+ for (i = 0; i < Nimage; i++) {
+ fseek (f, header.size + Line[i]*sizeof(Image), SEEK_SET);
+ status = Fwrite (&image[i], sizeof(Image), 1, f, "image");
+ if (status != 1) {
+ fprintf (stderr, "ERROR: failure writing image data to catalog %s\n", ImageCat);
+ exit (0);
+ }
+ }
+
+ fclose (f);
+}
Index: /branches/elixir/Ohana/src/relphot/src/misc.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/misc.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/misc.c (revision 2485)
@@ -0,0 +1,182 @@
+
+sortA (double *X, int N) {
+
+ int l,j,ir,i;
+ double tX;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ }
+}
+
+sortB (double *X, double *Y, int N) {
+
+ int l,j,ir,i;
+ double tX, tY;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ tY = Y[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ tY = Y[ir];
+ Y[ir] = Y[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ Y[0] = tY;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ Y[i] = Y[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ Y[i] = tY;
+ }
+}
+
+sortC (double *X, double *Y, double *F1, double *F2, int N) {
+
+ int l,j,ir,i;
+ double tX, tY;
+ double t1, t2;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ tY = Y[l];
+ t1 = F1[l];
+ t2 = F2[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ tY = Y[ir];
+ Y[ir] = Y[0];
+ t1 = F1[ir];
+ F1[ir] = F1[0];
+ t2 = F2[ir];
+ F2[ir] = F2[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ Y[0] = tY;
+ F1[0] = t1;
+ F2[0] = t2;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ Y[i] = Y[j];
+ F1[i] = F1[j];
+ F2[i] = F2[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ Y[i] = tY;
+ F1[i] = t1;
+ F2[i] = t2;
+ }
+}
+
+sortD (double *X, double *Y, double *Z, int N) {
+
+ int l,j,ir,i;
+ double tX, tY, tZ;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ tY = Y[l];
+ tZ = Z[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ tY = Y[ir];
+ Y[ir] = Y[0];
+ tZ = Z[ir];
+ Z[ir] = Z[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ Y[0] = tY;
+ Z[0] = tZ;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ Y[i] = Y[j];
+ Z[i] = Z[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ Y[i] = tY;
+ Z[i] = tZ;
+ }
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/name_region.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/name_region.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/name_region.c (revision 2485)
@@ -0,0 +1,24 @@
+# include "relphot.h"
+
+GSCRegion *name_region (char *name, int *Nregions) {
+
+ GSCRegion *region;
+ Header header;
+
+ ALLOCATE (region, GSCRegion, 1);
+
+ sprintf (region[0].filename, "%s/%s", CATDIR, name);
+ if (!fits_read_header (region[0].filename, &header)) {
+ if (VERBOSE) fprintf (stderr, "catalog file not yet available: %s\n", region[0].filename);
+ exit (0);
+ }
+ fits_scan (&header, "RA0", "%lf", 1, ®ion[0].RA[0]);
+ fits_scan (&header, "RA1", "%lf", 1, ®ion[0].RA[1]);
+ fits_scan (&header, "DEC0", "%lf", 1, ®ion[0].DEC[0]);
+ fits_scan (&header, "DEC1", "%lf", 1, ®ion[0].DEC[1]);
+ fits_free_header (&header);
+
+ *Nregions = 1;
+ return (region);
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/plot_scatter.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/plot_scatter.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/plot_scatter.c (revision 2485)
@@ -0,0 +1,61 @@
+# include "relphot.h"
+
+plot_scatter (Catalog *catalog, int Ncatalog) {
+
+ char c;
+ int i, j, k, m, N, Ntot;
+ short Mrel, Mcal, Mmos, Mgrid;
+ double *xlist, *ylist, *ilist;
+ StatType stats;
+ Graphdata graphdata;
+
+ Ntot = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ Ntot += catalog[i].average[j].Nm;
+ }
+ }
+ ALLOCATE (xlist, double, Ntot);
+ ALLOCATE (ylist, double, Ntot);
+ ALLOCATE (ilist, double, Ntot);
+
+ N = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+
+ /* calculate the average value for a single star */
+ if (catalog[i].average[j].code & STAR_BAD) continue;
+ m = catalog[i].average[j].offset;
+
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+ if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+ if ((Mcal = getMcal (m, i)) == NO_MAG) continue;
+ if ((Mmos = getMmos (m, i)) == NO_MAG) continue;
+ if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+
+ Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M;
+ xlist[N] = Mrel;
+ ylist[N] = iPhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;
+ ilist[N] = PhotInst (&catalog[i].measure[m]);
+ xlist[N] *= 0.001;
+ ylist[N] *= 0.001;
+ N++;
+ }
+ }
+ }
+
+ plot_defaults (&graphdata);
+ graphdata.xmin = PlotMmin;
+ graphdata.xmax = PlotMmax;
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "Mag.png");
+
+ plot_defaults (&graphdata);
+ graphdata.ymin = PlotdMmin;
+ graphdata.ymax = PlotdMmax;
+ plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", NULL);
+ free (xlist);
+ free (ylist);
+ free (ilist);
+}
Index: /branches/elixir/Ohana/src/relphot/src/plotstuff.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/plotstuff.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/plotstuff.c (revision 2485)
@@ -0,0 +1,224 @@
+# include "relphot.h"
+# include
+
+static int Xgraph[5] = {0,0,0,0,0};
+static int active;
+
+void XDead () {
+ signal (SIGPIPE, XDead);
+ fprintf (stderr, "kapa is dead, must restart\n");
+ Xgraph[active] = -1;
+}
+
+int open_graph (int N) {
+
+# ifdef ANSI
+# define F_SETFL 4
+# define O_NONBLOCK 0200000
+# define AF_UNIX 1
+# define SOCK_STREAM 1
+#define ENOENT 2 /* No such file or directory */
+# endif /* ANSI */
+
+ int i, InitSocket, status, addreslen;
+ struct sockaddr_un Address;
+ char temp[100], *display_name;
+ char socket_name[100];
+ FILE *f;
+
+ active = N;
+ sprintf (socket_name, "/tmp/kapa%d", N);
+ sprintf (temp, "rm -f %s", socket_name);
+ system (temp);
+
+ strcpy (Address.sun_path, socket_name);
+ Address.sun_family = AF_UNIX;
+ InitSocket = socket (AF_UNIX, SOCK_STREAM, 0);
+ status = bind (InitSocket, &Address, sizeof (Address));
+ status = listen (InitSocket, 1);
+
+ sprintf (temp, "kapa %s &", socket_name);
+# ifndef DEBUG
+ system (temp);
+# else
+ fprintf (stderr, "start kapa, press return: %s\n", temp);
+ fscanf (stdin, "%d", &i);
+# endif
+
+ addreslen = sizeof (Address);
+ Xgraph[N] = accept (InitSocket, &Address, &addreslen);
+ if (Xgraph[N] < 0) {
+ fprintf (stderr, "error starting kapa\n");
+ return (FALSE);
+ }
+ else {
+ return (TRUE);
+ }
+
+}
+
+void DonePlotting (Graphdata *graphmode, int N) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "DBOX", 4);
+ sprintf (buffer, "%f %f %f %f", graphmode[0].xmin, graphmode[0].xmax, graphmode[0].ymin, graphmode[0].ymax);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+ sprintf (buffer, "%s %s %s", "2222", "2222", "2222");
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+}
+
+void JpegPlot (Graphdata *graphmode, int N, char *filename) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "PNGF", 4);
+ sprintf (buffer, "LEN: %11d", strlen (filename));
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], filename, strlen (filename));
+ read (Xgraph[N], buffer, 4);
+}
+
+void PSPlot (Graphdata *graphmode, int N, char *filename) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "PSIT", 4);
+ sprintf (buffer, "LEN: %11d", strlen (filename));
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], filename, strlen (filename));
+ read (Xgraph[N], buffer, 4);
+}
+
+void PrepPlotting (int Npts, Graphdata *graphmode, int N) {
+
+ int i, status;
+ char buffer[128], buffer2[128];
+
+ active = N;
+ if (Npts < 1) return;
+
+ if (Xgraph[N] < 1) if (!open_graph(N)) return;
+
+ /* test Xgraph[N], flush junk from pipe */
+ signal (SIGPIPE, XDead);
+ fcntl (Xgraph[N], F_SETFL, O_NONBLOCK);
+ for (i = 0; (read (Xgraph[N], buffer, 64) > 0) && (i < 20); i++);
+ fcntl (Xgraph[N], F_SETFL, !O_NONBLOCK);
+
+ write (Xgraph[N], "ERAS", 4);
+
+ /* tell kapa to look for the incoming image */
+ status = write (Xgraph[N], "PLOT", 4);
+
+ /* send Xgraph[N] the plot details */
+ sprintf (buffer, "%8d %8d %d %d %d %d %f %f",
+ Npts, graphmode[0].style,
+ graphmode[0].ptype, graphmode[0].ltype,
+ graphmode[0].etype, graphmode[0].color,
+ graphmode[0].lweight, graphmode[0].size);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+ sprintf (buffer, "%f %f %f %f",
+ graphmode[0].xmin, graphmode[0].xmax,
+ graphmode[0].ymin, graphmode[0].ymax);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+}
+
+void PlotLabel (char *string, int N) {
+
+ char buffer[32];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "LABL", 4);
+ sprintf (buffer, " %6d %6d", strlen (string), 2);
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], string, strlen (string));
+}
+
+void PlotVector (int Npts, double *vect, int mode, int N) {
+
+ float *values;
+ int i, Nbytes;
+
+ active = N;
+ if (Npts < 1) return;
+
+ ALLOCATE (values, float, Npts);
+ for (i = 0; i < Npts; i++) {
+ values[i] = vect[i];
+ }
+
+ Nbytes = Npts * sizeof (float);
+ write (Xgraph[N], values, Nbytes);
+ free (values);
+
+}
+
+plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file) {
+
+ int i;
+ char c;
+ StatType stats;
+
+ stats.min = stats.max = xlist[0];
+ for (i = 0; i < N; i++) {
+ stats.min = MIN (stats.min, xlist[i]);
+ stats.max = MAX (stats.max, xlist[i]);
+ }
+ if (graphdata[0].xmin == dUNDEF) graphdata[0].xmin = 1.05*stats.min - 0.05*stats.max;
+ if (graphdata[0].xmax == dUNDEF) graphdata[0].xmax = 1.05*stats.max - 0.05*stats.min;
+
+ stats.min = stats.max = ylist[0];
+ for (i = 0; i < N; i++) {
+ stats.min = MIN (stats.min, ylist[i]);
+ stats.max = MAX (stats.max, ylist[i]);
+ }
+ if (graphdata[0].ymin == dUNDEF) graphdata[0].ymin = 1.05*stats.min - 0.05*stats.max;
+ if (graphdata[0].ymax == dUNDEF) graphdata[0].ymax = 1.05*stats.max - 0.05*stats.min;
+
+ PrepPlotting (N, graphdata, 0);
+ PlotVector (N, xlist, 0, 0);
+ PlotVector (N, ylist, 1, 0);
+ if (label != NULL) PlotLabel (label, 0);
+ DonePlotting (graphdata, 0);
+
+ if ((file != NULL) && SAVEPLOT) JpegPlot (graphdata, 0, file);
+ if (PLOTDELAY > 500000) {
+ fprintf (stdout, "press return\n");
+ fscanf (stdin, "%*c");
+ } else {
+ usleep (PLOTDELAY);
+ }
+}
+
+plot_defaults (Graphdata *graphdata) {
+
+ graphdata[0].style = 2;
+ graphdata[0].ptype = 2;
+ graphdata[0].ltype = 0;
+ graphdata[0].etype = 0;
+ graphdata[0].color = black;
+ graphdata[0].lweight = 0;
+ graphdata[0].size = 0.5;
+
+ graphdata[0].xmin = dUNDEF;
+ graphdata[0].xmax = dUNDEF;
+ graphdata[0].ymin = dUNDEF;
+ graphdata[0].ymax = dUNDEF;
+
+}
Index: /branches/elixir/Ohana/src/relphot/src/reload_catalogs.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/reload_catalogs.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/reload_catalogs.c (revision 2485)
@@ -0,0 +1,46 @@
+# include "relphot.h"
+
+void reload_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion) {
+
+ int i;
+ char filename[256];
+ Catalog catalog;
+
+ if (VERBOSE) fprintf (stderr, "re-loading catalog data\n");
+
+ /* load data from each region file */
+ for (i = 0; i < Nregion; i++) {
+ sprintf (filename, "%s/%s", CATDIR, region[i].filename);
+ catalog.filename = filename; /* don't free region before catalog! */
+ switch (lock_catalog (&catalog, LCK_XCLD)) {
+ case 0:
+ fprintf (stderr, "ERROR: can't lock file %s\n", region[i].filename);
+ exit (1);
+ case 1:
+ gcatalog (&catalog, TRUE); /* load from disk */
+ break;
+ case 2:
+ if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+ unlock_catalog (&catalog);
+ continue;
+ break;
+ default:
+ fprintf (stderr, "weird lock_catalog exit state\n");
+ exit (1);
+ }
+ initImageBins (&catalog, 1);
+ initMosaicBins (&catalog, 1);
+ initGridBins (&catalog, 1);
+
+ findImages (&catalog, 1);
+ findMosaics (&catalog, 1);
+
+ setMrelFinal (&catalog);
+ wcatalog (&catalog);
+ unlock_catalog (&catalog);
+
+ freeImageBins (1);
+ freeMosaicBins (1);
+ freeGridBins (1);
+ }
+}
Index: /branches/elixir/Ohana/src/relphot/src/relphot.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/relphot.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/relphot.c (revision 2485)
@@ -0,0 +1,85 @@
+# include "relphot.h"
+/** I'm not currently setting the lockfiles. this should be fixed! **/
+
+int main (int argc, char **argv) {
+
+ char c;
+ int i, Ncatalog;
+ Catalog *catalog, tcatalog;
+ GSCRegion fullregion, *region;
+
+ /* get configuration info, args */
+ initialize (argc, argv);
+
+ /* load regions and images based on seed region */
+ region = load_images (argv[1], &Ncatalog, &fullregion);
+
+ /* load catalog data from region files */
+ catalog = load_catalogs (region, Ncatalog, &fullregion);
+
+ /* match measurements with images, mosaics */
+ initImageBins (catalog, Ncatalog);
+ initMosaicBins (catalog, Ncatalog);
+ initGridBins (catalog, Ncatalog);
+ initMrel (catalog, Ncatalog);
+
+ findImages (catalog, Ncatalog);
+ findMosaics (catalog, Ncatalog); /* also sets Grid values */
+ SAVEPLOT = FALSE;
+
+ setExclusions (catalog, Ncatalog);
+
+ if (PLOTSTUFF) {
+ plot_star_coords (catalog, Ncatalog);
+ plot_mosaic_fields (catalog);
+ }
+
+ /* determine fit values */
+ for (i = 0; i < NLOOP; i++) {
+ setMrel (catalog, Ncatalog);
+ setMcal (catalog, FALSE);
+ setMmos (catalog, FALSE);
+ setMgrid (catalog);
+
+ if (PLOTSTUFF) {
+ plot_scatter (catalog, Ncatalog);
+ plot_grid (catalog);
+ plot_mosaics ();
+ plot_images ();
+ plot_stars (catalog, Ncatalog);
+ plot_chisq (catalog, Ncatalog);
+ }
+ if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
+ if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
+ if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (catalog, Ncatalog);
+ if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
+ global_stats (catalog, Ncatalog);
+ }
+
+ SAVEPLOT = TRUE;
+ plot_scatter (catalog, Ncatalog);
+ plot_grid (catalog);
+ plot_mosaics ();
+ plot_images ();
+ plot_stars (catalog, Ncatalog);
+ plot_chisq (catalog, Ncatalog);
+
+ if (USE_GRID) dump_grid ();
+ if (!UPDATE) exit (0);
+
+ /* set Mcal & Mmos for bad images */
+ setMcal (catalog, TRUE);
+ setMmos (catalog, TRUE);
+
+ /* at this point, we have correct cal coeffs in the image/mosaic structures */
+ free_catalogs (catalog, Ncatalog);
+ freeImageBins (Ncatalog);
+ freeMosaicBins (Ncatalog);
+ freeGridBins (Ncatalog);
+
+ /* load catalog data from region files, update Mrel include all data */
+ reload_catalogs (region, Ncatalog, &fullregion);
+ wimages ();
+
+}
+
Index: /branches/elixir/Ohana/src/relphot/src/setExclusions.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/setExclusions.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/setExclusions.c (revision 2485)
@@ -0,0 +1,56 @@
+# include "relphot.h"
+
+Coords *getCoords (int meas, int cat);
+
+int setExclusions (Catalog *catalog, int Ncatalog) {
+
+ int i, j, k, m, Narea, Nnocal, ecode;
+ Coords *coords;
+ double r, d, x, y;
+
+ Nnocal = Narea = 0;
+ for (i = 0; i < Ncatalog; i++) {
+ for (j = 0; j < catalog[i].Naverage; j++) {
+ m = catalog[i].average[j].offset;
+ for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+
+ /* select measurements by photcode */
+ ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].source);
+ if (ecode != photcode[0].code) goto mark_nocal;
+
+ /* select measurements by time */
+ if (TimeSelect) {
+ if (catalog[i].measure[m].t < TSTART) goto mark_nocal;
+ if (catalog[i].measure[m].t > TSTOP) goto mark_nocal;
+ }
+
+ /* select measurements by mag limit */
+ if (AreaSelect) {
+ r = catalog[i].average[j].R + catalog[i].measure[m].dR / 360000.0;
+ d = catalog[i].average[j].D + catalog[i].measure[m].dD / 360000.0;
+ if ((coords = getCoords (m, i)) == NULL) goto markbad;
+ RD_to_XY (&x, &y, r, d, coords);
+ if (x < AreaXmin) goto markbad;
+ if (x > AreaXmax) goto markbad;
+ if (y < AreaYmin) goto markbad;
+ if (y > AreaYmax) goto markbad;
+ }
+ continue;
+
+ markbad:
+ catalog[i].measure[m].flags |= ID_MEAS_AREA;
+ Narea ++;
+ continue;
+
+ mark_nocal:
+ catalog[i].measure[m].flags |= ID_MEAS_NOCAL;
+ Nnocal ++;
+ continue;
+ }
+ }
+ }
+ if (VERBOSE) fprintf (stderr, "%d measurements marked by area\n", Narea);
+ if (VERBOSE) fprintf (stderr, "%d measurements marked nocal\n", Nnocal);
+}
+
+
Index: /branches/elixir/Ohana/src/relphot/src/setMrelFinal.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/setMrelFinal.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/setMrelFinal.c (revision 2485)
@@ -0,0 +1,157 @@
+# include "relphot.h"
+
+void skip_measurements (Catalog *catalog, int pass);
+Image *getimages (int *N);
+
+void setMrelFinal (Catalog *catalog) {
+
+ int i, j, m, ecode;
+ short *p;
+
+ /* if we reset the catalog, reset all the current measurements */
+ if (RESET) {
+
+ for (i = 0; i < catalog[0].Naverage; i++) {
+ p = (PhotPrimary) ? &catalog[0].average[i].M : &catalog[0].secfilt[PhotNsec*i+PhotSec].M;
+ *p = NO_MAG;
+ p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM;
+ *p = NO_MAG;
+ p = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm;
+ *p = NO_MAG;
+
+ m = catalog[0].average[i].offset;
+ for (j = 0; j < catalog[0].average[i].Nm; j++, m++) {
+
+ /* select measurements by photcode */
+ ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
+ if (ecode != photcode[0].code) continue;
+
+ /* select measurements by time */
+ if (TimeSelect) {
+ if (catalog[0].measure[m].t < TSTART) continue;
+ if (catalog[0].measure[m].t > TSTOP) continue;
+ }
+
+ catalog[0].measure[m].Mcal = 0;
+ catalog[0].measure[m].flags &= 0xff00;
+ catalog[0].measure[m].flags &= ~ID_MEAS_POOR;
+ catalog[0].measure[m].flags &= ~ID_MEAS_SKIP;
+ catalog[0].measure[m].flags &= ~ID_MEAS_AREA;
+ catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL;
+ }
+ }
+ }
+
+ setExclusions (catalog, 1); /* mark by area */
+
+ /* set catalog[0].found[i] = FALSE */
+ ALLOCATE (catalog[0].found, int, MAX (1, catalog[0].Naverage));
+ for (i = 0; i < catalog[0].Naverage; i++) {
+ catalog[0].found[i] = FALSE;
+ }
+
+ for (i = 0; i < 5; i++) {
+ skip_measurements (catalog, i); /* set ID_MEAS_SKIP for measures to be skipped */
+ setMrelOutput (catalog, 1, FALSE); /* set Mrel using all measures not skipped */
+ clean_measures (catalog, 1, TRUE); /* mark outliers ID_MEAS_POOR */
+ setMrelOutput (catalog, 1, TRUE); /* set Mrel using remaining measures */
+ }
+ setMcalOutput (catalog, 1);
+
+ /* clear ID_STAR_POOR, ID_STAR_FEW, ID_MEAS_NOCAL values before writing ??? */
+ for (i = 0; i < catalog[0].Naverage; i++) {
+ catalog[0].average[i].code &= ~ID_STAR_FEW;
+ catalog[0].average[i].code &= ~ID_STAR_POOR;
+ m = catalog[0].average[i].offset;
+ for (j = 0; j < catalog[0].average[i].Nm; j++, m++) {
+ catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL;
+ }
+ }
+}
+
+/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
+void skip_measurements (Catalog *catalog, int pass) {
+
+ int i, j, k, m, Nimage, ecode, d1, d2;
+ int Ntot, Ntry, Nkeep, Nskip;
+ short mag;
+ Image *image;
+
+ Ntot = Ntry = Nskip = Nkeep = 0;
+
+ /* allow measures from images marked POOR and FEW */
+ if (pass >= 3) IMAGE_BAD = ID_IMAGE_NOCAL;
+
+ /* allow measures marked as outliers (POOR) and off image region (AREA) */
+ if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP;
+
+ /* mark measures which should be ignored on second pass */
+ for (i = 0; i < catalog[0].Naverage; i++) {
+ Ntot += catalog[0].average[i].Nm;
+ if (catalog[0].found[i]) continue;
+
+ m = catalog[0].average[i].offset;
+ for (k = 0; k < catalog[0].average[i].Nm; k++, m++) {
+ Ntry++;
+
+ /* clear SKIP for all measures at first */
+ catalog[0].measure[m].flags &= ~ID_MEAS_SKIP;
+
+ /** never use these measurements (wrong photcode, bad time range) */
+ /* skipped via NOCAL, don't mark as skipped */
+ ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
+ if (ecode != photcode[0].code) continue;
+
+ /* skip measurements by time range */
+ if (TimeSelect) {
+ if (catalog[0].measure[m].t < TSTART) goto skip;
+ if (catalog[0].measure[m].t > TSTOP) goto skip;
+ }
+
+ /* skip measurements with sat. dophot values */
+ if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip;
+
+ /* skip measurements from BAD images and mosaics */
+ if ((getMcal (m, 0)) == NO_MAG) goto skip;
+ if ((getMmos (m, 0)) == NO_MAG) goto skip;
+
+ /* skip measurements by inst mag limit */
+ if ((pass < 4) && ImagSelect) {
+ mag = iPhotInst (&catalog[0].measure[m]);
+ if (mag < ImagMin) goto skip;
+ if (mag > ImagMax) goto skip;
+ }
+
+ d1 = (catalog[0].measure[m].dophot == 1);
+ d2 = (catalog[0].measure[m].dophot == 2);
+
+ /* skip measurements by measurement error */
+ if ((pass < 2) && !(d1 || d2)) goto skip;
+
+ /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */
+ if ((pass < 1) && !d1) goto skip;
+ Nkeep ++;
+
+ continue;
+
+ skip:
+ catalog[0].measure[m].flags |= ID_MEAS_SKIP;
+ Nskip ++;
+ }
+ }
+ fprintf (stderr, "pass %d, Ntot: %d, Ntry: %d, Nskip: %d, Nkeep: %d\n",
+ pass, Ntot, Ntry, Nskip, Nkeep);
+}
+
+/* dophot types:
+
+ 1 - star
+ 2 - galaxy
+ 3 - star
+ 4 - weak
+ 5 - edge
+ 7 - weak
+ 9 - weak
+ 10 - sat
+*/
+
Index: /branches/elixir/Ohana/src/relphot/src/time.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/time.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/time.c (revision 2485)
@@ -0,0 +1,321 @@
+# include
+
+int dms_to_ddd (double *Value, char *string);
+int str_to_radec (double *ra, double *dec, char *str1, char *str2);
+int chk_time (char *line);
+double sec_to_jd (unsigned long second);
+unsigned long int jd_to_sec (double jd);
+char *sec_to_date (unsigned long second);
+unsigned long date_to_sec (char *date);
+
+/***** convert [-]00:00:00 to 0.0000 ****/
+int dms_to_ddd (double *Value, char *string) {
+
+ int valid, neg, status;
+ double tmp, value;
+ char *p1, *p2, *px;
+
+ valid = FALSE;
+ neg = FALSE;
+ stripwhite (string);
+ p1 = string;
+ px = string + strlen(string);
+
+ if (string[0] == '-') {
+ valid = TRUE;
+ neg = TRUE;
+ p1 = &string[1];
+ }
+ if (string[0] == '+') {
+ valid = TRUE;
+ neg = FALSE;
+ p1 = &string[1];
+ }
+ if (isdigit(string[0])) {
+ valid = TRUE;
+ p1 = &string[0];
+ }
+ if (!valid) { return (FALSE); }
+
+ status = 1;
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) return (FALSE); /* entry not a number: +fred */
+ value = tmp;
+ if (p2 == px) goto escape; /* entry only number: +1.0 */
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape; /* entry not a number: +1:fred */
+ status = 2;
+ value += tmp / 60.0;
+ if (p2 == px) goto escape; /* entry only number: +1:1 */
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */
+ value += tmp / 3600.0;
+
+ escape:
+ if (neg) {
+ value *= -1;
+ }
+ *Value = value;
+
+ return (status);
+}
+
+/**********/
+int str_to_radec (double *ra, double *dec, char *str1, char *str2) {
+
+ double Ra, Dec;
+
+ *ra = *dec = 0;
+ switch (dms_to_ddd (&Ra, str1)) {
+ case 0:
+ fprintf (stderr, "syntax error in RA\n");
+ return (FALSE);
+ case 1:
+ break;
+ case 2:
+ Ra = Ra * 15;
+ break;
+ }
+ switch (dms_to_ddd (&Dec, str2)) {
+ case 0:
+ fprintf (stderr, "syntax error in DEC\n");
+ return (FALSE);
+ case 1:
+ case 2:
+ break;
+ }
+ *ra = Ra;
+ *dec = Dec;
+ return (TRUE);
+}
+
+# define FORMAT_DAYS 1
+# define FORMAT_HOURS 2
+# define FORMAT_MINUTES 3
+# define FORMAT_SECONDS 4
+# define FORMAT_JD 5
+# define FORMAT_DATE 6
+
+/**********/
+int chk_time (char *line) {
+
+ char *p1, *p2;
+ double tmp;
+ int mode;
+
+ p1 = line;
+ tmp = strtod (p1, &p2);
+ if (p2 == p1 + strlen (p1) - 1) {
+ if (*p2 == 'd') {
+ mode = FORMAT_DAYS;
+ }
+ if (*p2 == 'h') {
+ mode = FORMAT_HOURS;
+ }
+ if (*p2 == 'm') {
+ mode = FORMAT_MINUTES;
+ }
+ if (*p2 == 's') {
+ mode = FORMAT_SECONDS;
+ }
+ if (*p2 == 'j') {
+ mode = FORMAT_JD;
+ }
+ } else {
+ mode = FORMAT_DATE;
+ }
+ return (mode);
+}
+
+/**********/
+int str_to_time (char *line, unsigned int *second) {
+
+ struct timeval now;
+ double jd;
+
+ if (!strcasecmp (line, "NOW")) {
+ gettimeofday (&now, (struct timezone *) NULL);
+ *second = now.tv_sec;
+ return (TRUE);
+ }
+
+ if (!strcasecmp (line, "TODAY")) {
+ gettimeofday (&now, (struct timezone *) NULL);
+ *second = 86400 * ((int)(now.tv_sec / 86400));
+ return (TRUE);
+ }
+
+ switch (chk_time (line)) {
+ case 0:
+ return (FALSE);
+ case FORMAT_DAYS:
+ *second = strtod (line, 0) * 86400.0;
+ return (TRUE);
+ case FORMAT_HOURS:
+ *second = strtod (line, 0) * 3600.0;
+ return (TRUE);
+ case FORMAT_MINUTES:
+ *second = strtod (line, 0) * 60.0;
+ return (TRUE);
+ case FORMAT_SECONDS:
+ *second = strtod (line, 0);
+ return (TRUE);
+ case FORMAT_JD:
+ jd = strtod (line, 0);
+ *second = jd_to_sec (jd);
+ return (TRUE);
+ case FORMAT_DATE:
+ *second = date_to_sec (line);
+ return (TRUE);
+ }
+ return (FALSE);
+}
+
+/**********/
+int str_to_dtime (char *line, double *second) {
+
+ switch (chk_time (line)) {
+ case 0:
+ case FORMAT_JD:
+ case FORMAT_DATE:
+ return (FALSE);
+ case FORMAT_DAYS:
+ *second = strtod (line, 0) * 86400.0;
+ return (TRUE);
+ case FORMAT_HOURS:
+ *second = strtod (line, 0) * 3600.0;
+ return (TRUE);
+ case FORMAT_MINUTES:
+ *second = strtod (line, 0) * 60.0;
+ return (TRUE);
+ case FORMAT_SECONDS:
+ *second = strtod (line, 0);
+ return (TRUE);
+ }
+ return (FALSE);
+}
+
+/**********/
+double sec_to_jd (unsigned long second) {
+
+ double jd;
+
+ jd = second/86400.0 + 2440587.5;
+ return (jd);
+}
+
+/**********/
+unsigned long int jd_to_sec (double jd) {
+
+ unsigned long int second;
+
+ second = (jd - 2440587.5)*86400;
+ return (second);
+}
+
+/**********/
+char *sec_to_date (unsigned long second) {
+
+ struct tm *gmt;
+ char *line;
+
+ ALLOCATE (line, char, 64);
+ gmt = gmtime (&second);
+ sprintf (line, "%4d/%02d/%02d,%02d:%02d:%02d", 1900+gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
+ return (line);
+
+}
+
+/***** date in format yyyy/mm/dd,hh:mm:ss *****/
+unsigned long date_to_sec (char *date) {
+
+ unsigned long second;
+ double tmp, jd;
+ struct tm now;
+ char *p1, *p2, *px;
+
+ p1 = date;
+ px = date + strlen(date);
+ bzero (&now, sizeof(now));
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ now.tm_year = tmp;
+ if (now.tm_year > 1000) now.tm_year -= 1900;
+ if (now.tm_year < 50) now.tm_year += 100;
+ if (p2 == px) goto escape;
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ now.tm_mon = tmp - 1; /* mon runs from 0 - 11 */
+ if (p2 == px) goto escape;
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ now.tm_mday = tmp;
+ if (p2 == px) goto escape;
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ p1 = p2 + 1;
+ now.tm_hour = tmp;
+ if (p2 == px) goto escape;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ now.tm_min = tmp;
+ if (p2 == px) goto escape;
+ p1 = p2 + 1;
+
+ tmp = strtod (p1, &p2);
+ if (p2 == p1) goto escape;
+ now.tm_sec = tmp;
+ if (p2 == px) goto escape;
+ p1 = p2 + 1;
+
+ escape:
+ jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4)
+ + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12)
+ - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5;
+
+ second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec;
+
+ return (second);
+}
+
+
+/* times may be in forms as:
+ 20040200450s (N seconds since 1970.0)
+ 2440900.232j (julian date)
+ 99/02/23,03:22:18 (date string)
+ (separators may be anything except space, +, -)
+ 99:02:15:12:23:30
+ 99:02:15:12h23m30s
+ */
+
+
+/* fseek with timeout - 0.5 sec */
+int Fseek (FILE *f, long offset, int whence) {
+
+ int status, k;
+
+ status = fseek (f, offset, whence);
+ if (status == -1) {
+ int k;
+ /* fprintf (stderr, "problem seeking position: %d\n", errno); */
+ for (k = 0; (k < 10) && ((status = fseek (f, 0, SEEK_SET)) == -1); k++) usleep (50000);
+ if (status == -1) {
+ /* fprintf (stderr, "ERROR: serious problem seeking position: %d\n", errno); */
+ /* clear lock at this point? */
+ return (0);
+ }
+ }
+ return (1);
+}
Index: /branches/elixir/Ohana/src/relphot/src/wcatalog.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/wcatalog.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/wcatalog.c (revision 2485)
@@ -0,0 +1,25 @@
+# include "relphot.h"
+
+void wcatalog (Catalog *catalog) {
+
+ fits_modify (&catalog[0].header, "NRPHOT", "%t", 1, TRUE);
+
+ save_catalog (catalog, VERBOSE);
+
+ if (catalog[0].Naverage) {
+ free (catalog[0].average);
+ catalog[0].Naverage = 0;
+ }
+ if (catalog[0].Nmeasure) {
+ free (catalog[0].measure);
+ catalog[0].Nmeasure = 0;
+ }
+ if (catalog[0].Nmissing) {
+ free (catalog[0].missing);
+ catalog[0].Nmissing = 0;
+ }
+ if (catalog[0].Nsecfilt) {
+ free (catalog[0].secfilt);
+ catalog[0].Nsecfilt = 0;
+ }
+}
Index: /branches/elixir/Ohana/src/relphot/src/write_coords.c
===================================================================
--- /branches/elixir/Ohana/src/relphot/src/write_coords.c (revision 2485)
+++ /branches/elixir/Ohana/src/relphot/src/write_coords.c (revision 2485)
@@ -0,0 +1,29 @@
+# include "relphot.h"
+
+# define CD_COORDS 1
+
+void write_coords (Header *header, Coords *coords) {
+
+ fits_modify (header, "CTYPE1", "%s", 1, "RA---TAN");
+ fits_modify (header, "CTYPE2", "%s", 1, "DEC--TAN");
+
+ fits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1);
+ fits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2);
+
+ fits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1);
+ fits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2);
+
+# if (CD_COORDS)
+ fits_modify (header, "CD1_1", "%le", 1, coords[0].pc1_1 * coords[0].cdelt1);
+ fits_modify (header, "CD2_1", "%le", 1, coords[0].pc2_1 * coords[0].cdelt1);
+ fits_modify (header, "CD1_2", "%le", 1, coords[0].pc1_2 * coords[0].cdelt2);
+ fits_modify (header, "CD2_2", "%le", 1, coords[0].pc2_2 * coords[0].cdelt2);
+# else
+ fits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1);
+ fits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2);
+ fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
+ fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
+ fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
+ fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
+# endif
+}