Index: /trunk/Ohana/src/libohana/Makefile =================================================================== --- /trunk/Ohana/src/libohana/Makefile (revision 5236) +++ /trunk/Ohana/src/libohana/Makefile (revision 5237) @@ -22,13 +22,10 @@ libohana: $(LIB)/libohana.$(ARCH).a -INC1 = \ +INCS = \ $(DESTINC)/Xohana.h \ $(DESTINC)/ohana_allocate.h \ $(DESTINC)/ohana.h -INC2 = \ -$(DESTINC)/dvo.h - -OBJ1 = \ +OBJS = \ $(SRC)/ohana_allocate.$(ARCH).o \ $(SRC)/string.$(ARCH).o \ @@ -37,59 +34,8 @@ $(SRC)/time.$(ARCH).o \ $(SRC)/config.$(ARCH).o \ +$(SRC)/Fread.$(ARCH).o \ $(SRC)/version.$(ARCH).o -OBJ2 = \ -$(SRC)/Fread.$(ARCH).o \ -$(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 - -ACOM = $(ADEF)/autocode.c $(ADEF)/autocode.h $(ADEF)/common.h - -$(OBJ1): $(INC1) -$(OBJ2): $(INC2) -$(AOBJ): $(ASRC)/%.$(ARCH).o : $(ADEF)/%.d $(ACOM) - cd $(AUTO) && make install +$(OBJS): $(INCS) $(DESTLIB)/libohana.a: $(LIB)/libohana.$(ARCH).a @@ -98,8 +44,8 @@ cp $(LIB)/libohana.$(ARCH).a $(DESTLIB)/libohana.a -$(LIB)/libohana.$(ARCH).a: $(OBJ1) $(AOBJ) $(OBJ2) +$(LIB)/libohana.$(ARCH).a: $(OBJS) @if [ ! -d $(LIB) ]; then mkdir -p $(LIB); fi rm -f $(LIB)/libohana.$(ARCH).a - ar rcv $(LIB)/libohana.$(ARCH).a $(OBJ1) $(AOBJ) $(OBJ2) + ar rcv $(LIB)/libohana.$(ARCH).a $(OBJS) $(RANLIB) $(LIB)/libohana.$(ARCH).a @@ -114,5 +60,4 @@ rm -f */*.o rm -f */*.a - cd $(AUTO) && make clean dist: clean Index: /trunk/Ohana/src/libohana/doc/ChangeLog.txt =================================================================== --- /trunk/Ohana/src/libohana/doc/ChangeLog.txt (revision 5236) +++ /trunk/Ohana/src/libohana/doc/ChangeLog.txt (revision 5237) @@ -1,2 +1,15 @@ + +2005.10.07 + + I was having some memory collision problems, and attempting to + use the ohana_allocate functions reminded me that the libFITS + functions were not supported under ohana_allocate. This was + unhelpful. I bit the bullet and split libohana into libohana + (base functions only, including ohana_allocate) and libdvo + (functions based on the libautocode structures). Doing this + allowed me to make libFITS depend on libohana (including + ohana_allocate). BUT, this forced me to change all LDFLAGS + entries in ohana to swap -lohana -lFITS for -lFITS -lohana, + and to add include in some cases. libohana-1-4 Index: /trunk/Ohana/src/libohana/include/ohana.h =================================================================== --- /trunk/Ohana/src/libohana/include/ohana.h (revision 5236) +++ /trunk/Ohana/src/libohana/include/ohana.h (revision 5237) @@ -20,6 +20,4 @@ /* OHANA included stuff */ -# include - # ifndef OHANA_H # define OHANA_H @@ -152,3 +150,12 @@ */ +/* +# ifndef ANSI +# include +# else +# include +# include # endif +*/ + +# endif Index: /trunk/Ohana/src/libohana/src/Fread.c =================================================================== --- /trunk/Ohana/src/libohana/src/Fread.c (revision 5236) +++ /trunk/Ohana/src/libohana/src/Fread.c (revision 5237) @@ -1,7 +1,5 @@ # include -# include int ByteSwap (char *ptr, int size, int nitems, char *type); -int ConvertStruct (char *buffer, int size, int Nbytes, char *type); # define BYTE_EXCHANGE(X,Y) tmp = byte[X]; byte[X] = byte[Y]; byte[Y] = tmp; @@ -17,31 +15,5 @@ tmp = byte[X+3]; byte[X+3] = byte[X+4]; byte[X+4] = tmp; -# define REGIMAGE_SIZE 360 -# define DETREG_SIZE 416 -# define PHOTPARS_SIZE 108 -# define MISSING_SIZE 4 -# define SECFILT_SIZE 6 -# define SPECTRUM_SIZE 216 /* where is this?? */ -# define IMAGE_SIZE 240 - -# if (PANSTARRS) -# define AVERAGE_SIZE 40 -# define MEASURE_SIZE 32 -# endif - -# if (ELIXIR) -# define AVERAGE_SIZE 32 -# define MEASURE_SIZE 32 -# endif - -# if (LONEOS) -# define AVERAGE_SIZE 28 -# define MEASURE_SIZE 20 -# endif - -# define OLDPHOTPARS_SIZE 100 -# define OLDSECFILT_SIZE 4 -# define OLDDETREG_SIZE 288 - +/** this is also used in libautocode/def/common.h -- consolidate **/ # ifdef linux # define BYTE_SWAP @@ -52,9 +24,10 @@ # endif -/* add other architectures here, ie dec, alpha, etc */ # ifdef dec # define BYTE_SWAP # endif +/* add other architectures here, ie dec, alpha, etc */ + int Fread (void *ptr, int size, int nitems, FILE *f, char *type) { @@ -81,17 +54,10 @@ } -int ConvertStruct (char *buffer, int size, int nitems, char *type) { - - int status; - status = ByteSwap (buffer, size, nitems, type); - return (status); -} - int ByteSwap (char *ptr, int size, int nitems, char *type) { # ifdef BYTE_SWAP - int i, j; - unsigned char *byte, *byte0, *byte1, *byte2, *byte3, *byte4, *byte5, *byte6, *byte7, tmp; + int i; + unsigned char *byte0, *byte1, *byte2, *byte3, *byte4, *byte5, *byte6, *byte7, tmp; if (!strcmp (type, "char")) return (TRUE); @@ -150,4 +116,5 @@ } +# if (0) /** now in autocode **/ if (!strcmp (type, "regimage")) { if (size != REGIMAGE_SIZE) { @@ -352,5 +319,5 @@ } -# endif /** ELIXIR **/ +# endif /** PANSTARRS **/ @@ -494,4 +461,5 @@ } # endif /** others **/ +# endif /** now in autocode **/ fprintf (stderr, "unknown type %s\n", type); Index: unk/Ohana/src/libohana/src/LoadPhotcodes.c =================================================================== --- /trunk/Ohana/src/libohana/src/LoadPhotcodes.c (revision 5236) +++ (revision ) @@ -1,799 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/coordops.c =================================================================== --- /trunk/Ohana/src/libohana/src/coordops.c (revision 5236) +++ (revision ) @@ -1,564 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_catalog.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_catalog.c (revision 5236) +++ (revision ) @@ -1,130 +1,0 @@ -# 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, char mode, 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, mode, VERBOSE); - break; - case DVO_MODE_MEF: - if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_MEF)\n"); - load_catalog_mef (catalog, mode, VERBOSE); - break; - case DVO_MODE_SPLIT: - if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_SPLIT)\n"); - load_catalog_split (catalog, mode, 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: unk/Ohana/src/libohana/src/dvo_catalog_mef.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_catalog_mef.c (revision 5236) +++ (revision ) @@ -1,237 +1,0 @@ -# include -# include - -int load_catalog_mef (Catalog *catalog, char mode, 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 (mode & 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 (mode & 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 (mode & 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 (mode & 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: unk/Ohana/src/libohana/src/dvo_catalog_raw.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_catalog_raw.c (revision 5236) +++ (revision ) @@ -1,581 +1,0 @@ -# include -# include - -/* read data from raw-style catalog file; set data format values based on file */ - -int load_catalog_raw (Catalog *catalog, char mode, 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 (mode & 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 (mode & 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 (mode & 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 (mode & 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 (mode & 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: unk/Ohana/src/libohana/src/dvo_catalog_split.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_catalog_split.c (revision 5236) +++ (revision ) @@ -1,637 +1,0 @@ -# include -# include - -int load_catalog_split (Catalog *catalog, char mode, 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 (mode & 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 (mode & 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 (mode & 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); - } - - /* 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 (mode & 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 (mode & 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: unk/Ohana/src/libohana/src/dvo_convert.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_convert.c (revision 5236) +++ (revision ) @@ -1,439 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_convert_elixir.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_convert_elixir.c (revision 5236) +++ (revision ) @@ -1,246 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_convert_loneos.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_convert_loneos.c (revision 5236) +++ (revision ) @@ -1,244 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_convert_panstarrs.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_convert_panstarrs.c (revision 5236) +++ (revision ) @@ -1,248 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_image.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_image.c (revision 5236) +++ (revision ) @@ -1,121 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/dvo_image_raw.c =================================================================== --- /trunk/Ohana/src/libohana/src/dvo_image_raw.c (revision 5236) +++ (revision ) @@ -1,154 +1,0 @@ -# include -# 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 */ - 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: unk/Ohana/src/libohana/src/fits_db.c =================================================================== --- /trunk/Ohana/src/libohana/src/fits_db.c (revision 5236) +++ (revision ) @@ -1,230 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/imreg_datatypes.c =================================================================== --- /trunk/Ohana/src/libohana/src/imreg_datatypes.c (revision 5236) +++ (revision ) @@ -1,51 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/mosaic_astrom.c =================================================================== --- /trunk/Ohana/src/libohana/src/mosaic_astrom.c (revision 5236) +++ (revision ) @@ -1,122 +1,0 @@ -# include -# 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: /trunk/Ohana/src/libohana/src/ohana_allocate.c =================================================================== --- /trunk/Ohana/src/libohana/src/ohana_allocate.c (revision 5236) +++ /trunk/Ohana/src/libohana/src/ohana_allocate.c (revision 5237) @@ -198,5 +198,4 @@ } - fprintf (stderr, "checking for memory corruption...\n"); header = FALSE; for (i = 0; i < Nmemlist; i++) { @@ -216,4 +215,5 @@ } fprintf (stderr, "corrupt: %3d %c %c %-15s %3d\n", i, bottom, top, memlist[i].file, memlist[i].line); + abort (); } return; Index: unk/Ohana/src/libohana/src/photfits.c =================================================================== --- /trunk/Ohana/src/libohana/src/photfits.c (revision 5236) +++ (revision ) @@ -1,259 +1,0 @@ -# include -# 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: unk/Ohana/src/libohana/src/skydb.c =================================================================== --- /trunk/Ohana/src/libohana/src/skydb.c (revision 5236) +++ (revision ) @@ -1,400 +1,0 @@ - -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 */ - -