Index: /tags/libdvo-1-0/Ohana/src/libdvo/Makefile =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/Makefile (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/Makefile (revision 5388) @@ -0,0 +1,113 @@ +include ../../Configure + +HOME = $(ROOT)/src/libdvo +AUTO = $(ROOT)/src/libautocode +ASRC = $(AUTO)/src +ADEF = $(AUTO)/def +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +CFLAGS = -I$(LINC) -D$(ARCH) + +default: install + +install: $(DESTLIB)/libdvo.a + +libdvo: $(LIB)/libdvo.$(ARCH).a + +INCS = $(DESTINC)/dvo.h $(DESTINC)/autocode.h + +OBJS = \ +$(SRC)/coordops.$(ARCH).o \ +$(SRC)/LoadPhotcodes.$(ARCH).o \ +$(SRC)/imreg_datatypes.$(ARCH).o \ +$(SRC)/mosaic_astrom.$(ARCH).o \ +$(SRC)/fits_db.$(ARCH).o \ +$(SRC)/photfits.$(ARCH).o \ +$(SRC)/dvo_image.$(ARCH).o \ +$(SRC)/dvo_image_raw.$(ARCH).o \ +$(SRC)/dvo_catalog.$(ARCH).o \ +$(SRC)/dvo_catalog_raw.$(ARCH).o \ +$(SRC)/dvo_catalog_mef.$(ARCH).o \ +$(SRC)/dvo_catalog_split.$(ARCH).o \ +$(SRC)/dvo_convert.$(ARCH).o \ +$(SRC)/dvo_convert_elixir.$(ARCH).o \ +$(SRC)/dvo_convert_loneos.$(ARCH).o \ +$(SRC)/dvo_convert_panstarrs.$(ARCH).o + +AOBJ = \ +$(ASRC)/coords.$(ARCH).o \ +$(ASRC)/average.$(ARCH).o \ +$(ASRC)/average-loneos.$(ARCH).o \ +$(ASRC)/average-elixir.$(ARCH).o \ +$(ASRC)/average-panstarrs.$(ARCH).o \ +$(ASRC)/secfilt.$(ARCH).o \ +$(ASRC)/secfilt-loneos.$(ARCH).o \ +$(ASRC)/secfilt-elixir.$(ARCH).o \ +$(ASRC)/secfilt-panstarrs.$(ARCH).o \ +$(ASRC)/measure.$(ARCH).o \ +$(ASRC)/measure-loneos.$(ARCH).o \ +$(ASRC)/measure-elixir.$(ARCH).o \ +$(ASRC)/measure-panstarrs.$(ARCH).o \ +$(ASRC)/missing.$(ARCH).o \ +$(ASRC)/photcode.$(ARCH).o \ +$(ASRC)/image.$(ARCH).o \ +$(ASRC)/image-loneos.$(ARCH).o \ +$(ASRC)/image-elixir.$(ARCH).o \ +$(ASRC)/image-panstarrs.$(ARCH).o \ +$(ASRC)/regimage.$(ARCH).o \ +$(ASRC)/detreg.$(ARCH).o \ +$(ASRC)/photreg.$(ARCH).o \ +$(ASRC)/photreg-old.$(ARCH).o \ +$(ASRC)/smpdata.$(ARCH).o \ +$(ASRC)/spectrum.$(ARCH).o \ +$(ASRC)/spectrum-ascii.$(ARCH).o \ +$(ASRC)/Stars.$(ARCH).o \ +$(ASRC)/GSCRegion.$(ARCH).o \ +$(ASRC)/AddstarClientOptions.$(ARCH).o + +ACOM = $(ADEF)/autocode.c $(ADEF)/autocode.h $(ADEF)/common.h + +$(OBJS): $(INCS) $(AOBJ) + +$(AOBJ): $(ASRC)/%.$(ARCH).o : $(ADEF)/%.d $(ACOM) + cd $(AUTO) && make install + +$(DESTINC)/autocode.h: + cd $(AUTO) && make install + +$(DESTLIB)/libdvo.a: $(LIB)/libdvo.$(ARCH).a + @if [ ! -d $(DESTLIB) ]; then mkdir -p $(DESTLIB); fi + rm -f $(DESTLIB)/libdvo.a + cp $(LIB)/libdvo.$(ARCH).a $(DESTLIB)/libdvo.a + +$(LIB)/libdvo.$(ARCH).a: $(AOBJ) $(OBJS) + @if [ ! -d $(LIB) ]; then mkdir -p $(LIB); fi + rm -f $(LIB)/libdvo.$(ARCH).a + ar rcv $(LIB)/libdvo.$(ARCH).a $(AOBJ) $(OBJS) + $(RANLIB) $(LIB)/libdvo.$(ARCH).a + +$(DESTINC)/%: $(INC)/% + @if [ ! -d $(DESTINC) ]; then mkdir -p $(DESTINC); fi + rm -f $(DESTINC)/$* + cp $(INC)/$* $(DESTINC)/ + +clean: + rm -f */*~ + rm -f */#* + rm -f */*.o + rm -f */*.a + cd $(AUTO) && make clean + +dist: clean + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $*.c -o $@ Index: /tags/libdvo-1-0/Ohana/src/libdvo/include/dvo.h =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/include/dvo.h (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/include/dvo.h (revision 5388) @@ -0,0 +1,324 @@ +# include +# include +# include + +# ifndef DVO_H +# define DVO_H + +# define ELIXIR 1 +# define PANSTARRS 0 +# define LONEOS 0 + +/*** named data values (convert all to enums?) ***/ + +/* DVO table modes */ +enum {DVO_MODE_UNDEF, DVO_MODE_RAW, DVO_MODE_MEF, DVO_MODE_SPLIT, DVO_MODE_MYSQL} DVOTableMode; + +/* DVO table formats */ +enum {DVO_FORMAT_UNDEF, DVO_FORMAT_INTERNAL, DVO_FORMAT_ELIXIR, DVO_FORMAT_LONEOS, DVO_FORMAT_PANSTARRS} DVOTableFormat; + +/* image data modes in RegImage */ +enum {T_UNDEF = -1, T_NONE, T_OBJECT, T_DARK, T_BIAS, T_FLAT, T_MASK, T_FRINGE, T_SCATTER, T_MODES, T_FRINGEPTS, T_ANY, N_TYPE}; +enum {M_UNDEF = -1, M_NONE, M_MEF, M_SPLIT, M_SINGLE, M_CUBE, M_SLICE, M_MODES, N_MODE}; + +/* RegImage.flag values */ +# define IMREG_DIST 0x01 /* image distributed, only imregister-3.0 */ + +/* catalog values to be loaded */ +# define LOAD_AVES 0x01 +# define LOAD_MEAS 0x02 +# define LOAD_MISS 0x04 +# define LOAD_SECF 0x08 +# define LOAD_MEAS_META 0x10 + +/* invalid mag value */ +# define NO_MAG 0x7fff +# define NO_ERR 0xff + +/* photometry code types */ +# define PHOT_PRI 0x01 +# define PHOT_SEC 0x02 +# define PHOT_DEP 0x03 +# define PHOT_REF 0x04 +# define PHOT_ALT 0x05 /* never stored, only for look-ups */ + +/* Image.code values. these are codes to note bad images */ +# define ID_IMAGE_NEW 0x0000 /* no nrphot attempted */ +# define ID_IMAGE_NOCAL 0x0001 /* used within nrphot to mean "don't apply fit" */ +# define ID_IMAGE_POOR 0x0002 /* relphot says image is bad */ +# define ID_IMAGE_SKIP 0x0004 /* external information image is bad */ +# define ID_IMAGE_FEW 0x0008 /* currently too few measurements for good value */ + +/* Measure.flags values */ +# define ID_MEAS_NOCAL 0x0001 +# define ID_MEAS_POOR 0x0002 +# define ID_MEAS_SKIP 0x0004 +# define ID_MEAS_AREA 0x0008 +# define BLEND_IMAGE 0x0100 +# define BLEND_CATALOG 0x0200 +# define BLEND_IMAGE_NEIGHBOR 0x1000 +# define ID_MEAS_TRAIL 0x2000 +# define ID_MEAS_GHOST 0x4000 + +/* some subtle distinctions between the blend flags: + BLEND_IMAGE: the star on an image is matched with more + than one star in the catalog (image has worse seeing than catalog) + BLEND_CATALOG: the star in the catalog is matched with more + than one star on the image (image has better seeing than catalog) + CALIBRATED: relative photometry has been performed on this measurement + BLEND_IMAGE_NEIGHBOR: the star on an image is matched with more + than one star in the catalog, but not in the same catalog file. +*/ + +/* Average.code values */ +# define ID_STAR_FEW 0x0001 /* used within relphot: skip star */ +# define ID_STAR_POOR 0x0002 /* used within relphot: skip star */ +# define ID_PROPER 0x0400 /* star with large proper motion */ +# define ID_TRANSIENT 0x1000 /* is this mutually exclusive with USNO? */ +# define ID_VARIABLE 0x2000 /* not currently set? */ +# define ID_ASTEROID 0x2000 /* identified with an asteroid */ +# define ID_BAD_OBJECT 0x4000 /* if all measurements are bad, set this bit */ +# define ID_MOVING 0x8000 +# define ID_ROCK 0xa000 /* 0x8000 + 0x2000 */ +# define ID_GHOST 0xc001 /* 0x8000 + 0x4000 + 0x0001 */ +# define ID_TRAIL 0xc002 /* 0x8000 + 0x4000 + 0x0002 */ +# define ID_BLEED 0xc003 /* 0x8000 + 0x4000 + 0x0003 */ +# define ID_COSMIC 0xc004 /* 0x8000 + 0x4000 + 0x0004 */ + +/*** general dvo structures (internal use only / not IO) ***/ + +/* FITS DB structure */ +typedef struct { + FILE *f; + char *filename; + int dbstate; + int lockstate; + double timeout; + Header header; + Matrix matrix; + Header theader; + FTable ftable; + VTable vtable; + int mode; /* what data storage mode is used for disk file? */ + int format; /* what data format is used for disk file? */ + int virtual; /* is table in ftable or vtable? */ + int swapped; /* is table in internal byte-order? */ +} FITS_DB; + +# if (0) +/* structure for data on a catalog region */ +typedef struct { + char filename[256]; + double DEC[2], RA[2]; +} GSCRegion; +# endif + +typedef struct { + int Ncode; + int Nsecfilt; + int hashcode[0x10000]; + int hashNsec[0x10000]; + PhotCode *code; +} PhotCodeData; + +/* a catalog contains this data */ +typedef struct Catalog { + char *filename; /* catalog file */ + FILE *f; /* file descriptor */ + Header header; + + /* data in the catalog file */ + Average *average; + Measure *measure; + Missing *missing; + SecFilt *secfilt; + int Naverage, Nmeasure, Nmissing, Nsecfilt; /* current number of each component */ + int Nave_disk, Nmeas_disk, Nmiss_disk; /* number of component on disk */ + int Nmeas_off; /* dist seq of first loaded data value */ + /* note the different counting for Nsecfilt */ + + /* pointers to split data files */ + struct Catalog *measure_catalog; /* measure catalog data (split) */ + struct Catalog *missing_catalog; /* missing catalog data (split) */ + struct Catalog *secfilt_catalog; /* secfilt catalog data (split) */ + + /* extra catalog information */ + int lockmode; + int catmode; /* storage mode (raw, mef, split, mysql) */ + int catformat; /* storage format (elixir, panstarrs, etc) */ + int catflags; /* choices to be loaded */ + int sorted; /* is measure table average-sorted? */ + + /* pointers for data manipulation */ + int *found; + int *image; + int *mosaic; + float *X; + float *Y; + +} Catalog; + +/*** prototypes ***/ + +/* in fits_db.c */ +int fits_db_init PROTO((FITS_DB *db)); +int fits_db_create PROTO((FITS_DB *db)); +int fits_db_lock PROTO((FITS_DB *db, char *filename)); +int fits_db_load PROTO((FITS_DB *db)); +int fits_db_load_segment PROTO((FITS_DB *db, int start, int Nrows)); +int fits_db_save PROTO((FITS_DB *db)); +int fits_db_update PROTO((FITS_DB *db)); +int fits_db_close PROTO((FITS_DB *db)); +int fits_db_free PROTO((FITS_DB *db)); + +/* in coords.c, using libautocode/def/coords.d */ +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords); +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords); +int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords); +int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords); +int GetCoords (Coords *coords, Header *header); +int PutCoords (Coords *coords, Header *header); +void RegisterMosaic (Coords *coords); +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); + +int lock_catalog (Catalog *catalog, int lockmode); +int unlock_catalog (Catalog *catalog); +int load_catalog (Catalog *catalog, int VERBOSE); +int save_catalog (Catalog *catalog, char VERBOSE); +int update_catalog (Catalog *catalog, char VERBOSE); + +int FindMosaicForImage (Image *images, int Nimages, int entry); +int FindMosaicForImage_TableSearch (Image *images, int Nimages, int entry); +int FindMosaicForImage_MatchSearch (Image *images, int Nimages, int entry); +int BuildChipMatch (Image *images, int Nimages); + +short int putMi (double value); +double getMi (short int value); +void returnMcal (Image *image, double *c); +void assignMcal (Image *image, double *c, int order); +double applyMcal (Image *image, double x, double y); +double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order); + +PhotCode *GetPhotcodebyName (char *name); +PhotCode *GetPhotcodeEquivbyName (char *name); +PhotCode *GetPhotcodebyCode (int code); +PhotCode *GetPhotcodebyNsec (int Nsec); +PhotCode *GetPhotcodeEquivbyCode (int code); + +char *GetPhotcodeNamebyCode (int code); + +float PhotInst (Measure *measure); +float PhotCat (Measure *measure); +float PhotSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt); + +# if (0) +short iPhotInst (Measure *measure); +short iPhotAbs (Measure *measure); +short iPhotCat (Measure *measure); +short iPhotSys (Measure *measure, Average *average, SecFilt *secfilt); +short iPhotRel (Measure *measure, Average *average, SecFilt *secfilt); +short iPhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +short iPhotAve (PhotCode *code, Average *average, SecFilt *secfilt); +short iPhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +short iPhotXm (PhotCode *code, Average *average, SecFilt *secfilt); +short iPhotdM (PhotCode *code, Average *average, SecFilt *secfilt); +# endif + +float iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color); + +int LoadPhotcodes (char *filename); +int GetPhotcodeCodebyName (char *name); +int GetPhotcodeEquivCodebyName (char *name); +int GetPhotcodeEquivCodebyCode (int code); +int GetPhotcodeNsec (int code); +int GetPhotcodeNsecfilt (); +void SetZeroPoint (double ZP); +int *GetPhotcodeEquivList (int code, int *nlist); +void ParseColorTerms (char *terms, float *X, int *N); + +int get_image_type (char *name); +char *get_type_name (int type); +int get_image_mode (char *name); +char *get_mode_name (int mode); + +/* these functions refer to the DVO structures defined in the includes above + * they are being replaced with autocode entries (drop when totally autocoded) + */ + +int Fread (void *ptr, int size, int nitems, FILE *f, char *type); +int Fwrite (void *ptr, int size, int nitems, FILE *f, char *type); +int ByteSwap (char *ptr, int size, int nitems, char *type); +int ConvertStruct (char *buffer, int size, int Nbytes, char *type); + +/*** conversion functions / I/O conversions ***/ + +Average *ReadRawAverage (FILE *f, int Naverage, int format); +Measure *ReadRawMeasure (FILE *f, int Nmeasure, int format); +SecFilt *ReadRawSecFilt (FILE *f, int Nsecfilt, int format); +int WriteRawAverage (FILE *f, Average *average, int Naverage, int format); +int WriteRawMeasure (FILE *f, Measure *measure, int Nmeasure, int format); +int WriteRawSecFilt (FILE *f, SecFilt *secfilt, int Nsecfilt, int format); + +Average *FtableToAverage (FTable *ftable, int *Naverage, int *format); +Average *AverageLoneosToInternal (AverageLoneos *in, int Nvalues); +Average *AverageElixirToInternal (AverageElixir *in, int Nvalues); +Average *AveragePanstarrsToInternal (AveragePanstarrs *in, int Nvalues); +AverageLoneos *AverageInternalToLoneos (Average *in, int Nvalues); +AverageElixir *AverageInternalToElixir (Average *in, int Nvalues); +AveragePanstarrs *AverageInternalToPanstarrs (Average *in, int Nvalues); + +Measure *FtableToMeasure (FTable *ftable, int *Nmeasure, int *format); +Measure *MeasureLoneosToInternal (MeasureLoneos *in, int Nvalues); +Measure *MeasureElixirToInternal (MeasureElixir *in, int Nvalues); +Measure *MeasurePanstarrsToInternal (MeasurePanstarrs *in, int Nvalues); +MeasureLoneos *MeasureInternalToLoneos (Measure *in, int Nvalues); +MeasureElixir *MeasureInternalToElixir (Measure *in, int Nvalues); +MeasurePanstarrs *MeasureInternalToPanstarrs (Measure *in, int Nvalues); + +SecFilt *FtableToSecFilt (FTable *ftable, int *Nsecfilt, int *format); +SecFilt *SecFiltLoneosToInternal (SecFiltLoneos *in, int Nvalues); +SecFilt *SecFiltElixirToInternal (SecFiltElixir *in, int Nvalues); +SecFilt *SecFiltPanstarrsToInternal (SecFiltPanstarrs *in, int Nvalues); +SecFiltLoneos *SecFiltInternalToLoneos (SecFilt *in, int Nvalues); +SecFiltElixir *SecFiltInternalToElixir (SecFilt *in, int Nvalues); +SecFiltPanstarrs *SecFiltInternalToPanstarrs (SecFilt *in, int Nvalues); + +int AverageToFtable (FTable *ftable, Average *average, int Naverage, int format); +int MeasureToFtable (FTable *ftable, Measure *measure, int Nmeasure, int format); +int SecFiltToFtable (FTable *ftable, SecFilt *secfilt, int Nsecfilt, int format); + +int load_catalog_raw (Catalog *catalog, int VERBOSE); +int save_catalog_raw (Catalog *catalog, char VERBOSE); +int load_catalog_mef (Catalog *catalog, int VERBOSE); +int save_catalog_mef (Catalog *catalog, char VERBOSE); +int load_catalog_split (Catalog *catalog, int VERBOSE); +int save_catalog_split (Catalog *catalog, char VERBOSE); +int update_catalog_split (Catalog *catalog, char VERBOSE); + +/*** DVO image db I/O Functions ***/ +int dvo_image_load (FITS_DB *db, int VERBOSE, int FORCE_READ); +int dvo_image_save (FITS_DB *db, int VERBOSE); +int dvo_image_update (FITS_DB *db, int VERBOSE); +int dvo_image_load_raw (FITS_DB *db, int VERBOSE, int FORCE_READ); +int dvo_image_update_raw (FITS_DB *db, int VERBOSE); +int dvo_image_save_raw (FITS_DB *db, int VERBOSE); +int dvo_image_addrows (FITS_DB *db, Image *new, int Nnew); + +int FtableToImage (FTable *ftable, Header *theader, int *format); +int ImageToFtable (FTable *ftable, Header *theader, int format); +int ImageToVtable (VTable *vtable, Header *theader, int format); +Image *ImageElixirToInternal (ImageElixir *in, int Nvalues); +ImageElixir *ImageInternalToElixir (Image *in, int Nvalues); +Image *ImageLoneosToInternal (ImageLoneos *in, int Nvalues); +ImageLoneos *ImageInternalToLoneos (Image *in, int Nvalues); +Image *ImagePanstarrsToInternal (ImagePanstarrs *in, int Nvalues); +ImagePanstarrs *ImageInternalToPanstarrs (Image *in, int Nvalues); + +# endif Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/LoadPhotcodes.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/LoadPhotcodes.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/LoadPhotcodes.c (revision 5388) @@ -0,0 +1,798 @@ +# include + +# define NCTERMS 4 +# define F_PS 0.001 +# define NO_MAG_PS 100.0 + +static PhotCodeData *photcodes = NULL; +static double ZERO_POINT; +static short int Nseclist[0x10000]; + +/* static short int iZERO_POINT; */ + +void SetZeroPoint (double ZP) { + ZERO_POINT = ZP; + /* iZERO_POINT = 1000 * ZP; */ +} + +int LoadPhotcodes (char *filename) { + + FILE *f; + int i, Ns, Np, NPHOTCODE, Npri, Nfield; + PhotCode *photcode; + int code; + char *c; + char line[256]; + char name[32], type[32], Zero[32], Airmass[32], Offset[32], + C1[32], C2[32], Slope[32], Color[32], Primary[32]; + int c1, c2; + + /* allocate space to photcode table, free existing data */ + if (photcodes == NULL) { + ALLOCATE (photcodes, PhotCodeData, 1); + photcodes[0].code = NULL; + } + if (photcodes[0].code != NULL) free (photcodes[0].code); + + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + photcodes[0].Ncode = 0; + photcodes[0].Nsecfilt = 0; + photcodes[0].code = (PhotCode *) NULL; + return (FALSE); + } + + Np = 0; + NPHOTCODE = 10; + ALLOCATE (photcode, PhotCode, NPHOTCODE); + photcodes[0].Nsecfilt = 0; + for (i = 0; i < 0x10000; i++) { + photcodes[0].hashcode[i] = -1; + photcodes[0].hashNsec[i] = -1; + } + + while (scan_line (f, line) != EOF) { + for (c = line; isspace (*c); c++); + if (*c == '#') continue; + Nfield = sscanf (c, "%d %s %s %s %s %s %s %s %s %s %s", + &code, name, type, Zero, Airmass, Offset, C1, C2, Slope, Color, Primary); + if (Nfield != 11) { continue; } + + c1 = atof (C1); + c2 = atof (C2); + if (!strcmp (C1, "-")) { c1 = 0; } + if (!strcmp (C2, "-")) { c2 = 0; } + + photcode[Np].type = 0; + photcode[Np].code = code; + strcpy (photcode[Np].name, name); + if (!strncasecmp (type, "pri", 3)) { + photcode[Np].type = PHOT_PRI; + photcode[Np].C = 1000*atof (Zero); + photcode[Np].K = atof (Airmass); + photcode[Np].dC = 1000*atof (Offset); + photcode[Np].dX = 1000*atof (Color); + photcode[Np].c1 = c1; + photcode[Np].c2 = c2; + photcode[Np].equiv = atoi (Primary); + ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc); + Nseclist[0] = Np; + } + if (!strncasecmp (type, "sec", 3)) { + photcode[Np].type = PHOT_SEC; + photcode[Np].C = 1000*atof (Zero); + photcode[Np].K = atof (Airmass); + photcode[Np].dC = 1000*atof (Offset); + photcode[Np].dX = 1000*atof (Color); + photcode[Np].c1 = c1; + photcode[Np].c2 = c2; + photcode[Np].equiv = atoi (Primary); + photcodes[0].Nsecfilt ++; + ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc); + Nseclist[photcodes[0].Nsecfilt] = Np; + } + if (!strncasecmp (type, "dep", 3)) { + photcode[Np].type = PHOT_DEP; + photcode[Np].C = 1000*atof (Zero); /* zero point in millimags */ + photcode[Np].K = atof (Airmass); /* airmass coeff (millimag / millimag) */ + photcode[Np].dC = 1000*atof (Offset); /* color ref z.p. (millimag) */ + photcode[Np].dX = 1000*atof (Color); /* average color (millimag) */ + photcode[Np].c1 = c1; + photcode[Np].c2 = c2; + photcode[Np].equiv = atoi (Primary); + ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc); + } + if (!strncasecmp (type, "ref", 3)) { + photcode[Np].type = PHOT_REF; + photcode[Np].C = 0; + photcode[Np].K = 0; + photcode[Np].dC = 0; + photcode[Np].dX = 0; + photcode[Np].c1 = 0; + photcode[Np].c2 = 0; + photcode[Np].equiv = 0; + photcode[Np].X[0] = 0; + photcode[Np].Nc = 0; + } + + /* alt photcodes are a little different: they have the SAME photcode as an existing + pri/sec photcode, but define an alternate transformation for that code */ + if (!strncasecmp (type, "alt", 3)) { + photcode[Np].type = PHOT_ALT; + photcode[Np].C = 1000*atof (Zero); /* zero point in millimags */ + photcode[Np].K = atof (Airmass); /* airmass coeff (millimag / millimag) */ + photcode[Np].dC = 1000*atof (Offset); /* color ref z.p. (millimag) */ + photcode[Np].dX = 1000*atof (Color); /* average color (millimag) */ + photcode[Np].c1 = c1; + photcode[Np].c2 = c2; + photcode[Np].equiv = atoi (Primary); + ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc); + } + if (!photcode[Np].type) { + fprintf (stderr, "error in Photfile: unknown type %s\n", type); + } + + Np++; + if (Np == NPHOTCODE) { + NPHOTCODE += 10; + REALLOCATE (photcode, PhotCode, NPHOTCODE); + } + } + fclose (f); + + /* set up hashcode for photcode refs: + * the hashcode gives the structure sequence for a given photcode: + * photcode[i].hashcode[photcode[i].code] == i + */ + + Ns = 0; + for (i = 0; i < Np; i++) { + if (photcode[i].type == PHOT_ALT) continue; /* no hashcode for ALT codes */ + if (photcodes[0].hashcode[photcode[i].code] != -1) { + fprintf (stderr, "duplicate photcodes in file\n"); + free (photcode); + return (FALSE); + } + photcodes[0].hashcode[photcode[i].code] = i; + if (photcode[i].type == PHOT_SEC) { + photcodes[0].hashNsec[photcode[i].code] = Ns; + Ns ++; + } + } + /* validity check for references */ + for (i = 0; i < Np; i++) { + if (photcode[i].type == PHOT_DEP) { + Npri = photcodes[0].hashcode[photcode[i].equiv]; + if ((Npri >= Np) || (Npri < 0)) { + fprintf (stderr, "reference for dependent photcode is not in photcodes\n"); + free (photcode); + return (FALSE); + } + if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) { + fprintf (stderr, "reference for dependent photcode is not a primary or secondary code\n"); + free (photcode); + return (FALSE); + } + } + if (photcode[i].type == PHOT_ALT) { + Npri = photcodes[0].hashcode[photcode[i].code]; + if ((Npri >= Np) || (Npri < 0)) { + fprintf (stderr, "reference for alternate photcode is not in photcodes\n"); + free (photcode); + return (FALSE); + } + if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) { + fprintf (stderr, "reference for alternate photcode is not a primary or secondary code\n"); + free (photcode); + return (FALSE); + } + } + } + photcodes[0].code = photcode; + photcodes[0].Ncode = Np; + + return (TRUE); + +} + +void ParseColorTerms (char *terms, float *X, int *N) { + + int i; + char *p; + + p = terms; + + for (i = 0; (p != NULL) && (i < NCTERMS); i++) { + X[i] = atof (p); + p = strchr (p, ','); + if (p == (char *) NULL) continue; + p ++; + } + *N = i; +} + +/********** photcode lookups **********/ + +/* return photcode for given name */ +PhotCode *GetPhotcodebyName (char *name) { + + int i; + + if (name == (char *) NULL ) return (NULL); + + for (i = 0; i < photcodes[0].Ncode; i++) { + if (!strcmp (photcodes[0].code[i].name, name)) { + return (&photcodes[0].code[i]); + } + } + return (NULL); +} +/* return photcode.code for given name */ +int GetPhotcodeCodebyName (char *name) { + + int i; + + if (name == (char *) NULL ) return (0); + + for (i = 0; i < photcodes[0].Ncode; i++) { + if (!strcmp (photcodes[0].code[i].name, name)) { + return (photcodes[0].code[i].code); + } + } + return (0); +} +/* return equivalent photcode for given name */ +PhotCode *GetPhotcodeEquivbyName (char *name) { + + int i, equiv; + + if (name == (char *) NULL ) return (NULL); + + for (i = 0; i < photcodes[0].Ncode; i++) { + if (!strcmp (photcodes[0].code[i].name, name)) { + if (photcodes[0].code[i].equiv == 0) return (NULL); + equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv]; + if (equiv == -1) return (NULL); + return (&photcodes[0].code[equiv]); + } + } + return (NULL); +} +/* return equivalent photcode.code for given name */ +int GetPhotcodeEquivCodebyName (char *name) { + + int i, equiv; + + if (name == (char *) NULL ) return (0); + + for (i = 0; i < photcodes[0].Ncode; i++) { + if (!strcmp (photcodes[0].code[i].name, name)) { + if (photcodes[0].code[i].equiv == 0) return (0); + equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv]; + if (equiv == -1) return (0); + return (photcodes[0].code[equiv].code); + } + } + return (0); +} + +/* return photcode for given code */ +PhotCode *GetPhotcodebyCode (int code) { + + int entry; + + if (code < 0) return (NULL); + if (code > 0x10000) return (NULL); + + entry = photcodes[0].hashcode[code]; + if (entry == -1) return (NULL); + + return (&photcodes[0].code[entry]); +} +/* return photcode.code for given code */ +char *GetPhotcodeNamebyCode (int code) { + + int entry; + + if (code < 0) return (NULL); + if (code > 0x10000) return (NULL); + + entry = photcodes[0].hashcode[code]; + if (entry == -1) return (NULL); + + return (photcodes[0].code[entry].name); +} +/* return equivalent photcode for given code */ +PhotCode *GetPhotcodeEquivbyCode (int code) { + + int entry, equiv; + + if (code < 0) return (NULL); + if (code > 0x10000) return (NULL); + + entry = photcodes[0].hashcode[code]; + if (entry == -1) return (NULL); + + if (photcodes[0].code[entry].equiv == 0) return (NULL); + equiv = photcodes[0].hashcode[photcodes[0].code[entry].equiv]; + + if (equiv == -1) return (NULL); + return (&photcodes[0].code[equiv]); +} +/* return equivalent photcode.code for given code */ +int GetPhotcodeEquivCodebyCode (int code) { + + int entry; + + if (code < 0) return (0); + if (code > 0x10000) return (0); + + entry = photcodes[0].hashcode[code]; + if (entry == -1) return (0); + return (photcodes[0].code[entry].equiv); +} + +int GetPhotcodeNsec (int code) { + + int Nsec; + + if (code < 0) return (-1); + if (code > 0x10000) return (-1); + + Nsec = photcodes[0].hashNsec[code]; + return (Nsec); +} + +/* Nsec of 0 is PRI */ +PhotCode *GetPhotcodebyNsec (int Nsec) { + + if (Nsec > photcodes[0].Nsecfilt) return (NULL); + if (Nsec < 0) return (NULL); + + return (&photcodes[0].code[Nseclist[Nsec]]); +} + +int GetPhotcodeNsecfilt () { + return (photcodes[0].Nsecfilt); +} + +/* ALLOCATE and return list of all photcodes + with photcode.equiv == code */ +int *GetPhotcodeEquivList (int code, int *nlist) { + + int i, Nlist; + int *list; + + ALLOCATE (list, int, MAX (1, photcodes[0].Ncode)); + Nlist = 0; + for (i = 0; i < photcodes[0].Ncode; i++) { + if (photcodes[0].code[i].equiv != code) continue; + list[Nlist] = photcodes[0].code[i].code; + Nlist ++; + } + REALLOCATE (list, int, MAX (1, Nlist)); + + *nlist = Nlist; + return (list); +} + +/**** conversion to INTERNAL vs TABLE types makes the iPhot versions irrelevant ****/ +# if (0) +/******** photometry conversions *********/ +double PhotInst (Measure *measure) { + + short Mi; + double M; + + Mi = iPhotInst (measure); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotCat (Measure *measure) { + + short Mi; + double M; + + Mi = iPhotCat (measure); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotSys (Measure *measure, Average *average, SecFilt *secfilt) { + + short Mi; + double M; + + Mi = iPhotSys (measure, average, secfilt); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotRel (Measure *measure, Average *average, SecFilt *secfilt) { + + short Mi; + double M; + + Mi = iPhotRel (measure, average, secfilt); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + short Mi; + double M; + + Mi = iPhotCal (thisone, average, secfilt, measure, code); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + short Mi; + double M; + + Mi = iPhotRef (code, average, secfilt, measure); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + short Mi; + double M; + + Mi = iPhotAve (code, average, secfilt); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) { + + short Mi; + double M; + + Mi = iPhotdM (code, average, secfilt); + M = 0.001*Mi; + return (M); +} + +/****/ +double PhotXm (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Mi; + double Xm; + + Mi = iPhotXm (code, average, secfilt); + Xm = (Mi == NO_MAG) ? -1.0 : pow (10.0, 0.01*Mi); + return (Xm); +} + +# endif + +/******** internal photometry conversions (keeps values in short int millimags) *********/ +float PhotInst (Measure *measure) { + + int Np; + float M; + + Np = photcodes[0].hashcode[measure[0].source]; + if (Np == -1) return (NO_MAG_PS); + + if (photcodes[0].code[Np].type == PHOT_REF) { + M = measure[0].M_PS; + return (M); + } + + M = measure[0].M_PS - measure[0].dt_PS - ZERO_POINT; + + return (M); + +} + +float PhotCat (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].source]; + if (Np == -1) return (NO_MAG_PS); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].M_PS; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C; + + return (Mcat); +} + +float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].source]; + if (Np == -1) return (NO_MAG_PS); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Msys = measure[0].M_PS; + return (Msys); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = iPhotColor (average, secfilt, NULL, code); + if (mc == NO_MAG_PS) return (Mcat); + mc = mc - F_PS*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotRel (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].source]; + if (Np == -1) return (NO_MAG_PS); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].M_PS; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].M_PS - ZERO_POINT + code[0].K*(measure[0].airmass_PS - 1.000) + F_PS*code[0].C - measure[0].Mcal_PS; + + /* for DEP, color must be made of PRI/SEC */ + mc = iPhotColor (average, secfilt, NULL, code); + if (mc == NO_MAG_PS) return (Mrel); + mc = mc - F_PS*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].source]; + if (Np == -1) { + return (NO_MAG_PS); + } + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].M_PS; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) { + return (NO_MAG_PS); + } + + Mcal = PhotRel (thisone, average, secfilt) + F_PS*code[0].C; + + mc = iPhotColor (average, secfilt, measure, code); + if (mc == NO_MAG_PS) return (Mcal); + mc = mc - F_PS*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +/* color term may not use DEP magnitude */ +float iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Ns1, Ns2, Ns; + float m1, m2, mc; + PhotCode *color; + + m1 = m2 = NO_MAG_PS; + + if (measure == NULL) { + Ns1 = photcodes[0].hashNsec[code[0].c1]; + Ns2 = photcodes[0].hashNsec[code[0].c2]; + + m1 = (Ns1 == -1) ? average[0].M_PS : secfilt[Ns1].M_PS; + m2 = (Ns2 == -1) ? average[0].M_PS : secfilt[Ns2].M_PS; + mc = ((m1 == NO_MAG_PS) || (m2 == NO_MAG_PS)) ? NO_MAG_PS : (m1 - m2); + return (mc); + } + + /* find magnitude matching first color term */ + color = GetPhotcodebyCode (code[0].c1); + if (color == NULL) return (NO_MAG_PS); + if (color[0].type == PHOT_REF) { + for (i = 0; (i < average[0].Nm) && (m1 == NO_MAG_PS); i++) { + if (measure[i].source == color[0].code) { + m1 = measure[i].M_PS; + } + } + } else { + Ns = photcodes[0].hashNsec[color[0].code]; + m1 = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + } + + /* find magnitude matching second color term */ + color = GetPhotcodebyCode (code[0].c2); + if (color == NULL) return (NO_MAG_PS); + if (color[0].type == PHOT_REF) { + for (i = 0; (i < average[0].Nm) && (m2 == NO_MAG_PS); i++) { + if (measure[i].source == color[0].code) { + m2 = measure[i].M_PS; + } + } + } else { + Ns = photcodes[0].hashNsec[color[0].code]; + m2 = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + } + mc = ((m1 == NO_MAG_PS) || (m2 == NO_MAG_PS)) ? NO_MAG_PS : (m1 - m2); + return (mc); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + Mref = Mave + F_PS*code[0].C; + + mc = iPhotColor (average, secfilt, measure, code); + if (mc == NO_MAG_PS) return (Mref); + mc = mc - F_PS*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +/***/ +float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mave; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + return (Mave); +} + +/*** note that this is NOT a wrapper around iPhotdM ***/ +float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dM; + + Ns = photcodes[0].hashNsec[code[0].code]; + dM = (Ns == -1) ? average[0].dM_PS : secfilt[Ns].dM_PS; + return (dM); +} + +/*** note that this is NOT a wrapper around iPhotXm ***/ +float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + short Mi; + float Xm; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mi = (Ns == -1) ? average[0].Xm : secfilt[Ns].Xm; + Xm = (Mi == NO_MAG_PS) ? -1.0 : pow (10.0, 0.01*Mi); + return (Xm); +} + +/* given a photcode pair c1 & c2, return the color of this star (NaN if not found) */ +int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color) { + + int i, Ns; + double M1, M2, dM; + PhotCode *code; + + code = GetPhotcodebyCode (c1); + if (code == NULL) return (FALSE); + if (code[0].type == PHOT_REF) { + for (i = 0; i < average[0].Nm; i++) { + if (measure[i].source == c1) { + M1 = measure[i].M_PS; + goto filter1; + } + } + return (FALSE); + } else { + Ns = photcodes[0].hashNsec[code[0].code]; + M1 = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + } + + filter1: + code = GetPhotcodebyCode (c2); + if (code == NULL) return (FALSE); + if (code[0].type == PHOT_REF) { + for (i = 0; i < average[0].Nm; i++) { + if (measure[i].source == c2) { + M2 = measure[i].M_PS; + goto filter2; + } + } + return (FALSE); + } else { + Ns = photcodes[0].hashNsec[code[0].code]; + M2 = (Ns == -1) ? average[0].M_PS : secfilt[Ns].M_PS; + } + + filter2: + + dM = M1 - M2; + *color = dM; + + return (TRUE); +} + +/* photcode table should have the following format: + +# code name type zero airmass offset c1 c2 slope primary + 1 B pri 24.0 0.15 - - - - - - + 2 B pri 24.0 0.15 - - - - - - + 3 B1 sec 22.5 0.18 0.15 1 2 0.10 0.50 1 + 1000 USNO_B ref - - - - - - - - + +*/ + + +/* + Nc1 = photcodes[0].code[Np].c1; + Ns1 = photcodes[0].hashNsec[Nc1]; + + Nc2 = photcodes[0].code[Np].c2; + Ns2 = photcodes[0].hashNsec[Nc2]; + + Xlam = photcodes[0].code[Np].X[0]; + Klam = photcodes[0].code[Np].K; + + m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M; + m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M; +*/ Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/coordops.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/coordops.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/coordops.c (revision 5388) @@ -0,0 +1,563 @@ +# include + +static Coords *mosaic = NULL; + +void RegisterMosaic (Coords *coords) { + mosaic = coords; +} + +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) { + + int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl; + char *type; + double L, M, X, Y, T, Z, Z2; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, sdp, cdp; + + *ra = 0; + *dec = 0; + stht = ctht = 1; + type = &coords[0].ctype[4]; + + /* DIS is equiv to TAN with higher order terms + PLY is cartesian with higher order terms + WRP is equiv to PLY, with implied mosaic */ + + Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP"); + Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); + PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR"); + Zenith1 = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG"); + Zenith2 = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM"); + Zenithal = !strcmp(type, "-ZEA") || Zenith1 || Zenith2; + if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE); + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + /** extra polynomial terms **/ + if (Polynomial) { + if (coords[0].Npolyterms > 1) { + L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; + M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; + } + if (coords[0].Npolyterms > 2) { + L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; + M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; + } + } + + /**** Locally Cartesian Projections ****/ + if (Cartesian) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + + /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */ + if (!strcmp(type, "-WRP")) { + if (mosaic == NULL) return (FALSE); + XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic); + } + return (TRUE); + } + + /**** Zenithal Projections ****/ + if (Zenithal) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } else { + sphi = L / R; + cphi = -M / R; + } + + /* this is wrong : STG is not TAN - need to put in correct relationships. but is a close approx */ + if (Zenith1) { + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + } + if (Zenith2) { + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + } + if (!strcmp(type, "-ZEA")) { + stht = 1 - 0.5*SQ(R*RAD_DEG); + ctht = sqrt (1 - stht*stht); + } + + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + return (TRUE); + } + + /**** Other Conventional Projections ****/ + if (PseudoCyl) { + if (!strcmp(type, "-AIT")) { + Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M)); + if (Z2 < 0) return (FALSE); + Z = sqrt (Z2); + alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0); + delta = DEG_RAD * asin (RAD_DEG*M*Z); + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } + if (!strcmp(type, "-GLS")) { + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / cos (RAD_DEG * M); + delta = M; + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } + if (!strcmp(type, "-PAR")) { + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / (1.0 - SQ(2.0*M/180)); + delta = 3 * DEG_RAD * asin (M/180.0); + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } + return (TRUE); + } + return (FALSE); +} + +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) { + + char *type; + int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl; + double phi, theta; + double determ; + double X, Y, L, M, Lo, Mo, dL, dM; + double sphi, cphi, stht; + double salp, calp, sdel, cdel, sdp, cdp; + double P, A, Rc; + + status = TRUE; + *x = 0; + *y = 0; + type = &coords[0].ctype[4]; + L = M = 0; + + Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP"); + Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE"); + PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR"); + Zenith1 = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG"); + Zenith2 = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM"); + Zenithal = !strcmp(type, "-ZEA") || Zenith1 || Zenith2; + if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE); + + /**** Locally Cartesian Projections ****/ + if (Cartesian) { + if (!strcmp(type, "-WRP")) { + if (mosaic == NULL) return (FALSE); + RD_to_XY (&Lo, &Mo, ra, dec, mosaic); + L = (Lo - coords[0].crval1); + M = (Mo - coords[0].crval2); + } else { + L = (ra - coords[0].crval1); + M = (dec - coords[0].crval2); + } + } + + /**** Zenithal Projections ****/ + if (Zenithal) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; /* sin(theta) */ + sphi = cdel*salp; /* = cos(theta)*sin(phi) */ + cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ + if (stht < 0) status = FALSE; + + if (Zenith1) { + L = DEG_RAD * sphi / stht; + M = -DEG_RAD * cphi / stht; + } + if (Zenith2) { + L = DEG_RAD * sphi; + M = -DEG_RAD * cphi; + } + if (!strcmp(type, "-ZEA")) { + Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht); + L = Rc * sphi; + M = -Rc * cphi; + status = TRUE; + } + } + + /**** Other Standard Projections ****/ + if (PseudoCyl) { + if (!strcmp(type, "-AIT")) { + phi = RAD_DEG*(ra - coords[0].crval1); + theta = RAD_DEG*(dec - coords[0].crval2); + P = 1.0 + cos (theta) * cos (0.5*phi); + if (P != 0.0) { + A = DEG_RAD * sqrt (2.0 / P); + L = 2.0 * A * cos (theta) * sin (0.5*phi); + M = A * sin (theta); + } else { + L = 0.0; + M = 0.0; + } + } + if (!strcmp(type, "-GLS")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + L = phi * cos(RAD_DEG * theta); + M = theta; + } + if (!strcmp(type, "-PAR")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1); + M = 180.0 * sin (RAD_DEG*theta/3.0); + } + } + + /* convert L,M to X,Y */ + determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M); + Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L); + + /** extra polynomial terms **/ + if (Polynomial) { + for (i = 0; i < 3; i++) { + Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + if (coords[0].Npolyterms > 1) { + Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; + Mo += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; + } + if (coords[0].Npolyterms > 2) { + Lo += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; + Mo += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; + } + dL = (L - Lo); + dM = (M - Mo); + + X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM); + Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dL); + } + } + /* check for correct size (iterate?) */ + + *x = X / coords[0].cdelt1 + coords[0].crpix1; + *y = Y / coords[0].cdelt2 + coords[0].crpix2; + + return (status); + +} + +int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords) { + + int status; + double tmpx, tmpy; + + status = RD_to_XY (&tmpx, &tmpy, ra, dec, coords); + *x = tmpx; + *y = tmpy; + + return (status); + +} + +int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords) { + + int status; + double tmpr, tmpd; + + status = XY_to_RD (&tmpr, &tmpd, x, y, coords); + *ra = tmpr; + *dec = tmpd; + + return (status); + +} + +int GetCoords (Coords *coords, Header *header) { + + int status, itmp, Polynomial; + double Lambda, rotate, scale; + double equinox; + char *ctype; + + rotate = 0.0; + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + coords[0].Npolyterms = 1; + strcpy (coords[0].ctype, "NONE"); + + status = FALSE; + if (fits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) { + status = fits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= fits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); + status &= fits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); + + if (fits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) { + status &= fits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2); + if (fits_scan (header, "CROTA2", "%f", 1, &rotate)) { + Lambda = coords[0].cdelt2 / coords[0].cdelt1; + coords[0].pc1_1 = cos(rotate*RAD_DEG); + coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda; + coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda; + coords[0].pc2_2 = cos(rotate*RAD_DEG); + } + if (fits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2); + status &= fits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1); + status &= fits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2); + } + + ctype = &coords[0].ctype[4]; + Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP"); + + if (Polynomial) { + if (fits_scan (header, "NPLYTERM", "%d", 1, &itmp)) { + coords[0].Npolyterms = itmp; + } else { + coords[0].Npolyterms = 3; + } + switch (coords[0].Npolyterms) { + case 3: + status &= fits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]); + status &= fits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]); + status &= fits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]); + status &= fits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]); + status &= fits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]); + status &= fits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]); + status &= fits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]); + status &= fits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]); + case 2: + status &= fits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]); + status &= fits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]); + status &= fits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]); + status &= fits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]); + status &= fits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]); + status &= fits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]); + case 0: + case 1: + break; + } + } + } else { + if (fits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2); + status &= fits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1); + status &= fits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2); + /* renormalize */ + scale = hypot (coords[0].pc1_1, coords[0].pc1_2); + coords[0].cdelt1 = coords[0].cdelt2 = scale; + coords[0].pc1_1 /= scale; + coords[0].pc1_2 /= scale; + coords[0].pc2_1 /= scale; + coords[0].pc2_2 /= scale; + } else { + status = FALSE; + } + } + } else { + /* some of my thesis data uses this simple linear model - convert on read? */ + if (fits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) { + status = fits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1); + status &= fits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2); + status &= fits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "DEC_X", "%f", 1, &coords[0].pc2_1); + status &= fits_scan (header, "DEC_Y", "%f", 1, &coords[0].pc2_2); + coords[0].crpix1 = coords[0].crpix2 = 0.0; + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + strcpy (coords[0].ctype, "GENE"); + } + } + if (status) { + if (!fits_scan (header, "EQUINOX", "%lf", 1, &equinox)) { + if (!fits_scan (header, "EPOCH", "%lf", 1, &equinox)) { + equinox = 2000.0; + } + } + if (fabs (equinox - 2000.0) > 0.1) { + coords_precess (&coords[0].crval1, &coords[0].crval2, equinox, 2000.0); + } + } + if (!status) { + fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype); + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + } + return (status); +} + +int PutCoords (Coords *coords, Header *header) { + + int OldAIPS; + char csys[16], ctype[16]; + double rotate, Lambda; + + /* modifications to the ctype? */ + OldAIPS = FALSE; + fits_modify (header, "CTYPE2", "%s", 1, coords[0].ctype); + if (!strcmp(coords[0].ctype, "MM")) { + fits_modify (header, "CTYPE1", "%s", 1, "LL"); + OldAIPS = TRUE; + } else { + strcpy (csys, "NONE"); + if (!strncmp (coords[0].ctype, "DEC-", 4)) strcpy (csys, "RA--"); + if (!strncmp (coords[0].ctype, "GLAT", 4)) strcpy (csys, "GLON"); + if (!strncmp (coords[0].ctype, "ELAT", 4)) strcpy (csys, "ELON"); + if (!strncmp (coords[0].ctype, "HLAT", 4)) strcpy (csys, "HLON"); + if (!strncmp (coords[0].ctype, "SLAT", 4)) strcpy (csys, "SLON"); + if (!strcmp (csys, "NONE")) return (FALSE); + sprintf (ctype, "%s-%s", csys, &coords[0].ctype[5]); + fits_modify (header, "CTYPE1", "%s", 1, ctype); + } + + fits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1); + fits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2); + fits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1); + fits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2); + fits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1); + fits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2); + + if (OldAIPS) { + Lambda = coords[0].cdelt2 / coords[0].cdelt1; + rotate = DEG_RAD*atan2 (coords[0].pc2_1*Lambda, coords[0].pc1_1); + fits_modify (header, "CROTA1", "%f", 1, rotate); + fits_modify (header, "CROTA2", "%f", 1, rotate); + return (TRUE); + } + + fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1); + fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2); + fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1); + fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2); + fits_modify (header, "NPLYTERM", "%d", 1, coords[0].Npolyterms); + + /* RA Terms */ + if (coords[0].Npolyterms > 1) { + fits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]); /* polyterms[0]); */ + fits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]); /* polyterms[1]); */ + fits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]); /* polyterms[2]); */ + } + if (coords[0].Npolyterms > 2) { + fits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]); /* polyterms[3]); */ + fits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]); /* polyterms[4]); */ + fits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]); /* polyterms[5]); */ + fits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]); /* polyterms[6]); */ + } + + /* Dec Terms */ + if (coords[0].Npolyterms > 1) { + fits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]); /* polyterms[7]); */ + fits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]); /* polyterms[8]); */ + fits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]); /* polyterms[9]); */ + } + if (coords[0].Npolyterms > 2) { + fits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]); /* polyterms[10]); */ + fits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]); /* polyterms[11]); */ + fits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]); /* polyterms[12]); */ + fits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]); /* polyterms[13]); */ + } + return (TRUE); +} + +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch) { + + double T; + double A, D, RA, DEC, zeta, z, theta; + double SA, CA, SD, CD; + + T = (out_epoch - in_epoch) / 100.0; + + zeta = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T); + theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T); + z = 0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T; + + A = *ra; + D = *dec; + SD = cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D); + CD = sqrt (1 - SD*SD); + SA = sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD; + CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD; + + DEC = DEG_RAD*asin(SD); + RA = DEG_RAD*atan2(SA, CA) + z; + + if (RA < 0) + RA += 360; + + *ra = RA; + *dec = DEC; +} + +/* -PLY projection is an extrapolation of the -TAN projection. + In addition to the usual linear terms of CRPIXi, PC00i00j, there are + higher order polynomial terms (up to 3rd order): + Axis 1 terms: + PCA1X2Y0 = coords.polyterm[0][0] = x^2 + PCA1X1Y1 = coords.polyterm[1][0] = xy + PCA1X0Y2 = coords.polyterm[2][0] = y^2 + PCA1X3Y0 = coords.polyterm[3][0] = x^3 + PCA1X2Y1 = coords.polyterm[4][0] = x^2 y + PCA1X1Y2 = coords.polyterm[5][0] = x y^2 + PCA1X0Y3 = coords.polyterm[6][0] = y^2 + Axis 2 terms: + PCA2X2Y0 = coords.polyterm[0][1] = x^2 + PCA2X1Y1 = coords.polyterm[1][1] = xy + PCA2X0Y2 = coords.polyterm[2][1] = y^2 + PCA2X3Y0 = coords.polyterm[3][1] = x^3 + PCA2X2Y1 = coords.polyterm[4][1] = x^2 y + PCA2X1Y2 = coords.polyterm[5][1] = x y^2 + PCA2X0Y3 = coords.polyterm[6][1] = y^2 +*/ + +# if (0) + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + if (Polynomi) { + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]); + Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]); + } + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]); + Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]); + } + } + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); +/** this code is the old method used for higher order terms. they + are essentially 6th order, with weird coupled terms. + I don't think any real data used these terms, but they should + be re-calculated, I would think +**/ + +# endif Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog.c (revision 5388) @@ -0,0 +1,130 @@ +# include +# include + +/* possible exit status for lock_catalog: + 0 - failure (including lock failure) + 1 - success + 2 - empty file (file may be open or closed!) +*/ + +int lock_catalog (Catalog *catalog, int lockmode) { + + int dbstate; + + /* record lockmode choice */ + catalog[0].lockmode = lockmode; + + /* set lock on database, create stream f */ + catalog[0].f = fsetlockfile (catalog[0].filename, 3600.0, catalog[0].lockmode, &dbstate); + + if (dbstate == LCK_MISSING) return (2); + if (dbstate == LCK_EMPTY) return (2); + if (catalog[0].f == NULL) return (0); + + fseek (catalog[0].f, 0, SEEK_SET); + return (1); +} + +int unlock_catalog (Catalog *catalog) { + + int dbstate; + + if (catalog[0].f == (FILE *) NULL) return (2); + fclearlockfile (catalog[0].filename, catalog[0].f, catalog[0].lockmode, &dbstate); + return (1); +} + +int load_catalog (Catalog *catalog, int VERBOSE) { + + int Naxis, split; + char measure[80]; + + /** we will have to check catmode for MYSQL before reading from the file **/ + + fits_fread_header (catalog[0].f, &catalog[0].header); + if (!fits_scan (&catalog[0].header, "NAXIS", "%d", 1, &Naxis)) { + if (VERBOSE) fprintf (stderr, "can't determine catalog db mode\n"); + return (FALSE); + } + split = fits_scan (&catalog[0].header, "MEASURE", "%s", 1, measure); + + /* currenly, the catmode information is carried per table. + the user-set value of catmode will set the output mode for new tables + an old table will keep its mode + */ + + catalog[0].catmode = DVO_MODE_MEF; + catalog[0].catformat = DVO_FORMAT_UNDEF; + if (split) catalog[0].catmode = DVO_MODE_SPLIT; + if (Naxis == 2) catalog[0].catmode = DVO_MODE_RAW; + /* don't use Naxis == 2 and split mode! */ + + switch (catalog[0].catmode) { + case DVO_MODE_RAW: + if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_RAW)\n"); + load_catalog_raw (catalog, VERBOSE); + break; + case DVO_MODE_MEF: + if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_MEF)\n"); + load_catalog_mef (catalog, VERBOSE); + break; + case DVO_MODE_SPLIT: + if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_SPLIT)\n"); + load_catalog_split (catalog, VERBOSE); + break; + default: + fprintf (stderr, "error getting catalog mode\n"); + exit (2); + } + return (TRUE); +} + +int save_catalog (Catalog *catalog, char VERBOSE) { + + switch (catalog[0].catmode) { + case DVO_MODE_RAW: + save_catalog_raw (catalog, VERBOSE); + break; + case DVO_MODE_MEF: + save_catalog_mef (catalog, VERBOSE); + break; + case DVO_MODE_SPLIT: + save_catalog_split (catalog, VERBOSE); + break; + default: + fprintf (stderr, "invalid catalog mode\n"); + exit (2); + } + return (TRUE); +} + +int update_catalog (Catalog *catalog, char VERBOSE) { + + /* update is only valid for catmode SPLIT */ + switch (catalog[0].catmode) { + case DVO_MODE_RAW: + save_catalog_raw (catalog, VERBOSE); + break; + case DVO_MODE_MEF: + save_catalog_mef (catalog, VERBOSE); + break; + case DVO_MODE_SPLIT: + /* new file needs to use save_catalog_split */ + if (catalog[0].Nave_disk == 0) { + save_catalog_split (catalog, VERBOSE); + } else { + update_catalog_split (catalog, VERBOSE); + } + break; + default: + fprintf (stderr, "invalid catalog mode\n"); + exit (2); + } + return (TRUE); +} + +/* + mode : items to read (LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF) + format : what structure on disk (INTERNAL, ELIXIR, LONEOS, PANSTARRS) + style : raw, mef, split, mysql +*/ Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_mef.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 5388) @@ -0,0 +1,236 @@ +# include + +int load_catalog_mef (Catalog *catalog, int VERBOSE) { + + int Nitems, Nbytes, Nexpect, Naverage, Nmeasure, Nmissing, Nsecfilt; + FILE *f; + + Header header; + Matrix matrix; + FTable ftable; + + f = catalog[0].f; + ftable.header = &header; + + /* move pointer past header -- must be already read (load_catalog) */ + fseek (f, catalog[0].header.size, SEEK_SET); + + /* matrix should be empty */ + if (!fits_fread_matrix (f, &matrix, &catalog[0].header)) { + if (VERBOSE) fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + /* get the components from the header */ + if (!fits_scan (&catalog[0].header, "NSTARS", "%d", 1, &Naverage)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMEAS", "%d", 1, &Nmeasure)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMISS", "%d", 1, &Nmissing)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NSECFILT", "%d", 1, &Nsecfilt)) Nsecfilt = 0; + + /** Nsecfilt is unusual: it does not list the number of data items in the table + instead, the number of items is Nsecfilt * Naverage. **/ + catalog[0].Nsecfilt = Nsecfilt; + + /* validate table mode */ + + /* read Average table header */ + if (!fits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table average header"); + return (FALSE); + } + /* read Average table data (or skip) */ + if (catalog[0].catflags & LOAD_AVES) { + if (!fits_fread_ftable_data (f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table average data"); + return (FALSE); + } + catalog[0].average = FtableToAverage (&ftable, &catalog[0].Naverage, &catalog[0].catformat); + if (Naverage != catalog[0].Naverage) { + fprintf (stderr, "Warning: mismatch between Naverage in PHU and Table headers (%d vs %d)\n", Naverage, catalog[0].Naverage); + } + } else { + Nbytes = fits_matrix_size (&header); + fseek (f, Nbytes, SEEK_CUR); + } + fits_free_header (&header); + /** free the ftable or not? data is being used still..? **/ + + /* read Measure table header */ + if (!fits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table measure header"); + return (FALSE); + } + /* read Measure table data */ + if (catalog[0].catflags & LOAD_MEAS) { + if (!fits_fread_ftable_data (f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table measure data"); + return (FALSE); + } + catalog[0].measure = FtableToMeasure (&ftable, &catalog[0].Nmeasure, &catalog[0].catformat); + if (Nmeasure != catalog[0].Nmeasure) { + fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers (%d vs %d)\n", Nmeasure, catalog[0].Nmeasure); + } + } else { + Nbytes = fits_matrix_size (&header); + fseek (f, Nbytes, SEEK_CUR); + } + + /* read Missing table header */ + if (!fits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table missing header"); + return (FALSE); + } + /* read Missing table data */ + if (catalog[0].catflags & LOAD_MISS) { + if (!fits_fread_ftable_data (f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table missing data"); + return (FALSE); + } + /* no conversions currently defined */ + catalog[0].missing = fits_table_get_Missing (&ftable, &catalog[0].Nmissing, NULL); + if (Nmissing != catalog[0].Nmissing) { + fprintf (stderr, "Warning: mismatch between Nmissing in PHU and Table headers (%d vs %d)\n", Nmissing, catalog[0].Nmissing); + } + } else { + Nbytes = fits_matrix_size (&header); + fseek (f, Nbytes, SEEK_CUR); + } + + // catalog[0].secfilt = fits_table_get_SecFilt (&ftable, &Nitems, NULL); + // catalog[0].Nsecfilt = Nitems / catalog[0].Naverage; + + /* read secfilt table header */ + if (!fits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table secfilt header"); + return (FALSE); + } + /* read secfilt table data */ + if (catalog[0].catflags & LOAD_SECF) { + if (!fits_fread_ftable_data (f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table secfilt data"); + return (FALSE); + } + Nexpect = catalog[0].Nsecfilt * catalog[0].Naverage; + catalog[0].secfilt = FtableToSecFilt (&ftable, &Nitems, &catalog[0].catformat); + if (Nexpect != Nitems) { + fprintf (stderr, "Warning: mismatch between Nsecfilt items in PHU and Table headers (%d vs %d)\n", Nexpect, Nitems); + } + } else { + /* no real need to skip the data array here... */ + Nbytes = fits_matrix_size (&header); + fseek (f, Nbytes, SEEK_CUR); + } + return (TRUE); +} + +/* save_catalog_mef writes a complete new file from scratch */ + +int save_catalog_mef (Catalog *catalog, char VERBOSE) { + + int Nitems; + FILE *f; + Matrix matrix; + Header header; + FTable ftable; + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (TRUE); + } + + /* make sure header is consistent with data */ + fits_modify (&catalog[0].header, "NSTARS", "%d", 1, catalog[0].Naverage); + fits_modify (&catalog[0].header, "NMEAS", "%d", 1, catalog[0].Nmeasure); + fits_modify (&catalog[0].header, "NMISS", "%d", 1, catalog[0].Nmissing); + fits_modify (&catalog[0].header, "NSECFILT", "%d", 1, catalog[0].Nsecfilt); + fits_modify (&catalog[0].header, "EXTEND", "%t", 1, TRUE); + + f = catalog[0].f; + /* rewind file pointers and truncate */ + fseek (f, 0, SEEK_SET); + ftruncate (fileno (catalog[0].f), 0); + ftable.header = &header; + + /* write table PHU header */ + if (!fits_fwrite_header (catalog[0].f, &catalog[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&catalog[0].header, &matrix); + if (!fits_fwrite_matrix (catalog[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out Average table (convert to FITS table format) */ + AverageToFtable (&ftable, catalog[0].average, catalog[0].Naverage, catalog[0].catformat); + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + + /* write out Measure table (convert to FITS table format) */ + MeasureToFtable (&ftable, catalog[0].measure, catalog[0].Nmeasure, catalog[0].catformat); + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + + /* write out Missing table (convert to FITS table format) */ + fits_table_set_Missing (&ftable, catalog[0].missing, catalog[0].Nmissing); + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + + /* write out SecFilt table (convert to FITS table format) */ + Nitems = catalog[0].Naverage * catalog[0].Nsecfilt; + SecFiltToFtable (&ftable, catalog[0].secfilt, Nitems, catalog[0].catformat); + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + + return (TRUE); +} + +/* + catalog data is: + header (bitpix == 8) + matrix (empty) + average header + average table + measure header + measure table + missing header + missing table + secfilt header + secfilt table +*/ + Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_raw.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 5388) @@ -0,0 +1,580 @@ +# include + +/* read data from raw-style catalog file; set data format values based on file */ + +int load_catalog_raw (Catalog *catalog, int VERBOSE) { + + int Nitems, nitems; + int i, Nmeas, Nmiss, size; + int NewMeasure, Nskip; + int AverageSize, MeasureSize, MissingSize, SecFiltSize; + FILE *f; + struct stat filestatus; + char format[80], telescope[80]; + + f = catalog[0].f; + + /* move pointer past header -- must be already read (load_catalog) */ + fseek (f, catalog[0].header.size, SEEK_SET); + + /* get the components from the header */ + catalog[0].Naverage = catalog[0].Nmeasure = catalog[0].Nmissing = catalog[0].Nsecfilt = 0; + if (!fits_scan (&catalog[0].header, "NSTARS", "%d", 1, &catalog[0].Naverage)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMEAS", "%d", 1, &catalog[0].Nmeasure)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMISS", "%d", 1, &catalog[0].Nmissing)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NSECFILT", "%d", 1, &catalog[0].Nsecfilt)) catalog[0].Nsecfilt = 0; + + /* determine catalog format */ + catalog[0].catformat = DVO_FORMAT_UNDEF; + if (fits_scan (&catalog[0].header, "FORMAT", "%s", 1, format)) { + if (!strcmp (format, "INTERNAL")) catalog[0].catformat = DVO_FORMAT_INTERNAL; + if (!strcmp (format, "LONEOS")) catalog[0].catformat = DVO_FORMAT_LONEOS; + if (!strcmp (format, "ELIXIR")) catalog[0].catformat = DVO_FORMAT_ELIXIR; + if (!strcmp (format, "PANSTARRS")) catalog[0].catformat = DVO_FORMAT_PANSTARRS; + if (catalog[0].catformat != DVO_FORMAT_UNDEF) goto got_format; + } + if (fits_scan (&catalog[0].header, "NEWMEAS", "%t", 1, &NewMeasure)) { + catalog[0].catformat = DVO_FORMAT_ELIXIR; + goto got_format; + } + if (fits_scan (&catalog[0].header, "TELESCOP", "%s", 1, telescope)) { + if (!strncmp (telescope, "LONEOS", strlen("LONEOS"))) { + catalog[0].catformat = DVO_FORMAT_LONEOS; + goto got_format; + } + if (!strncmp (telescope, "1.3m McGraw-Hill", strlen("1.3m McGraw-Hill"))) { + catalog[0].catformat = DVO_FORMAT_ELIXIR; + goto got_format; + } + } + if (VERBOSE) fprintf (stderr, "cannot determine catalog format\n"); + return (FALSE); + +got_format: + /* determine datatype sizes */ + switch (catalog[0].catformat) { + case DVO_FORMAT_INTERNAL: + AverageSize = sizeof(Average); + MeasureSize = sizeof(Measure); + SecFiltSize = sizeof (SecFilt); + break; + case DVO_FORMAT_LONEOS: + AverageSize = sizeof(AverageLoneos); + MeasureSize = sizeof(MeasureLoneos); + SecFiltSize = sizeof (SecFiltLoneos); + break; + case DVO_FORMAT_ELIXIR: + AverageSize = sizeof(AverageElixir); + MeasureSize = sizeof(MeasureElixir); + SecFiltSize = sizeof (SecFiltElixir); + break; + case DVO_FORMAT_PANSTARRS: + AverageSize = sizeof(AveragePanstarrs); + MeasureSize = sizeof(MeasurePanstarrs); + SecFiltSize = sizeof (SecFiltPanstarrs); + break; + default: + fprintf (stderr, "programming error in phot_catalog_raw\n"); + exit (2); + } + MissingSize = sizeof (Missing); + + /* predicted file size - for double checking data validity */ + size = catalog[0].header.size; + size += AverageSize * catalog[0].Naverage; + size += MeasureSize * catalog[0].Nmeasure; + size += MissingSize * catalog[0].Nmissing; + size += SecFiltSize * catalog[0].Nsecfilt * catalog[0].Naverage; + + /* check that file size makes sense */ + if (stat (catalog[0].filename, &filestatus) == -1) { + if (VERBOSE) fprintf (stderr, "failed to get status of catalog\n"); + return (FALSE); + } + if (size > filestatus.st_size) { + if (VERBOSE) { + fprintf (stderr, "star catalog has inconsistent size\n"); + fprintf (stderr, "average: %d = %d bytes\n", catalog[0].Naverage, catalog[0].Naverage*AverageSize); + fprintf (stderr, "measure: %d = %d bytes\n", catalog[0].Nmeasure, catalog[0].Nmeasure*MeasureSize); + fprintf (stderr, "missing: %d = %d bytes\n", catalog[0].Nmissing, catalog[0].Nmissing*MissingSize); + fprintf (stderr, "secfilt: %d = %d bytes\n", catalog[0].Nsecfilt, catalog[0].Nsecfilt*SecFiltSize*catalog[0].Naverage); + fprintf (stderr, "expect: %d, found: %d\n", size, (int)filestatus.st_size); + } + return (FALSE); + } + if (size < filestatus.st_size) { + if (VERBOSE) fprintf (stderr, "warning: file larger than expected\n"); + } + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars yet in catalog %s\n", catalog[0].filename); + return (TRUE); + } + + /* read and convert the averages (use a macro to clean this up?) */ + if (catalog[0].catflags & LOAD_AVES) { + catalog[0].average = ReadRawAverage (catalog[0].f, catalog[0].Naverage, catalog[0].catformat); + } else { + /* skip over averages */ + Nskip = catalog[0].Naverage * AverageSize; + fseek (f, Nskip, SEEK_CUR); + } + + /* read and convert the measures (use a macro to clean this up?) */ + if (catalog[0].catflags & LOAD_MEAS) { + catalog[0].measure = ReadRawMeasure (catalog[0].f, catalog[0].Nmeasure, catalog[0].catformat); + } else { + /* skip over measures */ + Nskip = catalog[0].Nmeasure * MeasureSize; + fseek (f, Nskip, SEEK_CUR); + } + + /* read and convert missing */ + if (catalog[0].catflags & LOAD_MISS) { + ALLOCATE (catalog[0].missing, Missing, MAX (catalog[0].Nmissing, 1)); + Nitems = catalog[0].Nmissing; + nitems = fread (catalog[0].missing, MissingSize, Nitems, f); + if (nitems != Nitems) { + if (VERBOSE) fprintf (stderr, "failed to read missing from catalog file %s (%d vs %d)\n", catalog[0].filename, nitems, Nitems); + return (FALSE); + } + fits_convert_Missing (catalog[0].missing, MissingSize, Nitems); + } else { + /* skip over missings */ + Nskip = catalog[0].Nmissing * MissingSize; + fseek (f, Nskip, SEEK_CUR); + } + + /* read and convert secfilt */ + if (catalog[0].catflags & LOAD_SECF) { + Nitems = catalog[0].Naverage * catalog[0].Nsecfilt; + catalog[0].secfilt = ReadRawSecFilt (catalog[0].f, Nitems, catalog[0].catformat); + } else { + /* skip over secfilts */ + Nskip = catalog[0].Nsecfilt * catalog[0].Naverage * SecFiltSize; + fseek (f, Nskip, SEEK_CUR); + } + + if (VERBOSE) fprintf (stderr, "read %d stars from catalog file %s (%d measurements, %d missing, %d secondary filters)\n", + catalog[0].Naverage, catalog[0].filename, catalog[0].Nmeasure, catalog[0].Nmissing, catalog[0].Nsecfilt); + + /* check data integrity */ + if (catalog[0].catflags & LOAD_AVES) { + for (i = Nmeas = Nmiss = 0; i < catalog[0].Naverage; i++) { + Nmeas += catalog[0].average[i].Nm; + Nmiss += catalog[0].average[i].Nn; + } + if ((Nmeas != catalog[0].Nmeasure) || (Nmiss != catalog[0].Nmissing)) { + if (VERBOSE) { + fprintf (stderr, "****** data in catalog %s is corrupt, sums don't check\n", catalog[0].filename); + fprintf (stderr, "****** Nmeas: %d, %d\n", Nmeas, catalog[0].Nmeasure); + fprintf (stderr, "****** Nmiss: %d, %d\n", Nmiss, catalog[0].Nmissing); + } + return (FALSE); + } + } + return (TRUE); +} + +int save_catalog_raw (Catalog *catalog, char VERBOSE) { + + int Nitems, nitems; + FILE *f; + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (TRUE); + } + + /* make sure header is consistent with data */ + fits_modify (&catalog[0].header, "NSTARS", "%d", 1, catalog[0].Naverage); + fits_modify (&catalog[0].header, "NMEAS", "%d", 1, catalog[0].Nmeasure); + fits_modify (&catalog[0].header, "NMISS", "%d", 1, catalog[0].Nmissing); + fits_modify (&catalog[0].header, "NSECFILT", "%d", 1, catalog[0].Nsecfilt); + + /* specify the appropriate data format */ + if (catalog[0].catformat == DVO_FORMAT_INTERNAL) fits_modify (&catalog[0].header, "FORMAT", "%s", 1, "INTERNAL"); + if (catalog[0].catformat == DVO_FORMAT_LONEOS) fits_modify (&catalog[0].header, "FORMAT", "%s", 1, "LONEOS"); + if (catalog[0].catformat == DVO_FORMAT_ELIXIR) fits_modify (&catalog[0].header, "FORMAT", "%s", 1, "ELIXIR"); + if (catalog[0].catformat == DVO_FORMAT_PANSTARRS) fits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PANSTARRS"); + + /* rewind file pointers and truncate file */ + f = catalog[0].f; + fseek (f, 0, SEEK_SET); + ftruncate (fileno (catalog[0].f), 0); + + /* write header data (use fits_write_header?) */ + nitems = fwrite (catalog[0].header.buffer, 1, catalog[0].header.size, f); + if (nitems != catalog[0].header.size) { + if (VERBOSE) fprintf (stderr, "failed to write header\n"); + return (FALSE); + } + + /* write averages and measures */ + WriteRawAverage (f, catalog[0].average, catalog[0].Naverage, catalog[0].catformat); + WriteRawMeasure (f, catalog[0].measure, catalog[0].Nmeasure, catalog[0].catformat); + + /* write missing data */ + Nitems = catalog[0].Nmissing; + fits_convert_Missing (catalog[0].missing, sizeof(Missing), Nitems); + nitems = fwrite (catalog[0].missing, sizeof(Missing), Nitems, f); + if (nitems != Nitems) { + if (VERBOSE) fprintf (stderr, "failed to write catalog file missing %s\n", catalog[0].filename); + return (FALSE); + } + + Nitems = catalog[0].Naverage * catalog[0].Nsecfilt; + WriteRawSecFilt (f, catalog[0].secfilt, Nitems, catalog[0].catformat); + return (TRUE); +} + +/*** + some warnings: + - need to be wary of exit condition errors +***/ + +/** Average / Raw Table conversions **/ + +Average *ReadRawAverage (FILE *f, int Naverage, int format) { + + int nitems; + Average *average; + AverageElixir *tmpAverageElixir; + AverageLoneos *tmpAverageLoneos; + AveragePanstarrs *tmpAveragePanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + ALLOCATE (average, Average, MAX (Naverage, 1)); + nitems = fread (average, sizeof(Average), Naverage, f); + if (nitems != Naverage) { + fprintf (stderr, "failed to read averages (%d vs %d)\n", nitems, Naverage); + return (NULL); + } + fits_convert_Average (average, sizeof(Average), Naverage); + break; + case DVO_FORMAT_ELIXIR: + ALLOCATE (tmpAverageElixir, AverageElixir, MAX (Naverage, 1)); + nitems = fread (tmpAverageElixir, sizeof(AverageElixir), Naverage, f); + if (nitems != Naverage) { + fprintf (stderr, "failed to read averages (%d vs %d)\n", nitems, Naverage); + return (NULL); + } + fits_convert_AverageElixir (tmpAverageElixir, sizeof(AverageElixir), Naverage); + average = AverageElixirToInternal (tmpAverageElixir, Naverage); + free (tmpAverageElixir); + break; + case DVO_FORMAT_LONEOS: + ALLOCATE (tmpAverageLoneos, AverageLoneos, MAX (Naverage, 1)); + nitems = fread (tmpAverageLoneos, sizeof(AverageLoneos), Naverage, f); + if (nitems != Naverage) { + fprintf (stderr, "failed to read averages (%d vs %d)\n", nitems, Naverage); + return (NULL); + } + fits_convert_AverageLoneos (tmpAverageLoneos, sizeof(AverageLoneos), Naverage); + average = AverageLoneosToInternal (tmpAverageLoneos, Naverage); + free (tmpAverageLoneos); + break; + case DVO_FORMAT_PANSTARRS: + ALLOCATE (tmpAveragePanstarrs, AveragePanstarrs, MAX (Naverage, 1)); + nitems = fread (tmpAveragePanstarrs, sizeof(AveragePanstarrs), Naverage, f); + if (nitems != Naverage) { + fprintf (stderr, "failed to read averages (%d vs %d)\n", nitems, Naverage); + return (NULL); + } + fits_convert_AveragePanstarrs (tmpAveragePanstarrs, sizeof(AveragePanstarrs), Naverage); + average = AveragePanstarrsToInternal (tmpAveragePanstarrs, Naverage); + free (tmpAveragePanstarrs); + break; + default: + fprintf (stderr, "error reading measures\n"); + return (NULL); + } + return (average); +} + +/* accepts and converts internal average formats and outputs + raw data in the specified format */ +int WriteRawAverage (FILE *f, Average *average, int Naverage, int format) { + + int nitems; + AverageElixir *tmpAverageElixir; + AverageLoneos *tmpAverageLoneos; + AveragePanstarrs *tmpAveragePanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_convert_Average (average, sizeof(Average), Naverage); + nitems = fwrite (average, sizeof(Average), Naverage, f); + if (nitems != Naverage) { + fprintf (stderr, "failed to write averages (%d vs %d)\n", nitems, Naverage); + return (FALSE); + } + break; + case DVO_FORMAT_ELIXIR: + tmpAverageElixir = AverageInternalToElixir (average, Naverage); + fits_convert_AverageElixir (tmpAverageElixir, sizeof(AverageElixir), Naverage); + nitems = fwrite (tmpAverageElixir, sizeof(AverageElixir), Naverage, f); + free (tmpAverageElixir); + if (nitems != Naverage) { + fprintf (stderr, "failed to write averages (%d vs %d)\n", nitems, Naverage); + return (FALSE); + } + break; + case DVO_FORMAT_LONEOS: + tmpAverageLoneos = AverageInternalToLoneos (average, Naverage); + fits_convert_AverageLoneos (tmpAverageLoneos, sizeof(AverageLoneos), Naverage); + nitems = fwrite (tmpAverageLoneos, sizeof(AverageLoneos), Naverage, f); + free (tmpAverageLoneos); + if (nitems != Naverage) { + fprintf (stderr, "failed to write averages (%d vs %d)\n", nitems, Naverage); + return (FALSE); + } + break; + case DVO_FORMAT_PANSTARRS: + tmpAveragePanstarrs = AverageInternalToPanstarrs (average, Naverage); + fits_convert_AveragePanstarrs (tmpAveragePanstarrs, sizeof(AveragePanstarrs), Naverage); + nitems = fwrite (tmpAveragePanstarrs, sizeof(AveragePanstarrs), Naverage, f); + free (tmpAveragePanstarrs); + if (nitems != Naverage) { + fprintf (stderr, "failed to write averages (%d vs %d)\n", nitems, Naverage); + return (FALSE); + } + break; + default: + fprintf (stderr, "error writing averages\n"); + return (FALSE); + } + return (TRUE); +} + +/** Average / Raw Table conversions **/ + +Measure *ReadRawMeasure (FILE *f, int Nmeasure, int format) { + + int nitems; + Measure *measure; + MeasureElixir *tmpMeasureElixir; + MeasureLoneos *tmpMeasureLoneos; + MeasurePanstarrs *tmpMeasurePanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + ALLOCATE (measure, Measure, MAX (Nmeasure, 1)); + nitems = fread (measure, sizeof(Measure), Nmeasure, f); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to read measures (%d vs %d)\n", nitems, Nmeasure); + return (NULL); + } + fits_convert_Measure (measure, sizeof(Measure), Nmeasure); + break; + case DVO_FORMAT_ELIXIR: + ALLOCATE (tmpMeasureElixir, MeasureElixir, MAX (Nmeasure, 1)); + nitems = fread (tmpMeasureElixir, sizeof(MeasureElixir), Nmeasure, f); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to read measures (%d vs %d)\n", nitems, Nmeasure); + return (NULL); + } + fits_convert_MeasureElixir (tmpMeasureElixir, sizeof(MeasureElixir), Nmeasure); + measure = MeasureElixirToInternal (tmpMeasureElixir, Nmeasure); + free (tmpMeasureElixir); + break; + case DVO_FORMAT_LONEOS: + ALLOCATE (tmpMeasureLoneos, MeasureLoneos, MAX (Nmeasure, 1)); + nitems = fread (tmpMeasureLoneos, sizeof(MeasureLoneos), Nmeasure, f); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to read measures (%d vs %d)\n", nitems, Nmeasure); + return (NULL); + } + fits_convert_MeasureLoneos (tmpMeasureLoneos, sizeof(MeasureLoneos), Nmeasure); + measure = MeasureLoneosToInternal (tmpMeasureLoneos, Nmeasure); + free (tmpMeasureLoneos); + break; + case DVO_FORMAT_PANSTARRS: + ALLOCATE (tmpMeasurePanstarrs, MeasurePanstarrs, MAX (Nmeasure, 1)); + nitems = fread (tmpMeasurePanstarrs, sizeof(MeasurePanstarrs), Nmeasure, f); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to read measures (%d vs %d)\n", nitems, Nmeasure); + return (NULL); + } + fits_convert_MeasurePanstarrs (tmpMeasurePanstarrs, sizeof(MeasurePanstarrs), Nmeasure); + measure = MeasurePanstarrsToInternal (tmpMeasurePanstarrs, Nmeasure); + free (tmpMeasurePanstarrs); + break; + default: + fprintf (stderr, "error reading measures\n"); + return (NULL); + } + return (measure); +} + +/* accepts and converts internal measure formats and outputs + raw data in the specified format */ +int WriteRawMeasure (FILE *f, Measure *measure, int Nmeasure, int format) { + + int nitems; + MeasureElixir *tmpMeasureElixir; + MeasureLoneos *tmpMeasureLoneos; + MeasurePanstarrs *tmpMeasurePanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_convert_Measure (measure, sizeof(Measure), Nmeasure); + nitems = fwrite (measure, sizeof(Measure), Nmeasure, f); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to write measures (%d vs %d)\n", nitems, Nmeasure); + return (FALSE); + } + break; + case DVO_FORMAT_ELIXIR: + tmpMeasureElixir = MeasureInternalToElixir (measure, Nmeasure); + fits_convert_MeasureElixir (tmpMeasureElixir, sizeof(MeasureElixir), Nmeasure); + nitems = fwrite (tmpMeasureElixir, sizeof(MeasureElixir), Nmeasure, f); + free (tmpMeasureElixir); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to write measures (%d vs %d)\n", nitems, Nmeasure); + return (FALSE); + } + break; + case DVO_FORMAT_LONEOS: + tmpMeasureLoneos = MeasureInternalToLoneos (measure, Nmeasure); + fits_convert_MeasureLoneos (tmpMeasureLoneos, sizeof(MeasureLoneos), Nmeasure); + nitems = fwrite (tmpMeasureLoneos, sizeof(MeasureLoneos), Nmeasure, f); + free (tmpMeasureLoneos); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to write measures (%d vs %d)\n", nitems, Nmeasure); + return (FALSE); + } + break; + case DVO_FORMAT_PANSTARRS: + tmpMeasurePanstarrs = MeasureInternalToPanstarrs (measure, Nmeasure); + fits_convert_MeasurePanstarrs (tmpMeasurePanstarrs, sizeof(MeasurePanstarrs), Nmeasure); + nitems = fwrite (tmpMeasurePanstarrs, sizeof(MeasurePanstarrs), Nmeasure, f); + free (tmpMeasurePanstarrs); + if (nitems != Nmeasure) { + fprintf (stderr, "failed to write measures (%d vs %d)\n", nitems, Nmeasure); + return (FALSE); + } + break; + default: + fprintf (stderr, "error writing measures\n"); + return (FALSE); + } + return (TRUE); +} + +/** SecFilt / Raw Table conversions **/ + +SecFilt *ReadRawSecFilt (FILE *f, int Nsecfilt, int format) { + + int nitems; + SecFilt *secfilt; + SecFiltElixir *tmpSecFiltElixir; + SecFiltLoneos *tmpSecFiltLoneos; + SecFiltPanstarrs *tmpSecFiltPanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + ALLOCATE (secfilt, SecFilt, MAX (Nsecfilt, 1)); + nitems = fread (secfilt, sizeof(SecFilt), Nsecfilt, f); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to read secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (NULL); + } + fits_convert_SecFilt (secfilt, sizeof(SecFilt), Nsecfilt); + break; + case DVO_FORMAT_ELIXIR: + ALLOCATE (tmpSecFiltElixir, SecFiltElixir, MAX (Nsecfilt, 1)); + nitems = fread (tmpSecFiltElixir, sizeof(SecFiltElixir), Nsecfilt, f); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to read secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (NULL); + } + fits_convert_SecFiltElixir (tmpSecFiltElixir, sizeof(SecFiltElixir), Nsecfilt); + secfilt = SecFiltElixirToInternal (tmpSecFiltElixir, Nsecfilt); + free (tmpSecFiltElixir); + break; + case DVO_FORMAT_LONEOS: + ALLOCATE (tmpSecFiltLoneos, SecFiltLoneos, MAX (Nsecfilt, 1)); + nitems = fread (tmpSecFiltLoneos, sizeof(SecFiltLoneos), Nsecfilt, f); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to read secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (NULL); + } + fits_convert_SecFiltLoneos (tmpSecFiltLoneos, sizeof(SecFiltLoneos), Nsecfilt); + secfilt = SecFiltLoneosToInternal (tmpSecFiltLoneos, Nsecfilt); + free (tmpSecFiltLoneos); + break; + case DVO_FORMAT_PANSTARRS: + ALLOCATE (tmpSecFiltPanstarrs, SecFiltPanstarrs, MAX (Nsecfilt, 1)); + nitems = fread (tmpSecFiltPanstarrs, sizeof(SecFiltPanstarrs), Nsecfilt, f); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to read secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (NULL); + } + fits_convert_SecFiltPanstarrs (tmpSecFiltPanstarrs, sizeof(SecFiltPanstarrs), Nsecfilt); + secfilt = SecFiltPanstarrsToInternal (tmpSecFiltPanstarrs, Nsecfilt); + free (tmpSecFiltPanstarrs); + break; + default: + fprintf (stderr, "error reading measures\n"); + return (NULL); + } + return (secfilt); +} + +/* accepts and converts internal secfilt formats and outputs + raw data in the specified format */ +int WriteRawSecFilt (FILE *f, SecFilt *secfilt, int Nsecfilt, int format) { + + int nitems; + SecFiltElixir *tmpSecFiltElixir; + SecFiltLoneos *tmpSecFiltLoneos; + SecFiltPanstarrs *tmpSecFiltPanstarrs; + + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_convert_SecFilt (secfilt, sizeof(SecFilt), Nsecfilt); + nitems = fwrite (secfilt, sizeof(SecFilt), Nsecfilt, f); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to write secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (FALSE); + } + break; + case DVO_FORMAT_ELIXIR: + tmpSecFiltElixir = SecFiltInternalToElixir (secfilt, Nsecfilt); + fits_convert_SecFiltElixir (tmpSecFiltElixir, sizeof(SecFiltElixir), Nsecfilt); + nitems = fwrite (tmpSecFiltElixir, sizeof(SecFiltElixir), Nsecfilt, f); + free (tmpSecFiltElixir); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to write secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (FALSE); + } + break; + case DVO_FORMAT_LONEOS: + tmpSecFiltLoneos = SecFiltInternalToLoneos (secfilt, Nsecfilt); + fits_convert_SecFiltLoneos (tmpSecFiltLoneos, sizeof(SecFiltLoneos), Nsecfilt); + nitems = fwrite (tmpSecFiltLoneos, sizeof(SecFiltLoneos), Nsecfilt, f); + free (tmpSecFiltLoneos); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to write secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (FALSE); + } + break; + case DVO_FORMAT_PANSTARRS: + tmpSecFiltPanstarrs = SecFiltInternalToPanstarrs (secfilt, Nsecfilt); + fits_convert_SecFiltPanstarrs (tmpSecFiltPanstarrs, sizeof(SecFiltPanstarrs), Nsecfilt); + nitems = fwrite (tmpSecFiltPanstarrs, sizeof(SecFiltPanstarrs), Nsecfilt, f); + free (tmpSecFiltPanstarrs); + if (nitems != Nsecfilt) { + fprintf (stderr, "failed to write secfilts (%d vs %d)\n", nitems, Nsecfilt); + return (FALSE); + } + break; + default: + fprintf (stderr, "error writing secfilts\n"); + return (FALSE); + } + return (TRUE); +} + Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_split.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 5388) @@ -0,0 +1,639 @@ +# include + +int load_catalog_split (Catalog *catalog, int VERBOSE) { + + int Nitems, Nexpect, Naverage, Nmeasure, Nmissing, Nsecfilt; + char *path, string[80]; + Header header; + Matrix matrix; + FTable ftable; + Catalog *measure, *missing, *secfilt; + + /* ftable header storage for below */ + ftable.header = &header; + + /* needed to find the split files below */ + path = pathname (catalog[0].filename); + + /* get the components from the header - these duplicate information in the split files (NAXIS2) */ + if (!fits_scan (&catalog[0].header, "NSTARS", "%d", 1, &Naverage)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMEAS", "%d", 1, &Nmeasure)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NMISS", "%d", 1, &Nmissing)) return (FALSE); + if (!fits_scan (&catalog[0].header, "NSECFILT", "%d", 1, &Nsecfilt)) Nsecfilt = 0; + + /** Nsecfilt is unusual: it does not list the number of data items in the table + instead, the number of items is Nsecfilt * Naverage. **/ + catalog[0].Nsecfilt = Nsecfilt; + + /* default values */ + catalog[0].average = NULL; + catalog[0].measure = NULL; + catalog[0].missing = NULL; + catalog[0].secfilt = NULL; + + /* XXX EAM : validate table mode */ + + /*** Average Table ***/ + + if (catalog[0].catflags & LOAD_AVES) { + /* move pointer past header -- must be already read (load_catalog) */ + fseek (catalog[0].f, catalog[0].header.size, SEEK_SET); + /* matrix should be empty */ + if (!fits_fread_matrix (catalog[0].f, &matrix, &catalog[0].header)) { + if (VERBOSE) fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + /* read Average table header */ + if (!fits_fread_header (catalog[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table average header"); + return (FALSE); + } + /* read Average table data */ + if (!fits_fread_ftable_data (catalog[0].f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table average data"); + return (FALSE); + } + catalog[0].average = FtableToAverage (&ftable, &catalog[0].Naverage, &catalog[0].catformat); + if (Naverage != catalog[0].Naverage) { + fprintf (stderr, "Warning: mismatch between Naverage in PHU and Table headers (%d vs %d)\n", Naverage, catalog[0].Naverage); + } + fits_free_header (&header); + } + + /*** Measure Table ***/ + + measure = NULL; + + /* (Full Load) */ + if (catalog[0].catflags & LOAD_MEAS) { + ALLOCATE (measure, Catalog, 1); + + /* get split filename from main header (paths relative to cpt file) */ + if (!fits_scan (&catalog[0].header, "MEASURE", "%s", 1, string)) return (FALSE); + ALLOCATE (measure[0].filename, char, strlen(path) + strlen(string) + 2); + sprintf (measure[0].filename, "%s/%s", path, string); + + /* lock & open catalog file */ + if (lock_catalog (measure, catalog[0].lockmode) != 1) { + fprintf (stderr, "cannot access measure file %s\n", measure[0].filename); + exit (2); + } + + /* read PHU */ + if (!fits_load_header (measure[0].f, &measure[0].header)) { + if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", measure[0].filename); + return (FALSE); + } + /* matrix should be empty */ + if (!fits_fread_matrix (measure[0].f, &matrix, &measure[0].header)) { + if (VERBOSE) fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + /* read Measure table header */ + if (!fits_fread_header (measure[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read measure PHU header"); + return (FALSE); + } + /* read Measure table data */ + if (!fits_fread_ftable_data (measure[0].f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table measure data"); + return (FALSE); + } + /* convert data format to internal */ + catalog[0].measure = FtableToMeasure (&ftable, &catalog[0].Nmeasure, &catalog[0].catformat); + if (Nmeasure != catalog[0].Nmeasure) { + fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers (%d vs %d)\n", Nmeasure, catalog[0].Nmeasure); + } + catalog[0].Nmeas_off = 0; + fits_free_header (&header); + fits_free_matrix (&matrix); + } + + /* (Meta Load) */ + if (catalog[0].catflags & LOAD_MEAS_META) { + ALLOCATE (measure, Catalog, 1); + + /* get split filename from main header (paths relative to cpt file) */ + if (!fits_scan (&catalog[0].header, "MEASURE", "%s", 1, string)) return (FALSE); + ALLOCATE (measure[0].filename, char, strlen(path) + strlen(string) + 2); + sprintf (measure[0].filename, "%s/%s", path, string); + + /* lock & open catalog file */ + if (lock_catalog (measure, catalog[0].lockmode) != 1) { + fprintf (stderr, "cannot access measure file %s\n", measure[0].filename); + exit (2); + } + + /* read PHU */ + if (!fits_load_header (measure[0].f, &measure[0].header)) { + if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", measure[0].filename); + return (FALSE); + } + + /* allocate dummy data for valid realloc */ + ALLOCATE (catalog[0].measure, Measure, 1); + + /* set up the table size information */ + catalog[0].Nmeasure = 0; /* no rows loaded */ + catalog[0].Nmeas_off = Nmeasure; /* start of new entries */ + } + catalog[0].measure_catalog = measure; + + /*** Missing Table ***/ + + missing = NULL; + if (catalog[0].catflags & LOAD_MISS) { + ALLOCATE (missing, Catalog, 1); + + /* get split filename from main header (paths relative to cpt file) */ + if (!fits_scan (&catalog[0].header, "MISSING", "%s", 1, string)) return (FALSE); + ALLOCATE (missing[0].filename, char, strlen(path) + strlen(string) + 2); + sprintf (missing[0].filename, "%s/%s", path, string); + + /* lock & open catalog file */ + if (lock_catalog (missing, catalog[0].lockmode) != 1) { + fprintf (stderr, "cannot access missing file %s\n", missing[0].filename); + exit (2); + } + + /* read PHU */ + if (!fits_load_header (missing[0].f, &missing[0].header)) { + if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", missing[0].filename); + return (FALSE); + } + /* matrix should be empty */ + if (!fits_fread_matrix (missing[0].f, &matrix, &missing[0].header)) { + if (VERBOSE) fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + /* read Missing table header */ + if (!fits_fread_header (missing[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table missing header"); + return (FALSE); + } + /* read Missing table data */ + if (!fits_fread_ftable_data (missing[0].f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table missing data"); + return (FALSE); + } + /* no conversions currently defined */ + catalog[0].missing = fits_table_get_Missing (&ftable, &catalog[0].Nmissing, NULL); + if (Nmissing != catalog[0].Nmissing) { + fprintf (stderr, "Warning: mismatch between Nmissing in PHU and Table headers (%d vs %d)\n", Nmissing, catalog[0].Nmissing); + } + fits_free_header (&header); + fits_free_matrix (&matrix); + } + catalog[0].missing_catalog = missing; + + /*** Secfilt Table ***/ + + secfilt = NULL; + if (catalog[0].catflags & LOAD_SECF) { + ALLOCATE (secfilt, Catalog, 1); + + /* get split filename from main header (paths relative to cpt file) */ + if (!fits_scan (&catalog[0].header, "SECFILT", "%s", 1, string)) return (FALSE); + ALLOCATE (secfilt[0].filename, char, strlen(path) + strlen(string) + 2); + sprintf (secfilt[0].filename, "%s/%s", path, string); + + /* lock & open catalog file */ + if (lock_catalog (secfilt, catalog[0].lockmode) != 1) { + fprintf (stderr, "cannot access secfilt file %s\n", secfilt[0].filename); + exit (2); + } + + /* read PHU */ + if (!fits_load_header (secfilt[0].f, &secfilt[0].header)) { + if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", secfilt[0].filename); + return (FALSE); + } + /* matrix should be empty */ + if (!fits_fread_matrix (secfilt[0].f, &matrix, &secfilt[0].header)) { + if (VERBOSE) fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + /* read secfilt table header */ + if (!fits_fread_header (secfilt[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table secfilt header"); + return (FALSE); + } + /* read secfilt table data */ + if (!fits_fread_ftable_data (secfilt[0].f, &ftable)) { + if (VERBOSE) fprintf (stderr, "can't read table secfilt data"); + return (FALSE); + } + Nexpect = catalog[0].Nsecfilt * catalog[0].Naverage; + catalog[0].secfilt = FtableToSecFilt (&ftable, &Nitems, &catalog[0].catformat); + if (Nexpect != Nitems) { + fprintf (stderr, "Warning: mismatch between Nsecfilt items in PHU and Table headers (%d vs %d)\n", Nexpect, Nitems); + } + fits_free_header (&header); + fits_free_matrix (&matrix); + } + catalog[0].secfilt_catalog = secfilt; + + /* save the current number so we can do partial updates */ + catalog[0].Nave_disk = Naverage; + catalog[0].Nmeas_disk = Nmeasure; + catalog[0].Nmiss_disk = Nmissing; + + return (TRUE); +} + +/* save_catalog_split writes complete new files from scratch */ + +int save_catalog_split (Catalog *catalog, char VERBOSE) { + + int Nitems; + Matrix matrix; + Header header; + FTable ftable; + Catalog *measure, *missing, *secfilt; + + ftable.header = &header; + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (TRUE); + } + + /* make sure header is consistent with data */ + fits_modify (&catalog[0].header, "NSTARS", "%d", 1, catalog[0].Naverage); + fits_modify (&catalog[0].header, "NMEAS", "%d", 1, catalog[0].Nmeasure); + fits_modify (&catalog[0].header, "NMISS", "%d", 1, catalog[0].Nmissing); + fits_modify (&catalog[0].header, "NSECFILT", "%d", 1, catalog[0].Nsecfilt); + fits_modify (&catalog[0].header, "EXTEND", "%t", 1, TRUE); + + /* rewind file pointers and truncate (file is still open) */ + fseek (catalog[0].f, 0, SEEK_SET); + + /* write table PHU header - always write this out */ + /* XXX EAM : check if disk file size has changed */ + if (!fits_fwrite_header (catalog[0].f, &catalog[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* in split mode, we can save only part of the data */ + + /*** Average Table ***/ + + if (catalog[0].average != NULL) { + ftruncate (fileno (catalog[0].f), catalog[0].header.size); + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&catalog[0].header, &matrix); + if (!fits_fwrite_matrix (catalog[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out Average table (convert to FITS table format) */ + AverageToFtable (&ftable, catalog[0].average, catalog[0].Naverage, catalog[0].catformat); + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Measure Table ***/ + + if (catalog[0].measure != NULL) { + + measure = catalog[0].measure_catalog; + + /* XXX EAM : warn about this condition; add code to handle? */ + if (catalog[0].Nmeas_off != 0) { + fprintf (stderr, "WARNING: LOAD_MEAS_META mixed with save??\n"); + fprintf (stderr, "WARNING: this should not be allowed to happen!\n"); + } + + /* rewind file pointers and truncate (file is still open) */ + fseek (measure[0].f, 0, SEEK_SET); + ftruncate (fileno (measure[0].f), 0); + + /* write table PHU header */ + if (!fits_fwrite_header (measure[0].f, &measure[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&measure[0].header, &matrix); + if (!fits_fwrite_matrix (measure[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out Measure table (convert to FITS table format) */ + MeasureToFtable (&ftable, catalog[0].measure, catalog[0].Nmeasure, catalog[0].catformat); + if (!fits_fwrite_Theader (measure[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (measure[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Missing Table ***/ + + if (catalog[0].missing != NULL) { + + missing = catalog[0].missing_catalog; + + /* rewind file pointers and truncate (file is still open) */ + fseek (missing[0].f, 0, SEEK_SET); + ftruncate (fileno (missing[0].f), 0); + + /* write table PHU header */ + if (!fits_fwrite_header (missing[0].f, &missing[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&missing[0].header, &matrix); + if (!fits_fwrite_matrix (missing[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out Missing table (convert to FITS table format) */ + fits_table_set_Missing (&ftable, catalog[0].missing, catalog[0].Nmissing); + if (!fits_fwrite_Theader (missing[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (missing[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Secfilt Table ***/ + + if (catalog[0].secfilt != NULL) { + + secfilt = catalog[0].secfilt_catalog; + + /* rewind file pointers and truncate (file is still open) */ + fseek (secfilt[0].f, 0, SEEK_SET); + ftruncate (fileno (secfilt[0].f), 0); + + /* write table PHU header */ + if (!fits_fwrite_header (secfilt[0].f, &secfilt[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&secfilt[0].header, &matrix); + if (!fits_fwrite_matrix (secfilt[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out SecFilt table (convert to FITS table format) */ + Nitems = catalog[0].Naverage * catalog[0].Nsecfilt; + SecFiltToFtable (&ftable, catalog[0].secfilt, Nitems, catalog[0].catformat); + if (!fits_fwrite_Theader (secfilt[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (secfilt[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + } + + return (TRUE); +} + +/* update_catalog_split only writes new lines to file. + * if file is empty, call save_catalog_split instead. + * XXX EAM : save_catalog SHOULD do this + */ + +int update_catalog_split (Catalog *catalog, char VERBOSE) { + + int i, Nx, Ny, Nlines; + int Nitems, Nskip, Nout, Ndisk, Nstart; + Matrix matrix; + Header header; + FTable ftable; + VTable vtable; + Catalog *measure, *missing, *secfilt; + + ftable.header = &header; + vtable.header = &header; + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (TRUE); + } + + /* make sure header is consistent with data */ + fits_modify (&catalog[0].header, "NSTARS", "%d", 1, catalog[0].Naverage); + fits_modify (&catalog[0].header, "NMEAS", "%d", 1, catalog[0].Nmeasure + catalog[0].Nmeas_off); + fits_modify (&catalog[0].header, "NMISS", "%d", 1, catalog[0].Nmissing); + fits_modify (&catalog[0].header, "NSECFILT", "%d", 1, catalog[0].Nsecfilt); + fits_modify (&catalog[0].header, "EXTEND", "%t", 1, TRUE); + + /* rewind file pointers and truncate (file is still open) */ + fseek (catalog[0].f, 0, SEEK_SET); + + /* write table PHU header - always write this out */ + /* XXX EAM : check if disk file size has changed */ + if (!fits_fwrite_header (catalog[0].f, &catalog[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /*** Average Table ***/ + + if (catalog[0].average != NULL) { + + /* skip past matrix (already at end of header) */ + Nskip = fits_matrix_size (&catalog[0].header); + fseek (catalog[0].f, Nskip, SEEK_CUR); + + /* how many lines to write out? */ + Nout = catalog[0].Naverage - catalog[0].Nave_disk; + + /* write out Average table (convert to FITS table format) */ + AverageToFtable (&ftable, catalog[0].average, catalog[0].Naverage, catalog[0].catformat); + /* convert only output rows to vtable */ + fits_table_to_vtable (&ftable, &vtable, catalog[0].Nave_disk, Nout); + + if (!fits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_vtable (catalog[0].f, &vtable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_vtable (&vtable); + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Measure Table ***/ + + if (catalog[0].measure != NULL) { + + measure = catalog[0].measure_catalog; + + /* skip past PHU header and matrix */ + Nskip = measure[0].header.size + fits_matrix_size (&measure[0].header); + fseek (measure[0].f, Nskip, SEEK_SET); + + Ndisk = catalog[0].Nmeas_disk; + Nstart = catalog[0].Nmeas_disk - catalog[0].Nmeas_off; /* where is first new line? */ + Nout = catalog[0].Nmeasure - Nstart; /* how many lines to write out? */ + Nlines = catalog[0].Nmeasure + catalog[0].Nmeas_off; /* how many lines total in file */ + + /* convert to output format FITS table (only rows for output : 0 - Nout) */ + MeasureToFtable (&ftable, &catalog[0].measure[Nstart], Nout, catalog[0].catformat); + + fits_scan (&header, "NAXIS1", "%d", 1, &Nx); + fits_scan (&header, "NAXIS2", "%d", 1, &Ny); + + /* convert all output rows to vtable */ + ALLOCATE (vtable.row, int, MAX (1, Nout)); + ALLOCATE (vtable.buffer, char *, MAX (1, Nout)); + for (i = 0; i < Nout; i++) { + ALLOCATE (vtable.buffer[i], char, MAX (1, Nx)); + memcpy (vtable.buffer[i], &ftable.buffer[i*Nx], Nx); + vtable.row[i] = i + Ndisk; + } + + /* modify vtable to represent full disk table */ + fits_modify (&header, "NAXIS2", "%d", 1, Nlines); + header.Naxis[1] = Nlines; + + vtable.size = fits_matrix_size (&header); + vtable.Nrow = Nout; + vtable.pad = vtable.size - Nx*Ny; + + if (!fits_fwrite_Theader (measure[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_vtable (measure[0].f, &vtable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_vtable (&vtable); + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Missing Table ***/ + /* missing table CANNOT be written unsorted, thus it is always written + out in full */ + + if (catalog[0].missing != NULL) { + + missing = catalog[0].missing_catalog; + + /* rewind file pointers and truncate (file is still open) */ + fseek (missing[0].f, 0, SEEK_SET); + ftruncate (fileno (missing[0].f), 0); + + /* write table PHU header */ + if (!fits_fwrite_header (missing[0].f, &missing[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + + /* this is probably a NOP, do I have to keep it in? */ + fits_create_matrix (&missing[0].header, &matrix); + if (!fits_fwrite_matrix (missing[0].f, &matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + fits_free_matrix (&matrix); + + /* write out Missing table (convert to FITS table format) */ + fits_table_set_Missing (&ftable, catalog[0].missing, catalog[0].Nmissing); + if (!fits_fwrite_Theader (missing[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (missing[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_table (&ftable); + fits_free_header (&header); + } + + /*** Secfilt Table ***/ + + if (catalog[0].secfilt != NULL) { + + secfilt = catalog[0].secfilt_catalog; + + /* skip past PHU header and matrix */ + Nskip = secfilt[0].header.size + fits_matrix_size (&secfilt[0].header); + fseek (secfilt[0].f, Nskip, SEEK_SET); + + /* how many lines to write out? */ + Nout = catalog[0].Nsecfilt * (catalog[0].Naverage - catalog[0].Nave_disk); + Ndisk = catalog[0].Nsecfilt * catalog[0].Nave_disk; + + /* convert to output format FITS table */ + Nitems = catalog[0].Naverage * catalog[0].Nsecfilt; + SecFiltToFtable (&ftable, catalog[0].secfilt, Nitems, catalog[0].catformat); + /* convert only output rows to vtable */ + fits_table_to_vtable (&ftable, &vtable, Ndisk, Nout); + + if (!fits_fwrite_Theader (secfilt[0].f, &header)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_vtable (secfilt[0].f, &vtable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + fits_free_vtable (&vtable); + fits_free_table (&ftable); + fits_free_header (&header); + } + + return (TRUE); +} + +/* in split mode, extra files are linked to catalog->measure_catalog, etc. Each + has a valid filename, f, header. The primary catalog data pointers are set + to point at the corresponding entries from the elements on the chain. An + unloaded entry has a null pointer here. +*/ + +/* XXX EAM : this file needs work on the error exit conditions and memory leaks, esp under errors */ + +/* XXX EAM : update is not efficient. MeasureToFtable should only + convert the new rows (Nmeas_disk to Nmeasure). the resulting + table represents the end rows of the vtable. we need to define + the vtable based on the ftable, but with Ny = Nmeasure */ + Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert.c (revision 5388) @@ -0,0 +1,438 @@ +# include + +/* The Ftable-TYPE conversion functions determine the format of table based on EXTNAME in header. + they convert the table to the internal format, and set 'format'. + + The TYPE-Ftable conversions functions create output tables in the format requested + by the 'format' function parameter. + */ + +/** this file might be more readable if I use macros for the repetative + constructions below **/ + +/*** Average / FTable conversion functions ***/ + +Average *FtableToAverage (FTable *ftable, int *Naverage, int *format) { + + Average *average; + char extname[80]; + + /* convert to the internal format */ + if (!fits_scan (ftable[0].header, "EXTNAME", "%s", 1, extname)) { + fprintf (stderr, "EXTNAME missing for average table\n"); + return (FALSE); + } + if (!strcmp (extname, "DVO_AVERAGE")) { + average = fits_table_get_Average (ftable, Naverage, NULL); + *format = DVO_FORMAT_INTERNAL; + return (average); + } + if (!strcmp (extname, "DVO_AVERAGE_ELIXIR")) { + AverageElixir *tmpAverage; + tmpAverage = fits_table_get_AverageElixir (ftable, Naverage, NULL); + average = AverageElixirToInternal (tmpAverage, *Naverage); + free (tmpAverage); + *format = DVO_FORMAT_ELIXIR; + return (average); + } + if (!strcmp (extname, "DVO_AVERAGE_LONEOS")) { + AverageLoneos *tmpAverage; + tmpAverage = fits_table_get_AverageLoneos (ftable, Naverage, NULL); + average = AverageLoneosToInternal (tmpAverage, *Naverage); + free (tmpAverage); + *format = DVO_FORMAT_LONEOS; + return (average); + } + if (!strcmp (extname, "DVO_AVERAGE_PANSTARRS")) { + AveragePanstarrs *tmpAverage; + tmpAverage = fits_table_get_AveragePanstarrs (ftable, Naverage, NULL); + average = AveragePanstarrsToInternal (tmpAverage, *Naverage); + free (tmpAverage); + *format = DVO_FORMAT_PANSTARRS; + return (average); + } + + fprintf (stderr, "table format unknown: %s\n", extname); + + *Naverage = 0; + return (NULL); +} + +int AverageToFtable (FTable *ftable, Average *average, int Naverage, int format) { + + AverageElixir *tmpAverageElixir; + AverageLoneos *tmpAverageLoneos ; + AveragePanstarrs *tmpAveragePanstarrs; + + /* convert from the internal format */ + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_table_set_Average (ftable, average, Naverage); + break; + case DVO_FORMAT_ELIXIR: + tmpAverageElixir = AverageInternalToElixir (average, Naverage); + fits_table_set_AverageElixir (ftable, tmpAverageElixir, Naverage); + free (tmpAverageElixir); + break; + case DVO_FORMAT_LONEOS: + tmpAverageLoneos = AverageInternalToLoneos (average, Naverage); + fits_table_set_AverageLoneos (ftable, tmpAverageLoneos , Naverage); + free (tmpAverageLoneos ); + break; + case DVO_FORMAT_PANSTARRS: + tmpAveragePanstarrs = AverageInternalToPanstarrs (average, Naverage); + fits_table_set_AveragePanstarrs (ftable, tmpAveragePanstarrs, Naverage); + free (tmpAveragePanstarrs); + break; + default: + fprintf (stderr, "table format unknown (average)\n"); + return (FALSE); + } + return (TRUE); +} + +/*** Measure / FTable conversion functions ***/ + +Measure *FtableToMeasure (FTable *ftable, int *Nmeasure, int *format) { + + Measure *measure; + char extname[80]; + + /* convert to the internal format */ + if (!fits_scan (ftable[0].header, "EXTNAME", "%s", 1, extname)) { + fprintf (stderr, "EXTNAME missing for measure table\n"); + return (FALSE); + } + if (!strcmp (extname, "DVO_MEASURE")) { + measure = fits_table_get_Measure (ftable, Nmeasure, NULL); + *format = DVO_FORMAT_INTERNAL; + return (measure); + } + if (!strcmp (extname, "DVO_MEASURE_ELIXIR")) { + MeasureElixir *tmpMeasure; + tmpMeasure = fits_table_get_MeasureElixir (ftable, Nmeasure, NULL); + measure = MeasureElixirToInternal (tmpMeasure, *Nmeasure); + free (tmpMeasure); + *format = DVO_FORMAT_ELIXIR; + return (measure); + } + if (!strcmp (extname, "DVO_MEASURE_LONEOS")) { + MeasureLoneos *tmpMeasure; + tmpMeasure = fits_table_get_MeasureLoneos (ftable, Nmeasure, NULL); + measure = MeasureLoneosToInternal (tmpMeasure, *Nmeasure); + free (tmpMeasure); + *format = DVO_FORMAT_LONEOS; + return (measure); + } + if (!strcmp (extname, "DVO_MEASURE_PANSTARRS")) { + MeasurePanstarrs *tmpMeasure; + tmpMeasure = fits_table_get_MeasurePanstarrs (ftable, Nmeasure, NULL); + measure = MeasurePanstarrsToInternal (tmpMeasure, *Nmeasure); + free (tmpMeasure); + *format = DVO_FORMAT_PANSTARRS; + return (measure); + } + + fprintf (stderr, "table format unknown: %s\n", extname); + + *Nmeasure = 0; + return (NULL); +} + +int MeasureToFtable (FTable *ftable, Measure *measure, int Nmeasure, int format) { + + MeasureElixir *tmpMeasureElixir; + MeasureLoneos *tmpMeasureLoneos; + MeasurePanstarrs *tmpMeasurePanstarrs; + + /* convert from the internal format */ + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_table_set_Measure (ftable, measure, Nmeasure); + break; + case DVO_FORMAT_ELIXIR: + tmpMeasureElixir = MeasureInternalToElixir (measure, Nmeasure); + fits_table_set_MeasureElixir (ftable, tmpMeasureElixir, Nmeasure); + free (tmpMeasureElixir); + break; + case DVO_FORMAT_LONEOS: + tmpMeasureLoneos = MeasureInternalToLoneos (measure, Nmeasure); + fits_table_set_MeasureLoneos (ftable, tmpMeasureLoneos, Nmeasure); + free (tmpMeasureLoneos); + break; + case DVO_FORMAT_PANSTARRS: + tmpMeasurePanstarrs = MeasureInternalToPanstarrs (measure, Nmeasure); + fits_table_set_MeasurePanstarrs (ftable, tmpMeasurePanstarrs, Nmeasure); + free (tmpMeasurePanstarrs); + break; + default: + fprintf (stderr, "table format unknown (measure)\n"); + return (FALSE); + } + return (TRUE); +} + +/*** SecFilt / FTable conversion functions ***/ + +SecFilt *FtableToSecFilt (FTable *ftable, int *Nsecfilt, int *format) { + + SecFilt *secfilt; + char extname[80]; + + /* convert to the internal format */ + if (!fits_scan (ftable[0].header, "EXTNAME", "%s", 1, extname)) { + fprintf (stderr, "EXTNAME missing for secfilt table\n"); + return (FALSE); + } + if (!strcmp (extname, "DVO_SECFILT")) { + secfilt = fits_table_get_SecFilt (ftable, Nsecfilt, NULL); + *format = DVO_FORMAT_INTERNAL; + return (secfilt); + } + if (!strcmp (extname, "DVO_SECFILT_ELIXIR")) { + SecFiltElixir *tmpSecFilt; + tmpSecFilt = fits_table_get_SecFiltElixir (ftable, Nsecfilt, NULL); + secfilt = SecFiltElixirToInternal (tmpSecFilt, *Nsecfilt); + free (tmpSecFilt); + *format = DVO_FORMAT_ELIXIR; + return (secfilt); + } + if (!strcmp (extname, "DVO_SECFILT_LONEOS")) { + SecFiltLoneos *tmpSecFilt; + tmpSecFilt = fits_table_get_SecFiltLoneos (ftable, Nsecfilt, NULL); + secfilt = SecFiltLoneosToInternal (tmpSecFilt, *Nsecfilt); + free (tmpSecFilt); + *format = DVO_FORMAT_LONEOS; + return (secfilt); + } + if (!strcmp (extname, "DVO_SECFILT_PANSTARRS")) { + SecFiltPanstarrs *tmpSecFilt; + tmpSecFilt = fits_table_get_SecFiltPanstarrs (ftable, Nsecfilt, NULL); + secfilt = SecFiltPanstarrsToInternal (tmpSecFilt, *Nsecfilt); + free (tmpSecFilt); + *format = DVO_FORMAT_PANSTARRS; + return (secfilt); + } + + fprintf (stderr, "table format unknown: %s\n", extname); + + *Nsecfilt = 0; + return (NULL); +} + +int SecFiltToFtable (FTable *ftable, SecFilt *secfilt, int Nsecfilt, int format) { + + SecFiltElixir *tmpSecFiltElixir; + SecFiltLoneos *tmpSecFiltLoneos; + SecFiltPanstarrs *tmpSecFiltPanstarrs; + + /* convert from the internal format */ + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_table_set_SecFilt (ftable, secfilt, Nsecfilt); + break; + case DVO_FORMAT_ELIXIR: + tmpSecFiltElixir = SecFiltInternalToElixir (secfilt, Nsecfilt); + fits_table_set_SecFiltElixir (ftable, tmpSecFiltElixir, Nsecfilt); + free (tmpSecFiltElixir); + break; + case DVO_FORMAT_LONEOS: + tmpSecFiltLoneos = SecFiltInternalToLoneos (secfilt, Nsecfilt); + fits_table_set_SecFiltLoneos (ftable, tmpSecFiltLoneos, Nsecfilt); + free (tmpSecFiltLoneos); + break; + case DVO_FORMAT_PANSTARRS: + tmpSecFiltPanstarrs = SecFiltInternalToPanstarrs (secfilt, Nsecfilt); + fits_table_set_SecFiltPanstarrs (ftable, tmpSecFiltPanstarrs, Nsecfilt); + free (tmpSecFiltPanstarrs); + break; + default: + fprintf (stderr, "table format unknown (secfilt)\n"); + return (FALSE); + } + return (TRUE); +} + +/*** Image Conversions ***/ + +int FtableToImage (FTable *ftable, Header *theader, int *format) { + + int Nimage; + char extname[80]; + + /* extname may be set from outside if the source is RAW not MEF */ + if (*format == DVO_FORMAT_ELIXIR) { + ImageElixir *tmpimage; + tmpimage = fits_table_get_ImageElixir (ftable, &Nimage, NULL); + ftable[0].buffer = (char *) ImageElixirToInternal (tmpimage, Nimage); + free (tmpimage); + goto set_header; + } + + /* convert to the internal format */ + if (!fits_scan (ftable[0].header, "EXTNAME", "%s", 1, extname)) { + fprintf (stderr, "EXTNAME missing for image table\n"); + return (FALSE); + } + if (!strcmp (extname, "DVO_IMAGE")) { + Image *image; + image = fits_table_get_Image (ftable, &Nimage, NULL); + *format = DVO_FORMAT_INTERNAL; + return (TRUE); + } + if (!strcmp (extname, "DVO_IMAGE_ELIXIR")) { + ImageElixir *tmpimage; + *format = DVO_FORMAT_ELIXIR; + tmpimage = fits_table_get_ImageElixir (ftable, &Nimage, NULL); + ftable[0].buffer = (char *) ImageElixirToInternal (tmpimage, Nimage); + free (tmpimage); + goto set_header; + } + if (!strcmp (extname, "DVO_IMAGE_LONEOS")) { + ImageLoneos *tmpimage; + *format = DVO_FORMAT_LONEOS; + tmpimage = fits_table_get_ImageLoneos (ftable, &Nimage, NULL); + ftable[0].buffer = (char *) ImageLoneosToInternal (tmpimage, Nimage); + free (tmpimage); + goto set_header; + } + if (!strcmp (extname, "DVO_IMAGE_PANSTARRS")) { + ImagePanstarrs *tmpimage; + *format = DVO_FORMAT_PANSTARRS; + tmpimage = fits_table_get_ImagePanstarrs (ftable, &Nimage, NULL); + ftable[0].buffer = (char *) ImagePanstarrsToInternal (tmpimage, Nimage); + free (tmpimage); + goto set_header; + } + fprintf (stderr, "table format unknown: %s\n", extname); + return (FALSE); + +set_header: + fits_free_header (theader); + fits_table_mkheader_Image (theader); + fits_modify (theader, "NAXIS2", "%d", 1, Nimage); + theader[0].Naxis[1] = Nimage; + ftable[0].size = fits_matrix_size (theader); + return (TRUE); +} + +int ImageToFtable (FTable *ftable, Header *theader, int format) { + + int Nimage; + ImageElixir *tmpImageElixir; + ImageLoneos *tmpImageLoneos; + ImagePanstarrs *tmpImagePanstarrs; + + Nimage = theader[0].Naxis[1]; + + /* convert from the internal format */ + switch (format) { + case DVO_FORMAT_INTERNAL: + fits_convert_Image ((Image *) ftable[0].buffer, sizeof(Image), Nimage); + break; + case DVO_FORMAT_ELIXIR: + tmpImageElixir = ImageInternalToElixir ((Image *) ftable[0].buffer, Nimage); + free (ftable[0].buffer); + fits_table_set_ImageElixir (ftable, tmpImageElixir, Nimage); + free (tmpImageElixir); + break; + case DVO_FORMAT_LONEOS: + tmpImageLoneos = ImageInternalToLoneos ((Image *) ftable[0].buffer, Nimage); + free (ftable[0].buffer); + fits_table_set_ImageLoneos (ftable, tmpImageLoneos, Nimage); + free (tmpImageLoneos); + break; + case DVO_FORMAT_PANSTARRS: + tmpImagePanstarrs = ImageInternalToPanstarrs ((Image *) ftable[0].buffer, Nimage); + free (ftable[0].buffer); + fits_table_set_ImagePanstarrs (ftable, tmpImagePanstarrs, Nimage); + free (tmpImagePanstarrs); + break; + default: + fprintf (stderr, "table format unknown (image ftable)\n"); + return (FALSE); + } + return (TRUE); +} + +int ImageToVtable (VTable *vtable, Header *theader, int format) { + + int i, Nrow, Nimage; + + Nrow = vtable[0].Nrow; + + /* convert from the internal format */ + if (format == DVO_FORMAT_INTERNAL) { + for (i = 0; i < Nrow; i++) { + fits_convert_Image ((Image *) vtable[0].buffer[i], sizeof(Image), 1); + } + return (TRUE); + } + if (format == DVO_FORMAT_ELIXIR) { + ImageElixir *tmpImage; + + /* convert table rows from internal to external format */ + for (i = 0; i < Nrow; i++) { + tmpImage = ImageInternalToElixir ((Image *) vtable[0].buffer[i], 1); + fits_convert_ImageElixir (tmpImage, sizeof(ImageElixir), 1); + free (vtable[0].buffer[i]); + vtable[0].buffer[i] = (char *) tmpImage; + } + + /* convert header from old format to new format */ + fits_scan (theader, "NAXIS2", "%d", 1, &Nimage); + + fits_free_header (theader); + fits_table_mkheader_ImageElixir (theader); + + fits_modify (theader, "NAXIS2", "%d", 1, Nimage); + theader[0].Naxis[1] = Nimage; + vtable[0].size = fits_matrix_size (theader); + return (TRUE); + } + if (format == DVO_FORMAT_LONEOS) { + ImageLoneos *tmpImage; + + /* convert table rows from internal to external format */ + for (i = 0; i < Nrow; i++) { + tmpImage = ImageInternalToLoneos ((Image *) vtable[0].buffer[i], 1); + fits_convert_ImageLoneos (tmpImage, sizeof(ImageLoneos), 1); + free (vtable[0].buffer[i]); + vtable[0].buffer[i] = (char *) tmpImage; + } + + /* convert header from old format to new format */ + fits_scan (theader, "NAXIS2", "%d", 1, &Nimage); + + fits_free_header (theader); + fits_table_mkheader_ImageLoneos (theader); + + fits_modify (theader, "NAXIS2", "%d", 1, Nimage); + theader[0].Naxis[1] = Nimage; + vtable[0].size = fits_matrix_size (theader); + return (TRUE); + } + if (format == DVO_FORMAT_PANSTARRS) { + ImagePanstarrs *tmpImage; + + /* convert table rows from internal to external format */ + for (i = 0; i < Nrow; i++) { + tmpImage = ImageInternalToPanstarrs ((Image *) vtable[0].buffer[i], 1); + fits_convert_ImagePanstarrs (tmpImage, sizeof(ImagePanstarrs), 1); + free (vtable[0].buffer[i]); + vtable[0].buffer[i] = (char *) tmpImage; + } + + /* convert header from old format to new format */ + fits_scan (theader, "NAXIS2", "%d", 1, &Nimage); + fits_free_header (theader); + fits_table_mkheader_ImagePanstarrs (theader); + fits_modify (theader, "NAXIS2", "%d", 1, Nimage); + theader[0].Naxis[1] = Nimage; + vtable[0].size = fits_matrix_size (theader); + return (TRUE); + } + fprintf (stderr, "table format unknown (image vtable)\n"); + return (FALSE); +} + Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_elixir.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 5388) @@ -0,0 +1,245 @@ +# include + +/* convert elixir-format measures to internal measures */ +Measure *MeasureElixirToInternal (MeasureElixir *in, int Nvalues) { + + int i; + Measure *out; + + ALLOCATE (out, Measure, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR_PS = in[i].dR * 0.01; + out[i].dD_PS = in[i].dD * 0.01; + out[i].M_PS = in[i].M * 0.001; + out[i].dM_PS = in[i].dM * 0.001; + out[i].dt_PS = in[i].dt * 0.001; + out[i].Mcal_PS = in[i].Mcal * 0.001; + out[i].Mgal_PS = in[i].Mgal * 0.001; + out[i].airmass_PS = in[i].airmass * 0.001; + out[i].FWx = in[i].FWx; + out[i].FWy = in[i].fwy * in[i].FWx * 0.01; + out[i].theta = in[i].theta; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].flags = in[i].flags; + } + return (out); +} + +/* convert internal measures to elixir-format measures */ +MeasureElixir *MeasureInternalToElixir (Measure *in, int Nvalues) { + + int i; + MeasureElixir *out; + + ALLOCATE (out, MeasureElixir, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR = in[i].dR_PS * 100.0; + out[i].dD = in[i].dD_PS * 100.0; + out[i].M = in[i].M_PS * 1000.0; + out[i].dM = in[i].dM_PS * 1000.0; + out[i].dt = in[i].dt_PS * 1000.0; + out[i].Mcal = in[i].Mcal_PS * 1000.0; + out[i].Mgal = in[i].Mgal_PS * 1000.0; + out[i].airmass = in[i].airmass_PS * 1000.0; + out[i].FWx = in[i].FWx; + out[i].fwy = 100.0 * in[i].FWy / in[i].FWx; + out[i].theta = in[i].theta; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].flags = in[i].flags; + } + return (out); +} + +/* convert elixir-format averages to internal averages */ +Average *AverageElixirToInternal (AverageElixir *in, int Nvalues) { + + int i; + Average *out; + + ALLOCATE (out, Average, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R_PS = in[i].R; + out[i].D_PS = in[i].D; + out[i].M_PS = in[i].M * 0.001; + out[i].dM_PS = in[i].dM * 0.001; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + out[i].Xg = in[i].Xg; + } + return (out); +} + +/* convert internal averages to elixir-format averages */ +AverageElixir *AverageInternalToElixir (Average *in, int Nvalues) { + + int i; + AverageElixir *out; + + ALLOCATE (out, AverageElixir, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R_PS; + out[i].D = in[i].D_PS; + out[i].M = in[i].M_PS * 1000.0; + out[i].dM = in[i].dM_PS * 1000.0; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + out[i].Xg = in[i].Xg; + } + return (out); +} + +/* convert elixir-format secfilts to internal secfilts */ +SecFilt *SecFiltElixirToInternal (SecFiltElixir *in, int Nvalues) { + + int i; + SecFilt *out; + + ALLOCATE (out, SecFilt, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M_PS = in[i].M * 0.001; + out[i].dM_PS = in[i].dM * 0.001; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert internal secfilts to elixir-format secfilts */ +SecFiltElixir *SecFiltInternalToElixir (SecFilt *in, int Nvalues) { + + int i; + SecFiltElixir *out; + + ALLOCATE (out, SecFiltElixir, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M = in[i].M_PS * 1000.0; + out[i].dM = in[i].dM_PS * 1000.0; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert elixir-format images to internal images */ +Image *ImageElixirToInternal (ImageElixir *in, int Nvalues) { + + int i; + Image *out; + + ALLOCATE (out, Image, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz_PS = in[i].secz * 0.001; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit_PS = in[i].apmifit * 0.001; + out[i].dapmifit_PS = in[i].dapmifit * 0.001; + out[i].Mcal_PS = in[i].Mcal * 0.001; + out[i].dMcal_PS = in[i].dMcal * 0.001; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} + +/* convert internal images to elixir-format images */ +ImageElixir *ImageInternalToElixir (Image *in, int Nvalues) { + + int i; + ImageElixir *out; + + ALLOCATE (out, ImageElixir, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz_PS * 1000.0; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit_PS * 1000.0; + out[i].dapmifit = in[i].dapmifit_PS * 1000.0; + out[i].Mcal = in[i].Mcal_PS * 1000.0; + out[i].dMcal = in[i].dMcal_PS * 1000.0; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_loneos.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 5388) @@ -0,0 +1,243 @@ +# include + +/* convert loneos-format measures to internal measures */ +Measure *MeasureLoneosToInternal (MeasureLoneos *in, int Nvalues) { + + int i; + Measure *out; + + ALLOCATE (out, Measure, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR_PS = in[i].dR * 0.01; + out[i].dD_PS = in[i].dD * 0.01; + out[i].M_PS = in[i].M * 0.001; + out[i].dM_PS = in[i].dM * 0.001; + out[i].Mcal_PS = in[i].Mcal * 0.001; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + + /* flags and averef are now split */ + out[i].averef = in[i].averef & 0x00ffffff; + out[i].flags = (in[i].averef & 0xff000000) >> 24; + + /* these values don't exist in the Loneos format */ + out[i].Mgal_PS = in[i].M; + out[i].dt_PS = 0; + out[i].airmass_PS = 0; + out[i].FWx = 0; + out[i].FWy = 0; + out[i].theta = 0; + } + return (out); +} + +/* convert internal measures to loneos-format measures */ +MeasureLoneos *MeasureInternalToLoneos (Measure *in, int Nvalues) { + + int i; + MeasureLoneos *out; + + ALLOCATE (out, MeasureLoneos, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR = in[i].dR_PS * 100.0; + out[i].dD = in[i].dD_PS * 100.0; + out[i].M = in[i].M_PS * 1000.0; + out[i].dM = in[i].dM_PS * 1000.0; + out[i].Mcal = in[i].Mcal_PS * 1000.0; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + + /* flags and averef are merged in Loneos */ + out[i].averef = (in[i].averef & 0x00ffffff) | (in[i].flags << 24); + } + return (out); +} + +/* convert loneos-format averages to internal averages */ +Average *AverageLoneosToInternal (AverageLoneos *in, int Nvalues) { + + int i; + Average *out; + + ALLOCATE (out, Average, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R_PS = in[i].R; + out[i].D_PS = in[i].D; + out[i].M_PS = in[i].M * 0.001; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + + /* these don't exist in Loneos */ + out[i].dM_PS = 0xffff; + out[i].Xg = 0xffff; + } + return (out); +} + +/* convert internal averages to loneos-format averages */ +AverageLoneos *AverageInternalToLoneos (Average *in, int Nvalues) { + + int i; + AverageLoneos *out; + + ALLOCATE (out, AverageLoneos, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R_PS; + out[i].D = in[i].D_PS; + out[i].M = in[i].M_PS * 1000.0; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + } + return (out); +} + +/* convert loneos-format secfilts to internal secfilts */ +SecFilt *SecFiltLoneosToInternal (SecFiltLoneos *in, int Nvalues) { + + int i; + SecFilt *out; + + ALLOCATE (out, SecFilt, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M_PS = in[i].M*0.001; + out[i].dM_PS = 0xffff; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert internal secfilts to loneos-format secfilts */ +SecFiltLoneos *SecFiltInternalToLoneos (SecFilt *in, int Nvalues) { + + int i; + SecFiltLoneos *out; + + ALLOCATE (out, SecFiltLoneos, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M = in[i].M_PS*1000.0; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert loneos-format images to internal images */ +Image *ImageLoneosToInternal (ImageLoneos *in, int Nvalues) { + + int i; + Image *out; + + ALLOCATE (out, Image, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz_PS = in[i].secz * 0.001; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit_PS = in[i].apmifit * 0.001; + out[i].dapmifit_PS = in[i].dapmifit * 0.001; + out[i].Mcal_PS = in[i].Mcal * 0.001; + out[i].dMcal_PS = in[i].dMcal * 0.001; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} + +/* convert internal images to loneos-format images */ +ImageLoneos *ImageInternalToLoneos (Image *in, int Nvalues) { + + int i; + ImageLoneos *out; + + ALLOCATE (out, ImageLoneos, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz_PS * 1000.0; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit_PS * 1000.0; + out[i].dapmifit = in[i].dapmifit_PS * 1000.0; + out[i].Mcal = in[i].Mcal_PS * 1000.0; + out[i].dMcal = in[i].dMcal_PS * 1000.0; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_panstarrs.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_panstarrs.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_convert_panstarrs.c (revision 5388) @@ -0,0 +1,247 @@ +# include + +/*** note that these structures are identical ***/ + +/* convert panstarrs-format measures to internal measures */ +Measure *MeasurePanstarrsToInternal (MeasurePanstarrs *in, int Nvalues) { + + int i; + Measure *out; + + ALLOCATE (out, Measure, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR_PS = in[i].dR; + out[i].dD_PS = in[i].dD; + out[i].M_PS = in[i].M; + out[i].dM_PS = in[i].dM; + out[i].Mcal_PS = in[i].Mcal; + out[i].Mgal_PS = in[i].Mgal; + out[i].airmass_PS = in[i].airmass; + out[i].dt_PS = in[i].dt; + out[i].FWx = in[i].FWx; + out[i].FWy = in[i].FWy; + out[i].theta = in[i].theta; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].flags = in[i].flags; + } + return (out); +} + +/* convert internal measures to panstarrs-format measures */ +MeasurePanstarrs *MeasureInternalToPanstarrs (Measure *in, int Nvalues) { + + int i; + MeasurePanstarrs *out; + + ALLOCATE (out, MeasurePanstarrs, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR = in[i].dR_PS; + out[i].dD = in[i].dD_PS; + out[i].M = in[i].M_PS; + out[i].dM = in[i].dM_PS; + out[i].Mcal = in[i].Mcal_PS; + out[i].Mgal = in[i].Mgal_PS; + out[i].airmass = in[i].airmass_PS; + out[i].dt = in[i].dt_PS; + out[i].FWx = in[i].FWx; + out[i].FWy = in[i].FWy; + out[i].theta = in[i].theta; + out[i].dophot = in[i].dophot; + out[i].source = in[i].source; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].flags = in[i].flags; + } + return (out); +} + +/* convert panstarrs-format averages to internal averages */ +Average *AveragePanstarrsToInternal (AveragePanstarrs *in, int Nvalues) { + + int i; + Average *out; + + ALLOCATE (out, Average, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R_PS = in[i].R; + out[i].D_PS = in[i].D; + out[i].M_PS = in[i].M; + out[i].dM_PS = in[i].dM; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Xg = in[i].Xg; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + } + return (out); +} + +/* convert internal averages to panstarrs-format averages */ +AveragePanstarrs *AverageInternalToPanstarrs (Average *in, int Nvalues) { + + int i; + AveragePanstarrs *out; + + ALLOCATE (out, AveragePanstarrs, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R_PS; + out[i].D = in[i].D_PS; + out[i].M = in[i].M_PS; + out[i].dM = in[i].dM_PS; + out[i].Xp = in[i].Xp; + out[i].Xm = in[i].Xm; + out[i].Xg = in[i].Xg; + out[i].Nm = in[i].Nm; + out[i].Nn = in[i].Nn; + out[i].code = in[i].code; + out[i].offset = in[i].offset; + out[i].missing = in[i].missing; + } + return (out); +} + +/* convert panstarrs-format secfilts to internal secfilts */ +SecFilt *SecFiltPanstarrsToInternal (SecFiltPanstarrs *in, int Nvalues) { + + int i; + SecFilt *out; + + ALLOCATE (out, SecFilt, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M_PS = in[i].M; + out[i].dM_PS = in[i].dM; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert internal secfilts to panstarrs-format secfilts */ +SecFiltPanstarrs *SecFiltInternalToPanstarrs (SecFilt *in, int Nvalues) { + + int i; + SecFiltPanstarrs *out; + + ALLOCATE (out, SecFiltPanstarrs, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M = in[i].M_PS; + out[i].dM = in[i].dM_PS; + out[i].Xm = in[i].Xm; + } + return (out); +} + +/* convert panstarrs-format images to internal images */ +Image *ImagePanstarrsToInternal (ImagePanstarrs *in, int Nvalues) { + + int i; + Image *out; + + ALLOCATE (out, Image, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz_PS = in[i].secz; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit_PS = in[i].apmifit; + out[i].dapmifit_PS = in[i].dapmifit; + out[i].Mcal_PS = in[i].Mcal; + out[i].dMcal_PS = in[i].dMcal; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} + +/* convert internal images to panstarrs-format images */ +ImagePanstarrs *ImageInternalToPanstarrs (Image *in, int Nvalues) { + + int i; + ImagePanstarrs *out; + + ALLOCATE (out, ImagePanstarrs, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + memcpy (out[i].dummy, in[i].dummy, 20); + strcpy (out[i].name, in[i].name); + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz_PS; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit_PS; + out[i].dapmifit = in[i].dapmifit_PS; + out[i].Mcal = in[i].Mcal_PS; + out[i].dMcal = in[i].dMcal_PS; + out[i].Xm = in[i].Xm; + out[i].source = in[i].source; + out[i].exptime = in[i].exptime; + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].code = in[i].code; + out[i].ccdnum = in[i].ccdnum; + out[i].order = in[i].order; + out[i].Mx = in[i].Mx; + out[i].My = in[i].My; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].Mxxx = in[i].Mxxx; + out[i].Mxxy = in[i].Mxxy; + out[i].Mxyy = in[i].Mxyy; + out[i].Myyy = in[i].Myyy; + out[i].Mxxxx = in[i].Mxxxx; + out[i].Mxxxy = in[i].Mxxxy; + out[i].Mxxyy = in[i].Mxxyy; + out[i].Mxyyy = in[i].Mxyyy; + out[i].Myyyy = in[i].Myyyy; + } + return (out); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image.c (revision 5388) @@ -0,0 +1,120 @@ +# include + +/* determine db mode (RAW/MEF) and read the complete table */ +int dvo_image_load (FITS_DB *db, int VERBOSE, int FORCE_READ) { + + int Naxis, Nread, status; + char buffer[241]; + + /* database name must be set first */ + if (db == NULL) return (FALSE); + if (db[0].f == NULL) return (FALSE); + + /* test for NAXIS == 2 in line 3 */ + bzero (buffer, 241); + fseek (db[0].f, 0, SEEK_SET); + Nread = fread (buffer, 1, 240, db[0].f); + if (Nread != 240) { + if (VERBOSE) fprintf (stderr, "can't determine image db mode\n"); + return (FALSE); + } + Nread = sscanf (&buffer[160], "NAXIS = %d", &Naxis); + if (Nread != 1) { + if (VERBOSE) fprintf (stderr, "can't determine image db mode\n"); + return (FALSE); + } + fseek (db[0].f, 0, SEEK_SET); + + /* default values */ + db[0].mode = DVO_MODE_MEF; + db[0].format = DVO_FORMAT_UNDEF; + if (Naxis == 2) db[0].mode = DVO_MODE_RAW; /* image table can only be RAW or MEF */ + + switch (db[0].mode) { + case DVO_MODE_MEF: + if (VERBOSE) fprintf (stderr, "reading images (mode DVO_MODE_MEF)\n"); + status = fits_db_load (db); + break; + case DVO_MODE_RAW: + if (VERBOSE) fprintf (stderr, "reading images (mode DVO_MODE_RAW)\n"); + status = dvo_image_load_raw (db, VERBOSE, FORCE_READ); + break; + default: + fprintf (stderr, "error getting image mode\n"); + exit (2); + } + FtableToImage (&db[0].ftable, &db[0].theader, &db[0].format); + db[0].swapped = TRUE; /* table has internal byte-order */ + return (TRUE); +} + +/* write the complete image table to disk in the requested mode */ +int dvo_image_save (FITS_DB *db, int VERBOSE) { + + int status; + + /* convert from internal to requested external format */ + ImageToFtable (&db[0].ftable, &db[0].theader, db[0].format); + db[0].swapped = FALSE; + + /* write data in appropriate mode */ + switch (db[0].mode) { + case DVO_MODE_MEF: + case DVO_MODE_SPLIT: + status = fits_db_save (db); + break; + case DVO_MODE_RAW: + status = dvo_image_save_raw (db, VERBOSE); + break; + default: + fprintf (stderr, "error writing image mode\n"); + exit (2); + } + return (TRUE); +} + +int dvo_image_addrows (FITS_DB *db, Image *new, int Nnew) { + + int Nimages; + + /* adjust header */ + Nimages = 0; + fits_scan (&db[0].header, "NIMAGES", "%d", 1, &Nimages); + Nimages += Nnew; + fits_modify (&db[0].header, "NIMAGES", "%d", 1, Nimages); + + fits_table_to_vtable (&db[0].ftable, &db[0].vtable, 0, 0); + fits_vadd_rows (&db[0].vtable, (char *) new, Nnew, sizeof(Image)); + + /* check that primary header and table header agree */ + if (Nimages != db[0].theader.Naxis[1]) { + fprintf (stderr, "header / table length mismatch!\n"); + } + return (TRUE); +} + +/* write the vtable data to disk in the requested mode */ +int dvo_image_update (FITS_DB *db, int VERBOSE) { + + int status; + + /* convert from internal to requested external format */ + /* theader is modified to match output format */ + ImageToVtable (&db[0].vtable, &db[0].theader, db[0].format); + db[0].swapped = FALSE; + + /* write data in appropriate mode */ + switch (db[0].mode) { + case DVO_MODE_MEF: + case DVO_MODE_SPLIT: + status = fits_db_update (db); + break; + case DVO_MODE_RAW: + status = dvo_image_update_raw (db, VERBOSE); + break; + default: + fprintf (stderr, "error writing image mode\n"); + exit (2); + } + return (TRUE); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image_raw.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image_raw.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/dvo_image_raw.c (revision 5388) @@ -0,0 +1,154 @@ +# include + +/* the old Image.dat files used a fake FITS header defining a finite data block + * this function reads in the header and data, and creates an ftable and theader to match + */ + +int dvo_image_load_raw (FITS_DB *db, int VERBOSE, int FORCE_READ) { + + int Nimage, Ndata, size, ImageSize, nbytes, Nbytes; + struct stat filestatus; + char format[80], telescope[80]; + + /* read fits header from file - return FALSE on error */ + if (!fits_fread_header (db[0].f, &db[0].header)) { + fprintf (stderr, "can't read primary header\n"); + return (FALSE); + } + + /* determine image table format */ + db[0].format = DVO_FORMAT_UNDEF; + if (fits_scan (&db[0].header, "FORMAT", "%s", 1, format)) { + if (!strcmp (format, "INTERNAL")) db[0].format = DVO_FORMAT_INTERNAL; + if (!strcmp (format, "LONEOS")) db[0].format = DVO_FORMAT_LONEOS; + if (!strcmp (format, "ELIXIR")) db[0].format = DVO_FORMAT_ELIXIR; + if (!strcmp (format, "PANSTARRS")) db[0].format = DVO_FORMAT_PANSTARRS; + if (db[0].format != DVO_FORMAT_UNDEF) goto got_format; + } + if (fits_scan (&db[0].header, "TELESCOP", "%s", 1, telescope)) { + if (!strncmp (telescope, "LONEOS", strlen("LONEOS"))) { + db[0].format = DVO_FORMAT_LONEOS; + goto got_format; + } + if (!strncmp (telescope, "1.3m McGraw-Hill", strlen("1.3m McGraw-Hill"))) { + db[0].format = DVO_FORMAT_ELIXIR; + goto got_format; + } + } + if (VERBOSE) fprintf (stderr, "cannot determine image table format\n"); + return (FALSE); + +got_format: + /* find number of images */ + Nimage = 0; + fits_scan (&db[0].header, "NIMAGES", "%d", 1, &Nimage); + if (stat (db[0].filename, &filestatus) == -1) { + if (VERBOSE) fprintf (stderr, "ERROR: failed to get status of image catalog\n"); + exit (1); + } + + /* get datatype size */ + ImageSize = 0; + if (db[0].format == DVO_FORMAT_INTERNAL) ImageSize = sizeof(Image); + if (db[0].format == DVO_FORMAT_LONEOS) ImageSize = sizeof(ImageLoneos); + if (db[0].format == DVO_FORMAT_ELIXIR) ImageSize = sizeof(ImageElixir); + if (db[0].format == DVO_FORMAT_PANSTARRS) ImageSize = sizeof(ImagePanstarrs); + + /* check that filesize makes sense */ + size = Nimage*ImageSize + db[0].header.size; + if (size != filestatus.st_size) { + Ndata = (filestatus.st_size - db[0].header.size) / ImageSize; + if (VERBOSE) fprintf (stderr, "ERROR: image catalog has inconsistent size\n"); + if (VERBOSE) fprintf (stderr, "header: %d, data: %d\n", Nimage, Ndata); + if (!FORCE_READ) exit (1); + Nimage = Ndata; + } + + /* create a dummy ftable, theader set for this data */ + /* (original table has NAXIS = 2, change to 0) */ + /* fits_modify (&db[0].header, "NAXIS", "%d", 1, 0); */ + fits_create_matrix (&db[0].header, &db[0].matrix); + db[0].ftable.header = &db[0].theader; + + if (db[0].format == DVO_FORMAT_INTERNAL) fits_table_mkheader_Image (&db[0].theader); + if (db[0].format == DVO_FORMAT_LONEOS) fits_table_mkheader_ImageLoneos (&db[0].theader); + if (db[0].format == DVO_FORMAT_ELIXIR) fits_table_mkheader_ImageElixir (&db[0].theader); + if (db[0].format == DVO_FORMAT_PANSTARRS) fits_table_mkheader_ImagePanstarrs (&db[0].theader); + + /* read data from file */ + Nbytes = ImageSize*Nimage; + ALLOCATE (db[0].ftable.buffer, char, Nbytes); + nbytes = fread (db[0].ftable.buffer, 1, Nbytes, db[0].f); + if (nbytes != Nbytes) { + if (VERBOSE) fprintf (stderr, "ERROR: problem loading image catalog\n"); + exit (1); + } + + fits_modify (&db[0].theader, "NAXIS2", "%d", 1, Nimage); + db[0].theader.Naxis[1] = Nimage; + db[0].ftable.size = fits_matrix_size (&db[0].theader); + db[0].swapped = FALSE; /* table does not have internal byte-order */ + return (TRUE); +} + +/* write out image db elements from vtable */ +int dvo_image_update_raw (FITS_DB *db, int VERBOSE) { + + int i, Nx, Ny, Nrow, Nbytes; + int status, start, offset; + int *row; + + if (VERBOSE) fprintf (stderr, "writing out %d images\n", db[0].vtable.Nrow); + + /* position to start of file */ + Fseek (db[0].f, 0, SEEK_SET); + + /* write out complete header (no check on disk size?) */ + status = fwrite (db[0].header.buffer, 1, db[0].header.size, db[0].f); + if (status != db[0].header.size) { + fprintf (stderr, "ERROR: failed writing data to image header\n"); + exit (1); + } + + /** this code is identical to fits_fwrite_vtable, except without padding */ + Nrow = db[0].vtable.Nrow; + row = db[0].vtable.row; + fits_scan (db[0].vtable.header, "NAXIS1", "%d", 1, &Nx); + fits_scan (db[0].vtable.header, "NAXIS2", "%d", 1, &Ny); + + /* file pointer is at beginning of desired table data */ + start = ftell (db[0].f); + + for (i = 0; i < Nrow; i++) { + offset = start + Nx*row[i]; + fseek (db[0].f, offset, SEEK_SET); + Nbytes = fwrite (db[0].vtable.buffer[i], sizeof (char), Nx, db[0].f); + if (Nbytes != Nx) { return (FALSE); } + } + return (TRUE); +} + +/* write out complete image db table from ftable */ +int dvo_image_save_raw (FITS_DB *db, int VERBOSE) { + + int status, Nx, Ny, size, Nbytes; + + if (VERBOSE) fprintf (stderr, "writing out %d images\n", db[0].theader.Naxis[1]); + + /* position to start of file */ + Fseek (db[0].f, 0, SEEK_SET); + + /* write out complete header (no check on disk size?) */ + status = fwrite (db[0].header.buffer, 1, db[0].header.size, db[0].f); + if (status != db[0].header.size) { + fprintf (stderr, "ERROR: failed writing data to image header\n"); + exit (1); + } + + fits_scan (db[0].ftable.header, "NAXIS1", "%d", 1, &Nx); + fits_scan (db[0].ftable.header, "NAXIS2", "%d", 1, &Ny); + size = Nx * Ny; + Nbytes = fwrite (db[0].ftable.buffer, sizeof(char), size, db[0].f); + if (Nbytes != size) return (FALSE); + return (TRUE); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/fits_db.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/fits_db.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/fits_db.c (revision 5388) @@ -0,0 +1,229 @@ +# include + +/* init the db structure */ +int fits_db_init (FITS_DB *db) { + + db[0].f = NULL; + db[0].filename = NULL; + db[0].header.buffer = NULL; + db[0].matrix.buffer = NULL; + db[0].theader.buffer = NULL; + db[0].ftable.buffer = NULL; + db[0].ftable.header = &db[0].theader; + + db[0].vtable.header = NULL; + db[0].vtable.buffer = NULL; + db[0].vtable.row = NULL; + return (TRUE); +} + +/* create an empty db */ +int fits_db_create (FITS_DB *db) { + fits_init_header (&db[0].header); + db[0].header.extend = TRUE; + fits_create_header (&db[0].header); + fits_create_matrix (&db[0].header, &db[0].matrix); + fits_print (&db[0].header, "NEXTEND", "%d", 1, 1); + db[0].ftable.header = &db[0].theader; + return (TRUE); +} + +int fits_db_lock (FITS_DB *db, char *filename) { + + /* database name must be set first */ + if (filename == NULL) { + fprintf (stderr, "db file is not set\n"); + return (FALSE); + } + + /* database handle must be set first */ + if (db == NULL) { + fprintf (stderr, "db handle is not set\n"); + return (FALSE); + } + + /* lock & open database */ + db[0].f = fsetlockfile (filename, db[0].timeout, db[0].lockstate, &db[0].dbstate); + if (db[0].f == NULL) { + if (db[0].dbstate == LCK_MISSING) { + fprintf (stderr, "no data in db %s\n", filename); + } else { + fprintf (stderr, "cannot set lock on db file %s\n", filename); + } + return (FALSE); + } + + db[0].filename = strcreate (filename); + return (TRUE); +} + +/* load the complete db table into memory - load first extension, do not validate EXTNAME */ +int fits_db_load (FITS_DB *db) { + + /* database name must be set first */ + if (db == NULL) { + fprintf (stderr, "db handle is not set\n"); + return (FALSE); + } + + /* init & load in FITS table data - return FALSE on error */ + if (!fits_fread_header (db[0].f, &db[0].header)) { + fprintf (stderr, "can't read primary header\n"); + return (FALSE); + } + if (!fits_fread_matrix (db[0].f, &db[0].matrix, &db[0].header)) { + fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + if (!fits_fread_header (db[0].f, &db[0].theader)) { + fprintf (stderr, "can't read table header"); + return (FALSE); + } + if (!fits_fread_ftable_data (db[0].f, &db[0].ftable)) { + fprintf (stderr, "can't read table data"); + return (FALSE); + } + db[0].swapped = FALSE; /* table does not have internal byte-order */ + return (TRUE); +} + +/* load the data from Nrows from table starting at start; load header, etc if needed */ +int fits_db_load_segment (FITS_DB *db, int start, int Nrows) { + + int Nskip; + + /* database name must be set first */ + if (db == NULL) { + fprintf (stderr, "db handle is not set\n"); + return (FALSE); + } + + /* database name must be opened */ + if (db[0].f == NULL) { + fprintf (stderr, "db is not opened\n"); + return (FALSE); + } + + /* position to start of file */ + Fseek (db[0].f, 0, SEEK_SET); + + /* load or skip header */ + if (db[0].header.buffer == NULL) { + if (!fits_fread_header (db[0].f, &db[0].header)) { + fprintf (stderr, "can't read primary header\n"); + return (FALSE); + } + } else { + Nskip = db[0].header.size; + Fseek (db[0].f, Nskip, SEEK_CUR); + } + + /* load or skip matrix */ + if (db[0].matrix.buffer == NULL) { + if (!fits_fread_matrix (db[0].f, &db[0].matrix, &db[0].header)) { + fprintf (stderr, "can't read primary matrix"); + return (FALSE); + } + } else { + Nskip = fits_matrix_size (&db[0].header); + Fseek (db[0].f, Nskip, SEEK_CUR); + } + + /* load or skip table header */ + if (db[0].theader.buffer == NULL) { + if (!fits_fread_header (db[0].f, &db[0].theader)) { + fprintf (stderr, "can't read table header"); + return (FALSE); + } + } else { + Nskip = db[0].header.size; + Fseek (db[0].f, Nskip, SEEK_CUR); + } + + /* read table segment into vtable */ + if (!fits_fread_vtable_range (db[0].f, &db[0].vtable, start, Nrows)) { + fprintf (stderr, "can't read table data"); + return (FALSE); + } + return (TRUE); +} + +/* write complete db file */ +int fits_db_save (FITS_DB *db) { + + /* write all data to file */ + make_backup (db[0].filename); + Fseek (db[0].f, 0, SEEK_SET); + + if (!fits_fwrite_header (db[0].f, &db[0].header)) { + fprintf (stderr, "can't write primary header"); + return (FALSE); + } + if (!fits_fwrite_matrix (db[0].f, &db[0].matrix)) { + fprintf (stderr, "can't write primary matrix"); + return (FALSE); + } + if (!fits_fwrite_Theader (db[0].f, &db[0].theader)) { + fprintf (stderr, "can't write table header"); + return (FALSE); + } + if (!fits_fwrite_table (db[0].f, &db[0].ftable)) { + fprintf (stderr, "can't write table data"); + return (FALSE); + } + return (TRUE); +} + +/* write vtable to db file (also appends rows to the end of the table) */ +int fits_db_update (FITS_DB *db) { + + /* this section is not valid if we have changed the size of header, matrix, theader */ + + /* write subset to file */ + make_backup (db[0].filename); /*** drop this!!?? ***/ + Fseek (db[0].f, 0, SEEK_SET); + + /* do we revert to the old version if this fails? */ + if (!fits_fwrite_header (db[0].f, &db[0].header)) { + fprintf (stderr, "can't update primary header"); + return (FALSE); + } + if (!fits_fwrite_matrix (db[0].f, &db[0].matrix)) { + fprintf (stderr, "can't update primary matrix"); + return (FALSE); + } + if (!fits_fwrite_Theader (db[0].f, &db[0].theader)) { + fprintf (stderr, "can't update table header"); + return (FALSE); + } + if (!fits_fwrite_vtable (db[0].f, &db[0].vtable)) { + fprintf (stderr, "can't update table data"); + return (FALSE); + } + return (TRUE); +} + +/* free memory associated with db handle */ +int fits_db_free (FITS_DB *db) { + fits_free_header (&db[0].header); + fits_free_matrix (&db[0].matrix); + fits_free_header (&db[0].theader); + fits_free_table (&db[0].ftable); + fits_free_vtable (&db[0].vtable); + if (db[0].filename != NULL) { + free (db[0].filename); + db[0].filename = NULL; + } + return (TRUE); +} + +/* close the db files (close open file & unlock) */ +int fits_db_close (FITS_DB *db) { + if (db[0].f == NULL) return (TRUE); + fclearlockfile (db[0].filename, db[0].f, db[0].lockstate, &db[0].dbstate); + return (TRUE); +} + +/* for this to be atomic, need to unlink before we unlock */ +/* unlink file here if resulting file is empty? */ +/* if (db[0].Nrow == 0) && (lockstate != LCK_SOFT)) unlink (dBFile); */ Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/imreg_datatypes.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/imreg_datatypes.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/imreg_datatypes.c (revision 5388) @@ -0,0 +1,50 @@ +# include + +char typename[N_TYPE][32] = {"none", "object", "dark", "bias", "flat", "mask", "fringe", "scatter", "modes", "frpts", "any"}; +char typecode[N_TYPE] = {'x', 'o', 'd', 'b', 'f', 'm', 'r', 's', 'M', 'F'}; +char modename[N_MODE][32] = {"none", "MEF", "SPLIT", "SINGLE", "CUBE", "SLICE", "MODES"}; + +int get_image_type (char *name) { + + int i; + + for (i = 0; i < N_TYPE; i++) { + if (!strncasecmp (name, typename[i], strlen(typename[i]))) { + return (i); + } + } + + /* special cases */ + if (!strncasecmp (name, "SKYFLAT", strlen("SKYFLAT"))) { + return (T_FLAT); + } + + return (T_UNDEF); +} + +char *get_type_name (int type) { + + if (type < T_NONE) return ("undef"); + if (type >= N_TYPE) return ("undef"); + return (typename[type]); +} + +int get_image_mode (char *name) { + + int i; + + for (i = 0; i < N_MODE; i++) { + if (!strncasecmp (name, modename[i], strlen(modename[i]))) { + return (i); + } + } + return (M_UNDEF); +} + +char *get_mode_name (int mode) { + + if (mode < M_NONE) return ("undef"); + if (mode >= N_MODE) return ("undef"); + return (modename[mode]); +} + Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/mosaic_astrom.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/mosaic_astrom.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/mosaic_astrom.c (revision 5388) @@ -0,0 +1,121 @@ +# include + +/* chip-match table: j = ChipMatch[i], images[j] is DIS for images[i] WRP */ + +static int Ndis = 0; +static int *DISentry = NULL; +static e_time *DIStzero = NULL; +static int *ChipMatch = NULL; + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage (Image *images, int Nimages, int entry) { + + int status; + + if (ChipMatch == NULL) { + status = FindMosaicForImage_TableSearch (images, Nimages, entry); + return (status); + } + status = FindMosaicForImage_MatchSearch (images, Nimages, entry); + return (status); +} + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage_TableSearch (Image *images, int Nimages, int entry) { + + e_time tzero; + int i; + + /* search backwards for image with same time and type DIS */ + + if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { + /* not a wrp image, do nothing */ + return (TRUE); /* error or not */ + } + + tzero = images[entry].tzero; + + for (i = entry - 1; i >= 0; i--) { + if (images[i].tzero != tzero) continue; + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + + /* found a valid match */ + RegisterMosaic (&images[i].coords); + } + + for (i = entry + 1; i < Nimages; i++) { + if (images[i].tzero != tzero) continue; + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + + /* found a valid match */ + RegisterMosaic (&images[i].coords); + return (TRUE); + } + return (FALSE); +} + +/* given an image array and a current entry of type WRP, find matching DIS & Register it */ +int FindMosaicForImage_MatchSearch (Image *images, int Nimages, int entry) { + + int N; + + if (strcmp(&images[entry].coords.ctype[4], "-WRP")) { + /* not a wrp image, do nothing */ + return (TRUE); /* error or not */ + } + + if (ChipMatch == NULL) return (FALSE); + + N = ChipMatch[entry]; + if (N < 0) return (FALSE); + if (N >= Nimages) return (FALSE); + + RegisterMosaic (&images[N].coords); + return (TRUE); +} + +int BuildChipMatch (Image *images, int Nimages) { + + int i, j, NDIS; + + if (DISentry != NULL) free (DISentry); + if (DIStzero != NULL) free (DIStzero); + if (ChipMatch != NULL) free (ChipMatch); + + /* find all DIS images, save tzero (& instrument) */ + Ndis = 0; + NDIS = 100; + ALLOCATE (DISentry, int, NDIS); + ALLOCATE (DIStzero, e_time, NDIS); + + for (i = 0; i < Nimages; i++) { + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + DISentry[Ndis] = i; + DIStzero[Ndis] = images[i].tzero; + Ndis ++; + if (Ndis >= NDIS) { + NDIS += 100; + REALLOCATE (DISentry, int, NDIS); + REALLOCATE (DIStzero, e_time, NDIS); + } + } + + /* find all matched WRP images */ + ALLOCATE (ChipMatch, int, Nimages); + for (i = 0; i < Nimages; i++) { + ChipMatch[i] = -1; + if (strcmp(&images[i].coords.ctype[4], "-WRP")) continue; + for (j = 0; j < Ndis; j++) { + if (DIStzero[j] != images[i].tzero) continue; + ChipMatch[i] = DISentry[j]; + goto got_match; + } + + fprintf (stderr, "WARNING: can't find matching mosaic \n"); + ChipMatch[i] = -2; + + got_match: + continue; + } + return (TRUE); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/photfits.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/photfits.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/photfits.c (revision 5388) @@ -0,0 +1,258 @@ +# include +# define MAX_ORDER 3 + +/**** XXXX warning: these functions have not been corrected to handle + the change of Mcal from millimags to mags + I don't think most functions are using the Mcal polynomial terms in + any case... ****/ + +/* all terms of order > 0 are stored as 16bit floats + the zero-order term is stored as a short int (-32k,+32k) */ + +/* convert double to low-precision short int */ +short int putMi (double value) { + + double T; + unsigned int I; + + if (value > 0) { + T = 1000 * log10 (value); + I = MAX (MIN (T, 16383), -16382) + 0x3fff; + return (I); + } + if (value < 0) { + T = -value; + T = 1000 * log10 (T); + I = MAX (MIN (T, 16383), -16382) + 0xbfff; + return (I); + } + if (value == 0) { + I = 0; + return (0); + } + return (NO_MAG); +} + +/* convert low-precision short int to double */ +double getMi (short int value) { + + double T, V, sign; + short int I; + + if (value == 0) { + return (0.0); + } + if (value & 0x8000) { + sign = -1; + I = value & 0x7fff; + } else { + sign = +1; + I = value; + } + + T = (I - 0x3fff) / 1000.0; + V = sign * pow (10.0, T); + return (V); +} + +/* convert image parameters to c[i] coeffs */ +void returnMcal (Image *image, double *c) { + + switch (image[0].order) { + case 0: + c[0] = image[0].Mcal_PS; + return; + case 1: + c[0] = image[0].Mcal_PS; + c[1] = getMi (image[0].Mx); + c[2] = getMi (image[0].My); + return; + case 2: + c[0] = image[0].Mcal_PS; + c[1] = getMi (image[0].Mx); + c[2] = getMi (image[0].Mxx); + c[3] = getMi (image[0].My); + c[4] = getMi (image[0].Mxy); + c[5] = getMi (image[0].Myy); + return; + case 3: + c[0] = image[0].Mcal_PS; + c[1] = getMi (image[0].Mx); + c[2] = getMi (image[0].Mxx); + c[3] = getMi (image[0].Mxxx); + c[4] = getMi (image[0].My); + c[5] = getMi (image[0].Mxy); + c[6] = getMi (image[0].Mxxy); + c[7] = getMi (image[0].Myy); + c[8] = getMi (image[0].Mxyy); + c[9] = getMi (image[0].Myyy); + return; + case 4: + c[0] = image[0].Mcal_PS; + c[1] = getMi (image[0].Mx); + c[2] = getMi (image[0].Mxx); + c[3] = getMi (image[0].Mxxx); + c[4] = getMi (image[0].Mxxxx); + c[5] = getMi (image[0].My); + c[6] = getMi (image[0].Mxy); + c[7] = getMi (image[0].Mxxy); + c[8] = getMi (image[0].Mxxxy); + c[9] = getMi (image[0].Myy); + c[10] = getMi (image[0].Mxyy); + c[11] = getMi (image[0].Mxxyy); + c[12] = getMi (image[0].Myyy); + c[13] = getMi (image[0].Mxyyy); + c[14] = getMi (image[0].Myyyy); + return; + default: + c[0] = 0; + return; + } +} + +void assignMcal (Image *image, double *c, int order) { + + image[0].order = order; + + switch (order) { + case 0: + image[0].Mcal_PS = c[0]; + return; + case 1: + image[0].Mcal_PS = c[0]; + image[0].Mx = putMi(c[1]); + image[0].My = putMi(c[2]); + return; + case 2: + image[0].Mcal_PS = c[0]; + image[0].Mx = putMi(c[1]); + image[0].Mxx = putMi(c[2]); + image[0].My = putMi(c[3]); + image[0].Mxy = putMi(c[4]); + image[0].Myy = putMi(c[5]); + return; + case 3: + image[0].Mcal_PS = c[0]; + image[0].Mx = putMi(c[1]); + image[0].Mxx = putMi(c[2]); + image[0].Mxxx = putMi(c[3]); + image[0].My = putMi(c[4]); + image[0].Mxy = putMi(c[5]); + image[0].Mxxy = putMi(c[6]); + image[0].Myy = putMi(c[7]); + image[0].Mxyy = putMi(c[8]); + image[0].Myyy = putMi(c[9]); + return; + case 4: + image[0].Mcal_PS = c[0]; + image[0].Mx = putMi(c[1]); + image[0].Mxx = putMi(c[2]); + image[0].Mxxx = putMi(c[3]); + image[0].Mxxxx = putMi(c[4]); + image[0].My = putMi(c[5]); + image[0].Mxy = putMi(c[6]); + image[0].Mxxy = putMi(c[7]); + image[0].Mxxxy = putMi(c[8]); + image[0].Myy = putMi(c[9]); + image[0].Mxyy = putMi(c[10]); + image[0].Mxxyy = putMi(c[11]); + image[0].Myyy = putMi(c[12]); + image[0].Mxyyy = putMi(c[13]); + image[0].Myyyy = putMi(c[14]); + return; + default: + image[0].Mcal_PS = 0.0; + image[0].order = 0; + return; + } +} + +/* return value of image fit at x,y */ +double applyMcal (Image *image, double x, double y) { + + double Mcal, c[15]; + + returnMcal (image, c); /* convert image parameters to c[i] coeffs */ + Mcal = 0.0; + switch (image[0].order) { + case 0: + Mcal = c[0]; + break; + case 1: + Mcal = c[0] + x*c[1] + y*c[2]; + break; + case 2: + Mcal = c[0] + x*(c[1] + x*c[2]) + y*(c[3] + x*c[4] + y*c[5]); + break; + case 3: + Mcal = c[0] + x*(c[1] + x*(c[2] + x*c[3])) + y*(c[4] + x*(c[5] + x*c[6]) + y*(c[7] + x*c[8] + y*c[9])); + break; + case 4: + Mcal = c[0] + x*(c[1] + x*(c[2] + x*(c[3] + x*c[4]))) + y*(c[5] + x*(c[6] + x*(c[7] + x*c[8])) + y*(c[9] + x*(c[10] + x*c[11]) + y*(c[12] + x*c[13] + y*c[14]))); + break; + } + return (Mcal); +} + +double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order) { + + int i; + double *x, *y, *m, *dm; + double S, S2, s, dS; + + S2 = S = 0; + if (N < 2) { + return (0.0); + } + + x = X; y = Y; m = M; dm = dM; + switch (order) { + case 0: + for (i = 0; i < N; i++, m++, dm++) { + dS = *m - c[0]; + *dm = dS; + S += dS; + S2 += dS*dS; + } + break; + case 1: + for (i = 0; i < N; i++, m++, x++, y++, dm++) { + s = c[0] + *x*c[1] + *y*c[2]; + dS = *m - s; + *dm = dS; + S += dS; + S2 += dS*dS; + } + break; + case 2: + for (i = 0; i < N; i++, m++, x++, y++, dm++) { + s = c[0] + *x*(c[1] + *x*c[2]) + *y*(c[3] + *x*c[4] + *y*c[5]); + dS = *m - s; + *dm = dS; + S += dS; + S2 += dS*dS; + } + break; + case 3: + for (i = 0; i < N; i++, m++, x++, y++, dm++) { + s = c[0] + *x*(c[1] + *x*(c[2] + *x*c[3])) + *y*(c[4] + *x*(c[5] + *x*c[6]) + *y*(c[7] + *x*c[8] + *y*c[9])); + dS = *m - s; + *dm = dS; + S += dS; + S2 += dS*dS; + } + break; + case 4: + for (i = 0; i < N; i++, m++, x++, y++, dm++) { + s = c[0] + *x*(c[1] + *x*(c[2] + *x*(c[3] + *x*c[4]))) + *y*(c[5] + *x*(c[6] + *x*(c[7] + *x*c[8])) + *y*(c[9] + *x*(c[10] + *x*c[11]) + *y*(c[12] + *x*c[13] + *y*c[14]))); + dS = *m - s; + *dm = dS; + S += dS; + S2 += dS*dS; + } + break; + } + S = S / N; + S2 = sqrt (S2 / N - S*S); + return (S2); +} Index: /tags/libdvo-1-0/Ohana/src/libdvo/src/skydb.c =================================================================== --- /tags/libdvo-1-0/Ohana/src/libdvo/src/skydb.c (revision 5388) +++ /tags/libdvo-1-0/Ohana/src/libdvo/src/skydb.c (revision 5388) @@ -0,0 +1,400 @@ + +typedef struct { + float r, d; +} SkyCoord; + +/* each region is bounded on the sky by lines of constant RA & DEC. + the region represents a specified depth, and it may or may not have + children. The start and end of the children in the array are childS and childE + If the region at the given depth is populated with an object table, then 'object' is TRUE. + If the region at the given depth is populated with an image table, then 'image' is TRUE. +*/ + +typedef struct { + char name[16]; + float Rmin, Rmax; + float Dmin, Dmax; + int childS, childE; + char depth, child; + char object, image; +} SkyRegion; /* 48 bytes */ + +/* I have defined no accelerators other than the table hierarchy */ + +/* find region which overlaps c at given depth (-1 : max depth) */ +SkyRegion *SkyFindPoint (SkyRegion *db, SkyCoord c, int depth) { + + Ns = 0; + Ne = 1; + + while (1) { + No = -1; + for (i = Ns; (No == -1) && (i < Ne); i++) { + if (c.r < db[i].Rmin) continue; + if (c.r > db[i].Rmax) continue; + if (c.d < db[i].Dmin) continue; + if (c.d > db[i].Dmax) continue; + // got the needed region at this depth + No = i; + } + if (No == -1) return ((SkyRegion *) NULL); + if (depth == db[No].depth) return (&db[No]); + if ((depth > db[No].depth) && !db[No].child) return ((SkyRegion *) NULL); + + /* need to check Ns, Ne, or guarantee valid range */ + Ns = db[No].childS; + Ne = db[No].childE; + + if ((depth > db[No].depth) && db[No].child) continue; + return (&db[No]); + } +} + +/* find regions at all levels which overlap c */ +SkyRegion **SkyFindLevels (SkyRegion *db, SkyCoord c, int *Nregion) { + + SkyRegion **region; + + Ns = 0; + Ne = 1; + + N = 0; + NREGION = 10; + ALLOCATE (region, SkyRegion *, NREGION); + + while (1) { + No = -1; + for (i = Ns; (No == -1) && (i < Ne); i++) { + if (c.r < db[i].Rmin) continue; + if (c.r > db[i].Rmax) continue; + if (c.d < db[i].Dmin) continue; + if (c.d > db[i].Dmax) continue; + // got the needed region at this depth + No = i; + } + if (No == -1) return ((SkyRegion *) NULL); + + region[N] = &db[No]; + N++; + if (N == NREGION) { + NREGION += 10; + REALLOCATE (region, SkyRegion *, NREGION); + } + + /* need to check Ns, Ne, or guarantee valid range */ + if (db[No].child) { + Ns = db[No].childS; + Ne = db[No].childE; + } else { + *Nregion = N; + return (region); + } + } +} + +/* find regions contained within rectangular region c1 - c2 */ +SkyRegion **SkyFindArea (SkyRegion *db, SkyCoord c1, SkyCoord c2, *nlist) { + + int Nlist; + SkyRegion *list, *new, *sub; + + while (c1.r > 360.0) c1.r -= 360.0; + while (c1.r < 0.0) c1.r += 360.0; + while (c2.r > 360.0) c2.r -= 360.0; + while (c2.r < 0.0) c2.r += 360.0; + + /* check on c1.r > c2.r : cross boundary */ + if (c1.r > c2.r) { + c0 = c2; c0.r = 360.0; + list1 = SkyFindArea (db, c1, c0, &Nlist1); + c0 = c1; c0.r = 0.0; + list2 = SkyFindArea (db, c0, c2, &Nlist2); + Nlist = Nlist1 + Nlist2; + ALLOCATE (list, SkyRegion *, Nlist); + memcpy (&list[0], list1, Nlist1*sizeof()); + memcpy (&list[Nlist1], list2, Nlist2*sizeof()); + free (list1); + free (list2); + *nlist = Nlist; + return (list); + } + + Nlist = 1; + ALLOCATE (list, SkyRegion *, 1); + list[0] = &db[0]; + getchild = db[0].child; + + while (getchild) { + + getchild = FALSE; + Nnew = 0; + NNEW = Nlist + 100; + ALLOCATE (new, SkyRegion *, NNEW); + + for (i = 0; i < Nlist; i++) { + if (list[i][0].child) { + sub = SkyFindAreaDB(db, list[i], c1, c2, &Nsub); + if (Nnew + Nsub == NNEW) { + NNEW += 100; + REALLOCATE (new, SkyRegion *, NNEW); + } + for (i = 0; i < Nsub; i++) { + getchild |= sub[i][0].child; + new[Nnew] = sub[i]; + Nnew++; + } + free (sub); + } else { + new[Nnew] = list[i]; + Nnew ++; + if (Nnew == NNEW) { + NNEW += 100; + REALLOCATE (new, SkyRegion *, NNEW); + } + } + } + + free (list); + list = new; + Nlist = Nnew; + } + + return (list); + *nlist = Nlist; + +} + +/* ?? */ +SkyRegion *SkyFindAreaDB (SkyRegion *db, SkyRegion *ref, SkyCoord c1, SkyCoord c2, int *Nregion) { + + SkyRegion **region; + + Ns = ref[0].childS; + Ne = ref[0].childE; + + N = 0; + NREGION = 100; + ALLOCATE (region, SkyRegion *, NREGION); + + /* c1 min, c2 max */ + + for (i = Ns; (No == -1) && (i < Ne); i++) { + if (c2.r < db[i].Rmin) continue; + if (c1.r > db[i].Rmax) continue; + if (c2.d < db[i].Dmin) continue; + if (c1.d > db[i].Dmax) continue; + + region[N] = &db[i]; + N++; + if (N == NREGION) { + NREGION += 100; + REALLOCATE (region, SkyRegion *, NREGION); + } + } + *Nregion = N; + return (region); +} + +SkyRegion *SkyBuildTable (SkyRegion *seed, int Nseed, int level, int depth) { + + /* given a table of Rmin, Dmin, Rmax, Dmax, name, at level 3, generate + all higher levels and 'depth' extra levels */ + + /* levels: + 0 - fullsky.cpt + 1 - n????.cpt / s????.cpt + 2 - r????.cpt + 3 - ????.cpt + 4 - ????.??.cpt + */ + + SkyRegion *db; + + /* allocate at least 30 for levels 0 & 1 */ + NREGION = 100; + ALLOCATE (db, SkyRegion, NREGION); + + /* full sky */ + strcpy (db[0].name, "fullsky.cpt"); + db[0].Rmin = 0; db[0].Rmax = 360; + db[0].Dmin = -90; db[0].Dmax = 90; + db[0].depth = 0; + db[0].child = FALSE; + N = 1; + + db[0].childS = N; + /* north dec bands */ + for (dec = 0; dec < 90; dec += 7.5) { + sprintf (db[N].name, "n%04d.cpt", (int) 100*dec); + db[N].Rmin = 0; db[N].Rmax = 360; + db[N].Dmin = dec; db[N].Dmax = dec + 7.5; + db[N].depth = 1; + db[N].child = FALSE; + N++; + } + /* south dec bands */ + for (dec = 0; dec > -90; dec -= 7.5) { + sprintf (db[N].name, "s%04d.cpt", (int) 100*dec); + db[N].Rmin = 0; db[N].Rmax = 360; + db[N].Dmin = dec - 7.5; db[N].Dmax = dec; + db[N].depth = 1; + db[N].child = FALSE; + N++; + } + db[0].childE = N; + + /* subdivide dec bands based on seed */ + Rnumber = 0; + childS = db[0].childS; + childE = db[0].childE; + for (i = childS; i < childE; i++) { + + Nnew = 0; + NNEW = 100; + ALLOCATE (new, SkyRegion, NNEW); + + for (j = 0; j < Nseed; j++) { + if (seed[j].Rmin > db[i].Rmax) continue; + if (seed[j].Rmax < db[i].Rmin) continue; + if (seed[j].Dmin > db[i].Dmax) continue; + if (seed[j].Dmax < db[i].Dmin) continue; + + for (k = 0; k < Nnew; k++) { + if ((seed[j].Rmin == new[k].Rmin) && + (seed[j].Rmax == new[k].Rmax)) { + goto next_seed; + } + if ((seed[j].Rmin == new[k].Rmin) ^^ + (seed[j].Rmax != new[k].Rmax)) { + fprintf (stderr, "inconsistent blocks in seed file\n"); + return (NULL); + } + } + new[Nnew].Rmin = seed[j].Rmin; + new[Nnew].Rmax = seed[j].Rmax; + new[Nnew].Dmin = db[i].Dmin; + new[Nnew].Dmax = db[i].Dmax; + new[Nnew].depth = 2; + new[Nnew].child = FALSE; + strncpy (root, db[i].name, 5); root[5] = 0; + sprintf (new[Nnew].name, "%s/r%04d.cpt", root, Rnumber); + Rnumber ++; + Nnew ++; + if (Nnew == NNEW) { + NNEW += 100; + REALLOCATE (new, SkyRegion, NNEW); + } + next_seed: + } + /* update db list */ + db[i].childS = N; + db[i].childE = N + Nnew; + NREGION = N + Nnew + 100; + REALLOCATE (db, SkyRegion, NREGION); + memcpy (&db[N], new, Nnew*sizeof(SkyRegion)); + free (new); + N += Nnew; + } + + /* subdivide ra strips based on seed */ + childS = db[db[0].childS].childS; + childE = db[db[0].childE - 1].childE; + nextS = N; + for (i = childS; i < childE; i++) { + + Nnew = 0; + NNEW = 100; + ALLOCATE (new, SkyRegion, NNEW); + + for (j = 0; j < Nseed; j++) { + if (seed[j].Rmin > db[i].Rmax) continue; + if (seed[j].Rmax < db[i].Rmin) continue; + if (seed[j].Dmin > db[i].Dmax) continue; + if (seed[j].Dmax < db[i].Dmin) continue; + + for (k = 0; k < Nnew; k++) { + if ((seed[j].Dmin == new[k].Dmin) && + (seed[j].Dmax == new[k].Dmax)) { + goto next_seed; + } + if ((seed[j].Dmin == new[k].Dmin) ^^ + (seed[j].Dmax != new[k].Dmax)) { + fprintf (stderr, "inconsistent blocks in seed file\n"); + return (NULL); + } + } + new[Nnew].Dmin = seed[j].Dmin; + new[Nnew].Dmax = seed[j].Dmax; + new[Nnew].Rmin = db[i].Rmin; + new[Nnew].Rmax = db[i].Rmax; + new[Nnew].depth = 3; + new[Nnew].child = FALSE; + strcpy (new[Nnew].name, seed[j].name); + Nnew ++; + if (Nnew == NNEW) { + NNEW += 100; + REALLOCATE (new, SkyRegion, NNEW); + } + next_seed: + } + + /* update db list */ + db[i].childS = N; + db[i].childE = N + Nnew; + NREGION = N + Nnew + 100; + REALLOCATE (db, SkyRegion, NREGION); + memcpy (&db[N], new, Nnew*sizeof(SkyRegion)); + free (new); + N += Nnew; + } + nextE = N; + + /* subdivide entries once more */ + childS = nextS; + childE = nextE; + for (i = childS; i < childE; i++) { + Rnumber = 0; + db[i].childS = N; + dDec = (db[i].Dmax - db[i].Dmin) / 5.0; + dRa = (db[i].Rmax - db[i].Rmin) / 5.0; + for (nx = 0; nx < 5; nx++) { + for (ny = 0; ny < 5; ny++) { + db[N].Dmin = db[i].Dmin + dDec * (nx + 0); + db[N].Dmax = db[i].Dmax + dDec * (nx + 1); + db[N].Rmin = db[i].Rmin + dDec * (nx + 0); + db[N].Rmax = db[i].Rmax + dDec * (nx + 1); + db[N].depth = 3; + db[N].child = FALSE; + strncpy (root, db[i].name, 10); root[10] = 0; + sprintf (db[N].name, "%s.%02d.cpt", root, Rnumber); + Rnumber ++; + N ++; + if (N == NREGION) { + NREGION += 100; + REALLOCATE (new, SkyRegion, NREGION); + } + } + } + } + +} + +/* region is pointer to entry in db */ +SkyRegion *SkyExtend (SkyRegion *db, SkyRegion *region) { + +} + +/* region is pointer to entry in db */ +SkyRegion *SkyContract (SkyRegion *db, SkyRegion *region) { + +} + +SkyRegion *SkyFindGCircle (SkyRegion *db, SkyCoord c1, SkyCoord c2) { + +} + + +/* is entire db in memory, or do we load from disk for each operation? + 10,000 at level 3, 100,000 @ 4 : 0.5MB, 5MB */ + +