Index: /tags/relphot-1-5-1/Ohana/src/relphot/Makefile =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/Makefile (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/Makefile (revision 11064) @@ -0,0 +1,45 @@ +default: relphot +help: + @echo "make options: relphot (default)" + +include ../../Configure +HOME = $(ROOT)/src/relphot +BIN = $(HOME)/bin +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +include ../../Makefile.Common + +LDFLAGS := -lkapa -ldvo -lFITS -lohana $(LDFLAGS) + +relphot: $(BIN)/relphot.$(ARCH) +install: $(DESTBIN)/relphot + +RELPHOT = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/GridOps.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/global_stats.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/plot_scatter.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/reload_catalogs.$(ARCH).o \ +$(SRC)/relphot.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ +$(SRC)/setExclusions.$(ARCH).o \ +$(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o \ +$(SRC)/write_coords.$(ARCH).o + +$(RELPHOT): $(INC)/relphot.h +$(BIN)/relphot.$(ARCH): $(RELPHOT) Index: /tags/relphot-1-5-1/Ohana/src/relphot/bin/.cvsignore =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/bin/.cvsignore (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/bin/.cvsignore (revision 11064) @@ -0,0 +1,2 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix +*.linrh Index: /tags/relphot-1-5-1/Ohana/src/relphot/doc/ChangeLog.txt =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/doc/ChangeLog.txt (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/doc/ChangeLog.txt (revision 11064) @@ -0,0 +1,22 @@ + +- relphot-1-4: + * converted to gfits APIs (forces libfits 1.6) + * converted to new DVO APIs (forces libdvo 1.3) + * allow seletion by RA/DEC limits + * removed old, unused code + * added time range to time selections + +- relphot-1-3: + * dropped IMAGE_CATALOG from config + * convert to use libkapa + * dropped _PS from average.R,D,M,dM + * convert to dvo_image_lock,unlock + +- relphot-1-2: released 2005.10.20 + * substantial changes to sync with addstar v1.3 + (added mode/format support, new dvo load functions, etc). + * various cleanups (functions moved to libohana) + * changes to use new libohana (v1.5) / libfits (v1.4) + +- relphot-1-1 + * minor cleanup, loneos.h -> dvo.h Index: /tags/relphot-1-5-1/Ohana/src/relphot/doc/notes.txt =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/doc/notes.txt (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/doc/notes.txt (revision 11064) @@ -0,0 +1,135 @@ + +2006.10.17 + + I am working on a partial upgrade of relphot to handle more general + problems than I have treated in the past. In particular, I would + like to be able to have it assign the average magnitudes, regardless + of whether any or all images and/or objects are well-treated for + relative photometry. In fact, I would like it to be able to + determine average magnitudes even for sources which have only + externally supplied data, and thus have no matching images. I would + eventually like to have more flexibility about the filtering which + is performed on the measurements to determine the average + magnitudes. Some of this work will require a bit of a more careful + treatment of the process. + + At the moment, we can consider the relphot process to consist of the + following steps: + + 1) load the complete catalog data + 2) select a subset of the sources: + - appropriate photcode (equiv == selected photcode) + - in time range, if specified + - appropriate dophot type, if specified + - bright enough: + - mag < MAG_LIM + - dMag > SIGMA_LIM + - ImagMin < iMag < ImagMax + - other restrictions + - within valid image region for chip + 3) iterate to a solution for Mcal (image offsets) + - mark specific images as bad + - mark specific sources as bad + - mark specific measurements as bad + - include desired restrictions on Mcal (eg, all chips of mosaic + matched) + - only keep the image offsets + + 4) if the results are to be kept, recalculate the average magnitudes + using the image exclusions and Mcal values determined above + +2006.05.03 + + I am upgrading relphot to accept as input an arbitrary region (using + the standard ohana SkyRegion concept of a RA & DEC bounded patch on + the sky). This modification has a few implications. First, the + analysis is in fact performed on the populated SkyRegions which + overlap the requested area on the sky. + +- 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: /tags/relphot-1-5-1/Ohana/src/relphot/include/relphot.h =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/include/relphot.h (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/include/relphot.h (revision 11064) @@ -0,0 +1,213 @@ +# include +# include +# include +# include + +/* # define GRID_V1 */ +# define GRID_V2 +# define NO_IMAGE -100 + +# if (0) +typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, color; + double lweight, size; +} Graphdata; +# endif + +typedef struct { + unsigned int start; + unsigned int stop; + float Mcal; + float dMcal; + short Xm; + float 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 CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char CameraConfig[256]; +char SKY_TABLE[256]; +int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +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; + +int DophotSelect, DophotValue; + +double PlotMmin, PlotMmax, PlotdMmin, PlotdMmax; +enum {black, white, red, orange, yellow, green, blue, indigo, violet}; + +int TimeSelect; +time_t TSTART, TSTOP; + +SkyRegion UserPatch; +int UserPatchSelect; + +# 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 + +/*** relphot prototypes ***/ +void ConfigInit PROTO((int *argc, char **argv)); +void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); +char *GetPhotnamebyCode PROTO((PhotCodeData *photcodes, int code)); +void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); +int *SelectRefMosaic PROTO((Mosaic **refmosaic, int *Nimage)); +int args PROTO((int argc, char **argv)); +int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); +void clean_images PROTO(()); +void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final)); +void clean_mosaics PROTO(()); +void clean_stars PROTO((Catalog *catalog, int Ncatalog)); +int corner_check PROTO((double *x1, double *y1, double *x2, double *y2)); +void dumpGrid PROTO(()); +void dump_grid PROTO(()); +int edge_check PROTO((double *x1, double *y1, double *x2, double *y2)); +void findImages PROTO((Catalog *catalog, int Ncatalog)); +int findMosaics PROTO((Catalog *catalog, int Ncatalog)); + +void set_db (FITS_DB *in); +int Shutdown (char *format, ...); +void TrapSignal (int sig); +void SetProtect (int mode); +int SetSignals (); + +void freeGridBins PROTO((int Ncatalog)); +void freeImageBins PROTO((int Ncatalog)); +void freeMosaicBins PROTO((int Ncatalog)); +void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); +int gcatalog PROTO((Catalog *catalog, int FINAL)); +Coords *getCoords PROTO((int meas, int cat)); +float getMcal PROTO((int meas, int cat)); +float getMgrid PROTO((int meas, int cat)); +float getMmos PROTO((int meas, int cat)); +float getMrel PROTO((Catalog *catalog, int meas, int cat)); +Image *getimage PROTO((int N)); +Image *getimages PROTO((int *N)); +void global_stats PROTO((Catalog *catalog, int Ncatalog)); +void initGrid PROTO((int dX, int dY)); +void initGridBins PROTO((Catalog *catalog, int Ncatalog)); +void initImageBins PROTO((Catalog *catalog, int Ncatalog)); +void initImages PROTO((Image *input, int N)); +void initMosaicBins PROTO((Catalog *catalog, int Ncatalog)); +void initMosaicGrid PROTO((Image *image, int Nimage)); +void initMosaics PROTO((Image *image, int Nimage)); +void initMrel PROTO((Catalog *catalog, int Ncatalog)); +void initialize PROTO((int argc, char **argv)); +void initstats PROTO((char *mode)); +int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog)); +SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect)); +Image *select_images PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage)); + +int main PROTO((int argc, char **argv)); +void mark_images PROTO((Image *image, int Nimage, Image *timage, int Ntimage)); +void matchImage PROTO((Catalog *catalog, int meas, int cat)); +void matchMosaics PROTO((Catalog *catalog, int meas, int cat)); +double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); +void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); +void plot_defaults PROTO((Graphdata *graphdata)); +void plot_grid PROTO((Catalog *catalog)); +void plot_images PROTO(()); +void plot_list PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file)); +void plot_mosaic_fields PROTO((Catalog *catalog)); +void plot_mosaics PROTO(()); +void plot_scatter PROTO((Catalog *catalog, int Ncatalog)); +void plot_star_coords PROTO((Catalog *catalog, int Ncatalog)); +void plot_stars PROTO((Catalog *catalog, int Ncatalog)); +void reload_catalogs PROTO((SkyList *skylist)); +int reload_images PROTO((FITS_DB *db)); +int setExclusions PROTO((Catalog *catalog, int Ncatalog)); +void setMcal PROTO((Catalog *catalog, int Poor)); +void setMcalFinal PROTO(()); +int setMcalOutput PROTO((Catalog *catalog, int Ncatalog)); +void setMgrid PROTO((Catalog *catalog)); +int setMmos PROTO((Catalog *catalog, int Poor)); +int setMrel PROTO((Catalog *catalog, int Ncatalog)); +void setMrelFinal PROTO((Catalog *catalog)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark)); +void set_ZP PROTO((double ZERO)); +int setrefcode PROTO((Image *image, int Nimage)); +void skip_measurements PROTO((Catalog *catalog, int pass)); +void sortA PROTO((double *X, int N)); +void sortB PROTO((double *X, double *Y, int N)); +void sortC PROTO((double *X, double *Y, double *F1, double *F2, int N)); +void sortD PROTO((double *X, double *Y, double *Z, int N)); +StatType statsImageM PROTO((Catalog *catalog)); +StatType statsImageN PROTO((Catalog *catalog)); +StatType statsImageX PROTO((Catalog *catalog)); +StatType statsImagedM PROTO((Catalog *catalog)); +StatType statsMosaicM PROTO((Catalog *catalog)); +StatType statsMosaicN PROTO((Catalog *catalog)); +StatType statsMosaicX PROTO((Catalog *catalog)); +StatType statsMosaicdM PROTO((Catalog *catalog)); +StatType statsStarN PROTO((Catalog *catalog, int Ncatalog)); +StatType statsStarS PROTO((Catalog *catalog, int Ncatalog)); +StatType statsStarX PROTO((Catalog *catalog, int Ncatalog)); +void wcatalog PROTO((Catalog *catalog)); +void wimages PROTO(()); +void write_coords PROTO((Header *header, Coords *coords)); Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/ConfigInit.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/ConfigInit.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/ConfigInit.c (revision 11064) @@ -0,0 +1,75 @@ +# include "relphot.h" + +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); + + 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, "GSCFILE", "%s", 0, GSCFILE); + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "CATMODE", "%s", 0, CATMODE); + ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { + SKY_DEPTH = 2; + } + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { + SKY_TABLE[0] = 0; + } + + 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 (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + 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) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/GridOps.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/GridOps.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/GridOps.c (revision 11064) @@ -0,0 +1,585 @@ +# include "relphot.h" + +static int Ngrid; +static float *gridM; +static float *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, float, Ngrid); + ALLOCATE (gridS, float, Ngrid); + ALLOCATE (gridN, int, Ngrid); + bzero (gridM, Ngrid*sizeof(float)); + bzero (gridS, Ngrid*sizeof(float)); + 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, float, Ngrid); + ALLOCATE (gridS, float, Ngrid); + ALLOCATE (gridN, int, Ngrid); + bzero (gridM, Ngrid*sizeof(float)); + bzero (gridS, Ngrid*sizeof(float)); + 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); + + fprintf (stderr, "error: star out of grid\n"); + exit (1); +} +# endif + +float getMgrid (int meas, int cat) { + + int i; + float 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; + float 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; + Mcal = getMcal (m, c); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, c); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue; + + n = catalog[c].measure[m].averef; + Msys = PhotSys (&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_PS, 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, Narea; + float 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); + + Narea = 0; + 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) { + Narea ++; + continue; + } + Mcal = getMcal (m, c); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, c); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mrel = getMrel (catalog, m, c)) == NO_MAG) continue; + + n = catalog[c].measure[m].averef; + Msys = PhotSys (&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]; + N++; + } + } + + fprintf (stderr, "skipped %d meas for area\n", Narea); + + 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); + +} + +void dump_grid () { + + int i, j, Nimage; + int *imlist; + FILE *f; + Header header, theader; + Matrix matrix; + Mosaic *refmosaic; + + /* 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 */ + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + gfits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3); + gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis"); + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_matrix (&matrix); + + /* save grid mag values */ + gfits_init_header (&theader); + theader.Naxes = 2; + theader.Naxis[0] = gridX; + theader.Naxis[1] = gridY; + theader.bitpix = -32; + gfits_create_Theader (&theader, "IMAGE"); + gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET"); + gfits_create_matrix (&theader, &matrix); + for (i = 0; i < gridX; i++) { + for (j = 0; j < gridY; j++) { + gfits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]); + } + } + write_coords (&theader, &refmosaic[0].coords); + gfits_fwrite_header (f, &theader); + gfits_fwrite_matrix (f, &matrix); + gfits_free_matrix (&matrix); + + /* save grid Nmeas values */ + gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS"); + gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_create_matrix (&theader, &matrix); + for (i = 0; i < gridX; i++) { + for (j = 0; j < gridY; j++) { + gfits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]); + } + } + write_coords (&theader, &refmosaic[0].coords); + gfits_fwrite_header (f, &theader); + gfits_fwrite_matrix (f, &matrix); + gfits_free_matrix (&matrix); + + /* save grid sigma values */ + gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA"); + gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_create_matrix (&theader, &matrix); + for (i = 0; i < gridX; i++) { + for (j = 0; j < gridY; j++) { + gfits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]); + } + } + write_coords (&theader, &refmosaic[0].coords); + gfits_fwrite_header (f, &theader); + gfits_fwrite_matrix (f, &matrix); + gfits_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); + + gfits_modify (&theader, "EXTNAME", "%s", 1, ccdname); + gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "PHOTCODE", "%s", 1, pname); + gfits_modify (&theader, "NX", "%d", 1, image[i].NX); + gfits_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; + gfits_modify (&theader, "NAXIS1", "%d", 1, Nx); + gfits_modify (&theader, "NAXIS2", "%d", 1, Ny); + gfits_create_matrix (&theader, &matrix); + + InterpolateGrid ((float *)matrix.buffer, Nx, Ny, &image[0].coords, &refmosaic[0].coords); + gfits_fwrite_header (f, &theader); + gfits_fwrite_matrix (f, &matrix); + gfits_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; + + gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[i]); + gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "NX", "%d", 1, camera.Nx); + gfits_modify (&theader, "NY", "%d", 1, camera.Ny); + + theader.Naxis[0] = RELPHOT_GRID_X; + theader.Naxis[1] = RELPHOT_GRID_Y; + gfits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X); + gfits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y); + gfits_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; + gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]); + } + } + gfits_fwrite_header (f, &theader); + gfits_fwrite_matrix (f, &matrix); + gfits_free_matrix (&matrix); + } +# endif + +} + +void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) { + + int i, j; + double x, y, r, d, X, Y, dx, dy; + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/ImageOps.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/ImageOps.c (revision 11064) @@ -0,0 +1,470 @@ +# 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; + + 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_PS / 3600.0; + dec = catalog[cat].average[ave].D - measure[0].dD_PS / 3600.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); */ +} + +float getMcal (int meas, int cat) { + + int i; + float value; + + i = bin[cat][meas]; + if (i == -1) return (NO_IMAGE); + + if (image[i].code & IMAGE_BAD) return (NO_MAG); + value = image[i].Mcal_PS; + return (value); +} + +Coords *getCoords (int meas, int cat) { + + int i; + + 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; + float 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; + Mmos = getMmos (m, c); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) 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 = PhotSys (&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_PS, 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_PS = stats.mean; + image[i].dMcal_PS = 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_PS); + slist[N] = image[i].dMcal_PS; + 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", MedOffset, stats.median, MaxScatter, MaxOffset); + + Nmark = 0; + for (i = 0; i < Nimage; i++) { + mark = FALSE; + image[i].code &= ~ID_IMAGE_POOR; + mark = (image[i].dMcal_PS > MaxScatter) || (fabs(image[i].Mcal_PS - 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, bin; + 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] = image[i].Mcal_PS; + dlist[i] = image[i].dMcal_PS; + xlist[i] = image[i].secz_PS; + } + + 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.00025*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Nimage; i++) { + bin = image[i].dMcal_PS / 0.00025; + 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, m, c, n, N; + double *list, *dlist; + float Mcal, Mmos; + 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]; + + Mcal = getMcal (m, c); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, c); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) 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, n; + 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, n; + 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] = image[i].Mcal_PS; + dlist[n] = 1; + n++; + } + + liststats (list, dlist, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsImagedM (Catalog *catalog) { + + int i, n; + 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] = image[i].dMcal_PS; + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/MosaicOps.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/MosaicOps.c (revision 11064) @@ -0,0 +1,678 @@ +# 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_PS; + + /* 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 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_PS; + dMcal += image[m].dMcal_PS; + Xm += image[m].Xm; + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ + image[m].Mcal_PS = 0.0; + image[m].dMcal_PS = 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 () { + + int i, j, im, Nimage; + Image *image; + + if (!MOSAICNAME[0]) return; + + image = getimages (&Nimage); + + for (i = 0; i < Nmosaic; i++) { + for (j = 0; j < Nimlist[i]; j++) { + im = imlist[i][j]; + image[im].Mcal_PS = mosaic[i].Mcal; + image[im].dMcal_PS = 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 (int 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); + } + } + return (TRUE); +} + +void matchMosaics (Catalog *catalog, int meas, int cat) { + + int i; + + 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_PS - catalog[cat].measure[meas].dR_PS / 3600.0; + dec = catalog[cat].average[ave].D_PS - catalog[cat].measure[meas].dD_PS / 3600.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; +} + +float getMmos (int meas, int cat) { + + int i; + float value; + + if (!MOSAICNAME[0]) return (0); + i = bin[cat][meas]; + if (i == -1) return (NO_IMAGE); + + 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; + float Msys, Mrel, Mcal, Mgrid; + double *list, *dlist, *Mlist, *dMlist; + StatType stats; + Image *image; + + 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]; + + Mcal = getMcal (m, c); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + if (catalog[c].measure[m].flags & MEAS_BAD) 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 = PhotSys (&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_PS, MIN_ERROR); + Mlist[N] = Msys; + dMlist[N] = 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; + 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] = 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; + 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] = 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, m, c, n, N; + double *list, *dlist; + float Mcal; + 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]; + + Mcal = getMcal (m, c); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) 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; + 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, N, mark, Nmark; + double *mlist, *slist, *dlist; + double MaxOffset, MedOffset, MaxScatter; + StatType stats; + + 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", MedOffset, stats.median, MaxScatter, 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; + 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]; + + if (catalog[c].measure[m].flags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + + ave = catalog[c].measure[m].averef; + xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR_PS / 3600.0; + ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD_PS / 3600.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, bin; + 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] = mosaic[i].Mcal; + dlist[i] = mosaic[i].dMcal; + xlist[i] = 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); + graphdata.size = 1.5; + graphdata.ptype = 7; + plot_list (&graphdata, Mlist, dlist, Nmosaic, "Mcal vs dMcal", "MdM.png"); + +# define NBIN 200 + REALLOCATE (xlist, double, NBIN); + REALLOCATE (Mlist, double, NBIN); + + /**** dMcal histgram ****/ + for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Nmosaic; i++) { + bin = mosaic[i].dMcal / 0.00005; + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/SetSignals.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/SetSignals.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/SetSignals.c (revision 11064) @@ -0,0 +1,118 @@ +# include "relphot.h" + +static int Protect = FALSE; +static int Trapped = FALSE; + +void TrapSignal (int sig) { + fprintf (stderr, "trapped signal %d\n", sig); + if (sig == 11) { + fprintf (stderr, "seg fault\n"); + exit (1); + } + if (Protect) { + Trapped = TRUE; + fprintf (stderr, "blocking until protected sections are clear\n"); + return; + } + Shutdown ("halted by signal (trapped)"); +} + +void SetProtect (int mode) { + Protect = mode; + if (Trapped && !Protect) Shutdown ("halted by signal (protect)"); +} + +int SetSignals () { + + int i; + + /* disable almost all signal interrupts */ + for (i = 0; i < 36; i++) { + switch (i) { + /* can't redirect this signals */ + case SIGKILL: /* kill -9: cannot be caught or ignored */ + case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ + case SIGCHLD: /* child halted: ignore */ + case SIGPWR: /* power failure - why ignore this? */ + case SIGWINCH: /* window resized */ + case SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ + break; + + default: + signal (i, TrapSignal); + } + } + return (TRUE); +} +/* + + Signal Value Action Comment + ------------------------------------------------------------------------- + SIGHUP 1 A Hangup detected on controlling terminal + or death of controlling process + SIGINT 2 A Interrupt from keyboard + SIGQUIT 3 A Quit from keyboard + SIGILL 4 A Illegal Instruction + SIGABRT 6 C Abort signal from abort(3) + SIGFPE 8 C Floating point exception + SIGKILL 9 AEF Kill signal + SIGSEGV 11 C Invalid memory reference + SIGPIPE 13 A Broken pipe: write to pipe with no readers + SIGALRM 14 A Timer signal from alarm(2) + SIGTERM 15 A Termination signal + SIGUSR1 30,10,16 A User-defined signal 1 + SIGUSR2 31,12,17 A User-defined signal 2 + SIGCHLD 20,17,18 B Child stopped or terminated + SIGCONT 19,18,25 Continue if stopped + SIGSTOP 17,19,23 DEF Stop process + SIGTSTP 18,20,24 D Stop typed at tty + SIGTTIN 21,21,26 D tty input for background process + SIGTTOU 22,22,27 D tty output for background process + + Next various other signals. + + Signal Value Action Comment + --------------------------------------------------------------------- + SIGTRAP 5 CG Trace/breakpoint trap + SIGIOT 6 CG IOT trap. A synonym for SIGABRT + SIGEMT 7,-,7 G + SIGBUS 10,7,10 AG Bus error + SIGSYS 12,-,12 G Bad argument to routine (SVID) + SIGSTKFLT -,16,- AG Stack fault on coprocessor + SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) + SIGIO 23,29,22 AG I/O now possible (4.2 BSD) + SIGPOLL AG A synonym for SIGIO (System V) + SIGCLD -,-,18 G A synonym for SIGCHLD + SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) + SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) + SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) + SIGPROF 27,27,29 AG Profile alarm clock + SIGPWR 29,30,19 AG Power failure (System V) + SIGINFO 29,-,- G A synonym for SIGPWR + SIGLOST -,-,- AG File lock lost + SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) + SIGUNUSED -,31,- AG Unused signal + (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and + sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a + sparc.) + + The letters in the "Action" column have the following meanings: + + A Default action is to terminate the process. + + B Default action is to ignore the signal. + + C Default action is to dump core. + + D Default action is to stop the process. + + E Signal cannot be caught. + + F Signal cannot be ignored. + + G Not a POSIX.1 conformant signal. + +*/ Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/Shutdown.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/Shutdown.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/Shutdown.c (revision 11064) @@ -0,0 +1,28 @@ +# include "relphot.h" + +static FITS_DB *db; + +void set_db (FITS_DB *in) { + db = in; +} + +/* clean up open / locked ImageCat before shutting down */ +int Shutdown (char *format, ...) { + va_list argp; + char *formatplus; + + ALLOCATE (formatplus, char, strlen(format)); + strcpy (formatplus, format); + strcat (formatplus, "\n"); + + va_start (argp, format); + vfprintf (stderr, formatplus, argp); + free (formatplus); + va_end (argp); + + SetProtect (TRUE); + gfits_db_close (db); + fprintf (stderr, "ERROR: addstar halted\n"); + exit (1); +} + Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/StarOps.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/StarOps.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/StarOps.c (revision 11064) @@ -0,0 +1,565 @@ +# 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)); +} + +float getMrel (Catalog *catalog, int meas, int cat) { + + int ave; + float 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_PS; + return (value); +} + +int setMrel (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, N; + float 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; + + // XXX allow REF stars to be included in the calculation + // this should be optionally set, and should allow for + // REF stars to be downweighted by more than their reported + // errors. how such info is carried is unclear... + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) { + Mcal = Mmos = Mgrid = 0; + } else { + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + } + + Msys = PhotSys (&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_PS, 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_PS = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = stats.sigma; + catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; + } + } + } + return (TRUE); +} + +int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) { + + int i, j, k, m, N, Nmax; + float 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; + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) { + Mcal = Mmos = Mgrid = 0; + } else { + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + } + + Msys = PhotSys (&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_PS, 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_PS = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = 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); + return (TRUE); +} + +/* 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; + float 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; + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + catalog[i].measure[m].Mcal_PS = Mcal + Mmos + Mgrid; + } + } + } + return (TRUE); +} + +void clean_stars (Catalog *catalog, int Ncatalog) { + + int i, j, Ndel, Nave, Ntot, mark; + float 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_PS; + 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", 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_PS; + 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; + float 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; */ + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&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_PS, 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; */ + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&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_PS, 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); + if (VERBOSE) 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; + float Mcal, Mmos; + 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++) { + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) 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, 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, n, Ntot; + double *list, *dlist; + float dM; + 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_PS; + list[n] = 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, bin; + float 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.00025*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_PS; + bin = dMrel / 0.00025; + 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 ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS; + 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, N; + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/args.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/args.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/args.c (revision 11064) @@ -0,0 +1,209 @@ +# include "relphot.h" +void usage (void); + +int args (int argc, char **argv) { + + int N; + double trange; + + /* 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_dtime (argv[N], &trange)) { + if (!str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + UserPatchSelect = FALSE; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatchSelect = 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 = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = 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 = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; + } + + if ( UserPatchSelect && (argc != 2)) usage (); + if (!UserPatchSelect && (argc != 3)) usage (); + + return TRUE; +} + +void usage () { + fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); + fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\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"); + exit (2); +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/bcatalog.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/bcatalog.c (revision 11064) @@ -0,0 +1,116 @@ +# include "relphot.h" + +int bcatalog (Catalog *subcatalog, Catalog *catalog) { + + int i, j, offset, ecode; + int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; + float 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; + + /* 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) { + float *p; + p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M_PS; + *p = NO_MAG; + p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM_PS; + *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 */ + // XXX ignore this criterion for REF measurements? + if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue; + + /* select measurements by mag limit */ + mag = PhotCat (&catalog[0].measure[offset]); + if (mag > MAG_LIM) continue; + + /* select measurements by measurement error */ + if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM_PS > SIGMA_LIM)) continue; + + /* select measurements by mag limit */ + if (ImagSelect) { + mag = PhotInst (&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_PS = 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); + } + } + + // XXXX test : what checks do I need to make elsewhere to avoid problems here? + # if 0 + if (Nm < 2) { /* enough measurements in band? */ + Nmeasure -= Nm; + continue; + } + # endif + 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 (VERBOSE) { + fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, + subcatalog[0].Naverage, subcatalog[0].Nmeasure); + } + return (TRUE); +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/global_stats.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/global_stats.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/global_stats.c (revision 11064) @@ -0,0 +1,44 @@ +# include "relphot.h" + +void 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: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas); + fprintf (stderr, "meas / mosaic: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas); + fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); + + fprintf (stderr, "chisq image: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas); + fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); + fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); + + fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); + fprintf (stderr, "dMcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas); + + fprintf (stderr, "Mcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas); + fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); + + fprintf (stderr, "dMrel star: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas); + + initstats (STATMODE); + +} + + Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/initialize.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/initialize.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/initialize.c (revision 11064) @@ -0,0 +1,54 @@ +# include "relphot.h" + +void initialize (int argc, char **argv) { + + int N; + + ConfigInit (&argc, argv); + args (argc, argv); + + N = UserPatchSelect ? 1 : 2; + if ((photcode = GetPhotcodebyName (argv[N])) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); + 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[N]); + 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, "GSCFILE %s\n", GSCFILE); + fprintf (stderr, "CATDIR %s\n", CATDIR); + exit (0); + } +} + Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/liststats.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/liststats.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/liststats.c (revision 11064) @@ -0,0 +1,119 @@ +# 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; + +void 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; + + ke = ks = dMo = 0; + + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/load_catalogs.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/load_catalogs.c (revision 11064) @@ -0,0 +1,52 @@ +# include "relphot.h" + +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) { + + int i, Nstar; + Catalog *catalog, tcatalog; + + if (VERBOSE) fprintf (stderr, "loading catalog data\n"); + + ALLOCATE (catalog, Catalog, skylist[0].Nregions); + + // load data from each region file, only use bright stars + for (i = 0; i < skylist[0].Nregions; i++) { + + // set up the basic catalog info + tcatalog.filename = skylist[0].filename[i]; + tcatalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + tcatalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + tcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + tcatalog.Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog + + if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename); + exit (1); + } + if (VERBOSE && !tcatalog.Nave_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); + + // select only the brighter stars + bcatalog (&catalog[i], &tcatalog); + dvo_catalog_unlock (&tcatalog); + dvo_catalog_free (&tcatalog); + } + + // XXX keep this test? + Nstar = 0; + for (i = 0; i < skylist[0].Nregions; i++) { + Nstar += catalog[i].Naverage; + } + if (Nstar < 2) { + fprintf (stderr, "warning: insufficient stars %d\n", Nstar); + } + + // XXX consider only returning the populated catalogs + *Ncatalog = skylist[0].Nregions; + return (catalog); +} + +/* this function loads all relevant catalog files for the first pass. it currently loads the data + read only (SOFT lock) since it assumes the image table has been locked. if we go to the new + addstar locking paradigm, in which the images and catalogs are updated independently, then we may + need to use an XCLD lock here. +*/ Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/load_images.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/load_images.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/load_images.c (revision 11064) @@ -0,0 +1,54 @@ +# include "relphot.h" + +SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) { + + Image *image, *subset; + int Nimage, Nsubset, Nchar; + int *LineNumber; + + SkyTable *sky = NULL; + SkyList *skylist = NULL; + + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area + if (RegionSelect) { + skylist = SkyListByPatch (sky, -1, region); + } else { + Nchar = strlen(regionName); + if (!strcmp (®ionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0; + skylist = SkyListByName (sky, regionName); + } + + // convert database table to internal structure + image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + + // select the images which overlap the selected sky regions + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + + gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); + + initImages (subset, Nsubset); + initMosaics (subset, Nsubset); + + return (skylist); +} + +int reload_images (FITS_DB *db) { + + Image *image; + int Nimage, Nx, i; + VTable *vtable; + + image = getimages (&Nimage); + + vtable = &db[0].vtable; + + gfits_scan (vtable[0].header, "NAXIS1", "%d", 1, &Nx); + for (i = 0; i < Nimage; i++) { + memcpy (vtable[0].buffer[i], &image[i], Nx); + } + return (TRUE); +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/plot_scatter.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/plot_scatter.c (revision 11064) @@ -0,0 +1,61 @@ +# include "relphot.h" + +void plot_scatter (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, N, Ntot; + float Mrel, Mcal, Mmos, Mgrid; + double *xlist, *ylist, *ilist; + 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; + Mcal = getMcal (m, i); + if (Mcal == NO_MAG) continue; + if (Mcal == NO_IMAGE) continue; + Mmos = getMmos (m, i); + if (Mmos == NO_MAG) continue; + if (Mmos == NO_IMAGE) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS; + xlist[N] = Mrel; + ylist[N] = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel; + ilist[N] = PhotInst (&catalog[i].measure[m]); + 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", "iMag.png"); + free (xlist); + free (ylist); + free (ilist); +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/plotstuff.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/plotstuff.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/plotstuff.c (revision 11064) @@ -0,0 +1,155 @@ +# include "relphot.h" +# include + +static int Xgraph[5] = {0,0,0,0,0}; +static int active; + +/* +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) + +void XDead () { + signal (SIGPIPE, XDead); + fprintf (stderr, "kapa is dead, must restart\n"); + Xgraph[active] = -1; +} + +int open_graph (int N) { + + Xgraph[N] = KiiOpen ("kapa", NULL); + if (Xgraph[N] < 0) { + fprintf (stderr, "error starting kapa\n"); + return (FALSE); + } + return (TRUE); +} + +void DonePlotting (Graphdata *graphmode, int N) { + + if (Xgraph[N] == 0) return; + + KapaBox (Xgraph[N], graphmode); + return; +} + +void JpegPlot (Graphdata *graphmode, int N, char *filename) { + + if (Xgraph[N] == 0) return; + + KapaPNG (Xgraph[N], filename); + return; +} + +void PSPlot (Graphdata *graphmode, int N, char *filename) { + + if (Xgraph[N] == 0) return; + + KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename, "default"); + return; +} + +void PrepPlotting (int Npts, Graphdata *graphmode, int N) { + + if (Npts < 1) return; + + if (Xgraph[N] < 1) if (!open_graph(N)) return; + + KapaClear (Xgraph[N], TRUE); + + KiiSendCommand (Xgraph[N], 4, "PLOT"); + + /* send Xgraph the plot details */ + KiiSendMessage (Xgraph[N], "%8d %8d %d %d %d %d %d %f %f", + Npts, graphmode[0].style, + graphmode[0].ptype, graphmode[0].ltype, + graphmode[0].etype, graphmode[0].ebar, graphmode[0].color, + graphmode[0].lweight, graphmode[0].size); + KiiSendMessage (Xgraph[N], "%g %g %g %g", + graphmode[0].xmin, graphmode[0].xmax, + graphmode[0].ymin, graphmode[0].ymax); + + return; +} + +void PlotLabel (char *string, int N) { + + if (Xgraph[N] == 0) return; + + KapaSendLabel (Xgraph[N], string, 2); +} + +void PlotVector (int Npts, double *vect, int mode, int N) { + + float *values; + int i, Nbytes; + + 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); + return; +} + +void plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file) { + + int i; + 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); + } +} + +void plot_defaults (Graphdata *graphdata) { + + KapaInitGraph (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: /tags/relphot-1-5-1/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/reload_catalogs.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/reload_catalogs.c (revision 11064) @@ -0,0 +1,55 @@ +# include "relphot.h" + +void reload_catalogs (SkyList *skylist) { + + int i; + int status; + struct stat filestat; + Catalog catalog; + + if (VERBOSE) fprintf (stderr, "re-loading catalog data\n"); + + /* load data from each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + catalog.filename = skylist[0].filename[i]; + + // only update existing db tables + status = stat (catalog.filename, &filestat); + if ((status == -1) && (errno == ENOENT)) { + if (VERBOSE) fprintf (stderr, "no file %s, skipping\n", catalog.filename); + continue; + } + + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (VERBOSE && (catalog.Nave_disk == 0)) { + fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + initImageBins (&catalog, 1); + initMosaicBins (&catalog, 1); + initGridBins (&catalog, 1); + + findImages (&catalog, 1); + findMosaics (&catalog, 1); + + setMrelFinal (&catalog); + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + freeImageBins (1); + freeMosaicBins (1); + freeGridBins (1); + } +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/relphot.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/relphot.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/relphot.c (revision 11064) @@ -0,0 +1,118 @@ +# include "relphot.h" + +int main (int argc, char **argv) { + + int i, status, Ncatalog; + Catalog *catalog; + FITS_DB db; + + SkyList *skylist = NULL; + + /* get configuration info, args */ + initialize (argc, argv); + + /* register database handle with shutdown procedure */ + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + // XXX get ZERO_POINT from config + dvo_image_create (&db, 25.0); + // Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + } else { + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + } + + /* load regions and images based on specified sky patch */ + // XXX need to mimic old-style load by passing patch name + // XXX need to reduce number of global variables in use. + // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true + skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect); + + /* unlock, if we can (else, unlocked below) */ + if (!UPDATE) dvo_image_unlock (&db); + + /* load catalog data from region files */ + catalog = load_catalogs (skylist, &Ncatalog); + + /* 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 */ + for (i = 0; i < Ncatalog; i++) dvo_catalog_free (&catalog[i]); + freeImageBins (Ncatalog); + freeMosaicBins (Ncatalog); + freeGridBins (Ncatalog); + + /* load catalog data from region files, update Mrel include all data */ + reload_catalogs (skylist); + setMcalFinal (); + + reload_images (&db); + + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + exit (0); +} + Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/select_images.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/select_images.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/select_images.c (revision 11064) @@ -0,0 +1,236 @@ +# include "relphot.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). 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 */ + +typedef struct { + double Xc[5]; + double Yc[5]; + double Rc; + double Dc; +} SkyRegionCoords; + +Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) { + + Image *image; + int i, j, k, m, found, nimage, NIMAGE; + int InRange, ecode; + double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + int *line_number; + Coords tcoords; + SkyRegionCoords *skycoords; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + // the comparison is made in the catalog local projection. below we set crval1,2 + 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"); + + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + skycoords[i].Yc[4] = skycoords[i].Yc[0]; + + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); + dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + } + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + nimage = 0; + NIMAGE = 100; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, int, NIMAGE); + + // go through the complete list of images, selecting ones which overlap any region + for (i = 0; i < Ntimage; i++) { + + /* exclude images by photcode */ + ecode = GetPhotcodeEquivCodebyCode (timage[i].source); + if (ecode != photcode[0].code) continue; + + /* exclude images by time */ + if (TimeSelect) { + if (timage[i].tzero < TSTART) continue; + if (timage[i].tzero > TSTOP) 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 ra,dec */ + for (j = 0; j < 5; j++) { + XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + } + + /* compare with each region file */ + for (m = 0; (m < skylist[0].Nregions) && !found; m++) { + + /* we make positional comparisons in the projection of catalog */ + tcoords.crval1 = skycoords[m].Rc; + tcoords.crval2 = skycoords[m].Dc; + + /* transform to ra,dec */ + InRange = TRUE; + for (j = 0; (j < 5) && InRange; j++) { + InRange = RD_to_XY (&Xi[j], &Yi[j], Ri[j], Di[j], &tcoords); + } + if (!InRange) continue; + + /* check if image corner inside catalog */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]); + } + /* check if catalog corner inside image */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&skycoords[m].Xc[j], &skycoords[m].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], &skycoords[m].Xc[k], &skycoords[m].Yc[k]); + } + } + if (!found) continue; + + image[nimage] = timage[i]; + /* always allow 'few' images 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_PS = 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); + } + } + } + + if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, int, MAX (nimage, 1)); + free (skycoords); + + *Nimage = nimage; + *LineNumber = line_number; + 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/setExclusions.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/setExclusions.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/setExclusions.c (revision 11064) @@ -0,0 +1,53 @@ +# include "relphot.h" + +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_PS / 3600.0; + d = catalog[i].average[j].D + catalog[i].measure[m].dD_PS / 3600.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); + return (TRUE); +} Index: /tags/relphot-1-5-1/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/setMrelFinal.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/setMrelFinal.c (revision 11064) @@ -0,0 +1,155 @@ +# include "relphot.h" + +void setMrelFinal (Catalog *catalog) { + + int i, j, m, ecode; + float *p; + short *q; + + /* 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_PS; + *p = NO_MAG; + p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM_PS; + *p = NO_MAG; + q = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm; + *q = 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_PS = 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, k, m, ecode, d1, d2; + int Ntot, Ntry, Nkeep, Nskip; + float mag; + + 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 */ + /* do NOT skip measurements without a matching image */ + 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 = PhotInst (&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 ++; + } + } + if (VERBOSE) 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/sort.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/sort.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/sort.c (revision 11064) @@ -0,0 +1,182 @@ +# include "relphot.h" + +void 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; + } +} + +void 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; + } +} + +void 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; + } +} + +void 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: /tags/relphot-1-5-1/Ohana/src/relphot/src/write_coords.c =================================================================== --- /tags/relphot-1-5-1/Ohana/src/relphot/src/write_coords.c (revision 11064) +++ /tags/relphot-1-5-1/Ohana/src/relphot/src/write_coords.c (revision 11064) @@ -0,0 +1,29 @@ +# include "relphot.h" + +# define CD_COORDS 1 + +void write_coords (Header *header, Coords *coords) { + + gfits_modify (header, "CTYPE1", "%s", 1, "RA---TAN"); + gfits_modify (header, "CTYPE2", "%s", 1, "DEC--TAN"); + + gfits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1); + gfits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2); + + gfits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1); + gfits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2); + +# if (CD_COORDS) + gfits_modify (header, "CD1_1", "%le", 1, coords[0].pc1_1 * coords[0].cdelt1); + gfits_modify (header, "CD2_1", "%le", 1, coords[0].pc2_1 * coords[0].cdelt1); + gfits_modify (header, "CD1_2", "%le", 1, coords[0].pc1_2 * coords[0].cdelt2); + gfits_modify (header, "CD2_2", "%le", 1, coords[0].pc2_2 * coords[0].cdelt2); +# else + gfits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1); + gfits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2); + gfits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1); + gfits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2); + gfits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1); + gfits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2); +# endif +}