Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/Makefile
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/Makefile (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/Makefile (revision 21775)
@@ -0,0 +1,61 @@
+include ../../Configure
+HOME = .
+
+LIB = $(HOME)/lib
+SRC = $(HOME)/src
+MAN = $(HOME)/doc
+INC = $(HOME)/include
+DESTBIN = $(LBIN)
+DESTLIB = $(LLIB)
+DESTINC = $(LINC)
+DESTMAN = $(LMAN)
+
+INCS = -I$(LINC)
+CFLAGS = -Wunused -o $*.$(ARCH).o $(INCS)
+LFLAGS = -c $(INCS)
+
+default: install
+
+install: $(DESTLIB)/libohana.a
+
+ohana: $(LIB)/libohana.$(ARCH).a
+
+HOBJ = $(DESTINC)/ohana.h $(DESTINC)/loneos.h $(DESTINC)/Xohana.h
+
+LOBJ = \
+$(SRC)/string.$(ARCH).o $(SRC)/findexec.$(ARCH).o \
+$(SRC)/config.$(ARCH).o $(SRC)/coordops.$(ARCH).o \
+$(SRC)/Fread.$(ARCH).o $(SRC)/glockfile.$(ARCH).o \
+$(SRC)/LoadPhotcodes.$(ARCH).o $(SRC)/photfits.$(ARCH).o \
+$(SRC)/phot_catalog.$(ARCH).o $(SRC)/imreg_datatypes.$(ARCH).o
+
+$(DESTLIB)/libohana.a: $(LIB)/libohana.$(ARCH).a
+ @if [ ! -d $(DESTLIB) ]; then mkdir -p $(DESTLIB); fi
+ rm -f $(DESTLIB)/libohana.a
+ cp $(LIB)/libohana.$(ARCH).a $(DESTLIB)/libohana.a
+
+$(LIB)/libohana.$(ARCH).a: $(HOBJ) $(LOBJ)
+ @if [ ! -d $(LIB) ]; then mkdir -p $(LIB); fi
+ rm -f $(LIB)/libohana.$(ARCH).a
+ ar rcv $(LIB)/libohana.$(ARCH).a $(LOBJ)
+ $(RANLIB) $(LIB)/libohana.$(ARCH).a
+
+$(LOBJ): $(DESTINC)/ohana.h $(DESTINC)/loneos.h
+
+$(DESTINC)/%: $(INC)/%
+ @if [ ! -d $(DESTINC) ]; then mkdir -p $(DESTINC); fi
+ rm -f $(DESTINC)/$*
+ cp $(INC)/$* $(DESTINC)/
+
+clean:
+ rm -f */*~
+ rm -f */#*
+ rm -f */*.o
+ rm -f */*.a
+
+dist: clean
+
+.SUFFIXES: .$(ARCH).o
+
+.c.$(ARCH).o:
+ $(CC) $(CFLAGS) -c $<
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ChangeLog.txt
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ChangeLog.txt (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ChangeLog.txt (revision 21775)
@@ -0,0 +1,8 @@
+
+libohana-1-1
+ dropped the configuration source ./.(name)rc : was just confusing
+ in glockfile : made closing a NULL file OK
+ added imreg_datatypes (part of imregister-1-1 datatype reorganization)
+
+libohana-1-0
+ import to CVS
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/VERSIONS
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/VERSIONS (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/VERSIONS (revision 21775)
@@ -0,0 +1,5 @@
+
+tag names used by libohana:
+
+TAG : Comment
+libohana-1-0 : first version under CVS
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coords_minimal.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coords_minimal.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coords_minimal.c (revision 21775)
@@ -0,0 +1,115 @@
+/*
+here is some C code to project from R,D to x,y (all in decimal
+degrees) in a SIN projection, with Ro, Do as projection center:
+*/
+
+# define FULLPROJECTION 0
+# define DEG_RAD 57.295779513082322
+# define RAD_DEG 0.017453292519943
+
+RD_to_XY (double *x, double *y, double R, double D, double Ro, double Do) {
+
+ double sdp, cdp, salp, calp, sdel, cdel, stht, sphi, cphi;
+
+ sdp = sin(RAD_DEG*Do);
+ cdp = cos(RAD_DEG*Do);
+ salp = sin(RAD_DEG*(ra - Ro));
+ calp = cos(RAD_DEG*(ra - Ro));
+ 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) { return 0; /* projection from the wrong side of the sphere */ }
+
+ X = DEG_RAD * sphi;
+ Y = -DEG_RAD * cphi;
+
+# if (FULLPROJECTION)
+
+ /*
+ these lines allow for a rotation / distortion 2x2 matrix (pci_j),
+ a (two direction) plate-scale shift (cdelt1, cdelt2),
+ and a reference center offset of Xo, Yo, if desired.
+ */
+
+ tmp_d = 1.0 / (pc_1_1*pc_2_2 - pc_1_2*pc_2_1);
+ *x = tmp_d * (pc_2_2*X - pc_1_2*Y) / cdelt1 + Xo;
+ *y = tmp_d * (pc_1_1*Y - pc_2_1*X) / cdelt2 + Yo;
+
+# else
+
+ *x = X;
+ *y = Y;
+
+# endif
+
+
+ return (1);
+
+}
+
+
+/*
+here is some C code to project from x,y to R,D (all in decimal
+degrees) in a SIN projection, with Ro, Do as projection center:
+*/
+
+XY_to_RD (double *ra, double *dec, double x, double y, double Ro, double Do) {
+
+
+ double L, M, X, Y, T, Z;
+ double R, sphi, cphi, stht, ctht;
+ double alpha, delta, salp, calp, sdel, sdp, cdp;
+
+ *ra = 0;
+ *dec = 0;
+ stht = ctht = 1;
+
+
+# if (FULLPROJECTION)
+ /*
+ these lines allow for a rotation / distortion 2x2 matrix (pci_j),
+ a (two direction) plate-scale shift (cdelt1, cdelt2),
+ and a reference center offset of Xo, Yo, if desired.
+ */
+ X = cdelt1 * (x - Xo);
+ Y = cdelt2 * (y - Yo);
+
+ L = (X*pc1_1 + Y*pc1_2);
+ M = (X*pc2_1 + Y*pc2_2);
+# else
+ L = x;
+ M = y;
+# endif
+
+ R = hypot (L,M);
+ if ((L == 0) && (M == 0)) {
+ sphi = 0;
+ cphi = 1;
+ }
+ else {
+ sphi = L / R;
+ cphi = -M / R;
+ }
+
+ ctht = RAD_DEG * R;
+ stht = sqrt (1 - ctht*ctht);
+
+ sdp = sin(RAD_DEG*Do);
+ cdp = cos(RAD_DEG*Do);
+
+ 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 + Ro;
+ *dec = DEG_RAD*delta;
+
+ return (1);
+
+}
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coordtrans.txt
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coordtrans.txt (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/coordtrans.txt (revision 21775)
@@ -0,0 +1,57 @@
+
+{
+
+ 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);
+ R = hypot (L,M);
+ if ((L == 0) && (M == 0)) {
+ sphi = 0;
+ cphi = 1;
+ } else {
+ sphi = L / R;
+ cphi = -M / R;
+ }
+ 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);
+ }
+ 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;
+}
+
+RD_to_XY () {
+
+ 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;
+
+ X = DEG_RAD * sphi / stht;
+ Y = -DEG_RAD * cphi / stht;
+ tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+ *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1;
+ *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2;
+
+}
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/libohana.html
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/libohana.html (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/libohana.html (revision 21775)
@@ -0,0 +1,148 @@
+
+ libohana API reference
+
+ functions in string.c
+
+ int get_argument (int argc, char **argv, char *arg);
+
+search argv, argc list (starting at argv[1]) for desired argument.
+returns entry number if found, 0 otherwise.
+
+
+ int remove_argument (int Narg, int &argc, char **argv);
+
+remove the given entry Narg from the argv list, adjusting argc
+appropriately.
+
+
+ int stripwhite (char *string);
+
+strip whitespace (space or tab) from start and end of string.
+
+
+ int strnumcmp (char *str1, char *str2);
+
+compare str1 and str2 as strings, or as numbers if both are pure
+numeric values (base 10 only). returns TRUE / FALSE.
+
+
+ char *strcreate (char *string);
+ char *strncreate (char *string, int Nbytes);
+ int scan_line (FILE *f, char *line);
+ int dparse (double *value, int Nentry, char *line);
+ int fparse (double *value, int Nentry, char *line);
+ char *_parse_nextword (char *string);
+
+ functions in findexec.c
+
+ char *pathname (char *file);
+ char *filerootname (char *file);
+ char *fileextname (char *file);
+ char *filebasename (char *file);
+ char *findexec (int argc, char **argv);
+ int mkdirhier (char *path);
+ char *getcwd_cfht (char *path, int size);
+
+ functions in glockfile.c
+
+ FILE *fsetlockfile (char *filename, double timeout, int type, int *state);
+ int fclearlockfile (char *filename, FILE *f, int type, int *state);
+ int setlockfile2 (char *filename, double timeout, int type, int *state);
+ int clearlockfile2 (char *filename, double timeout, int type, int *state);
+
+ deprecated functions : are they all removed?
+ int setlockfile (char *name, double timeout, int hard);
+ int clearlockfile (char *name, int fd, int hard);
+ int gnfsflush (char *filename);
+
+ functions in Fread.c
+
+ int Fread (char ptr, int size, int nitems, FILE *f, char *type);
+ int Fwrite (char ptr, int size, int nitems, FILE *f, char *type);
+ int ByteSwap (char *ptr, int size, int nitems, char *type);
+ int ConvertStruct (char *buffer, int size, int nitems, char *type);
+
+ functions in LoadPhotcodes.c
+
+ int LoadPhotcodes (char *filename, PhotCodeData *photcodes);
+ int FreePhotcodes (PhotCodeData *photcodes);
+ int GetPhotcodebyName (PhotCodeData *photcodes, char *name);
+ int GetPhotcodeRefs (PhotCodeData *photcodes, char *name, int *code, int *ref);
+ int GetPhotcodeEquiv (PhotCodeData *photcodes, char *name, int *code, int *ref);
+ char *GetPhotnamebyCode (PhotCodeData *photcodes, int code);
+
+ functions in config.c
+
+ char *SelectConfigFile (int *argc, char **argv, char *progname);
+ char *LoadConfigFile (char *filename);
+ char *ScanConfig (char *config, char *field, char *mode, int Nentry, ...);
+ char *expandline (char *line, char *config);
+ char *LoadRawConfigFile (char *filename, int options);
+
+ functions in coordops.c
+
+ int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords);
+ int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords);
+ int fXY_to_RD (float *ra, float *dec, float x, float y, Coords *coords);
+ int fRD_to_XY (float *x, float *y, float ra, float dec, Coords *coords);
+ int GetCoords (Coords *coords, Header *header);
+ int coords_precess (double *ra, double *dec, double in_epoch, double out_epoch);
+
+ functions in phot_catalog.c
+
+ int lock_catalog (Catalog *catalog, int lockmode);
+ int unlock_catalog (Catalog *catalog);
+ int load_catalog (Catalog *catalog, char mode, int VERBOSE);
+ int save_catalog (Catalog *catalog, char VERBOSE);
+ Measure *FixOldMeasure (OldMeasure *in, int Nvalues);
+ Average *FixOldAverage (OldAverage *in, int Nvalues);
+ SecFilt *FixOldSecFilt (OldSecFilt *in, int Nvalues);
+
+ functions in photfits.c
+
+ short int putMi (double value);
+ double getMi (short int value);
+ void returnMcal (Image *image, double *c);
+ void assignMcal (Image *image, double *c, int order);
+ double applyMcal (Image *image, double x, double y);
+ double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order);
+
+ defined macros
+
+the ohana library includes the following macros:
+
+
+ALLOCATE (void *ptr, type, int Nelement)
+REALLOCATE (void *ptr, type, int Nelement)
+
+SIGN(X)
+ROUND(X)
+SQR(X)
+SQ(X)
+MIN(X,Y)
+MAX(X,Y)
+SWAP(X,Y)
+
+whitespace (char c)
+
+
+ defined constants
+
+
+TRUE
+FALSE
+
+LCK_SOFT - block writing
+LCK_XCLD - block reading and writing
+LCK_HARD - block reading and writing, persistent
+
+LCK_UNLOCK - file is unlocked
+LCK_ACCESS - can't get access to file
+LCK_TIMEOUT - timeout setting lock
+LCK_HARDLCK - error setting hard lockfile
+LCK_HUNLOCK - error clearing hard lockfile
+
+LCK_EMPTY - locked file is empty
+LCK_FULL - locked file is not empty
+LCK_UNKNOWN - can't stat file to get size
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/load_catalog.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/load_catalog.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/load_catalog.c (revision 21775)
@@ -0,0 +1,263 @@
+# include
+# include
+
+/* the function in this file uses the old definition of Catalog,
+ and is deprecated by the new locking scheme */
+
+typedef struct {
+ Average *average;
+ Measure *measure;
+ Missing *missing;
+ SecFilt *secfilt;
+ Header header;
+ int Naverage, Nmeasure, Nmissing, Nsecfilt;
+ int *found;
+ int *image;
+ float *X, *Y;
+} Catalog;
+
+Measure *FixOldMeasure (OldMeasure *in, int Nvalues) {
+
+ int i;
+ Measure *out;
+
+ ALLOCATE (out, Measure, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].dR = in[i].dR;
+ out[i].dD = in[i].dD;
+ out[i].M = in[i].M;
+ out[i].dM = in[i].dM;
+ out[i].Mcal = in[i].Mcal;
+ out[i].dophot = in[i].dophot;
+ out[i].source = in[i].source;
+ out[i].t = in[i].t;
+
+ out[i].dt = 0xffff;
+ out[i].averef = in[i].average & 0x00ffffff;
+ out[i].flags = (in[i].average & 0xff000000) >> 24;
+ }
+
+ free (in);
+
+ return (out);
+}
+
+Average *FixOldAverage (OldAverage *in, int Nvalues) {
+
+ int i;
+ Average *out;
+
+ ALLOCATE (out, Average, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].R = in[i].R;
+ out[i].D = in[i].D;
+ out[i].M = in[i].M;
+ 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].dM = 0xffff;
+ }
+
+ free (in);
+
+ return (out);
+}
+
+SecFilt *FixOldSecFilt (OldSecFilt *in, int Nvalues) {
+
+ int i;
+ SecFilt *out;
+
+ ALLOCATE (out, SecFilt, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].M = in[i].M;
+ out[i].Xm = in[i].Xm;
+
+ out[i].dM = 0xffff;
+ }
+
+ free (in);
+
+ return (out);
+}
+
+load_catalog (char *filename, Catalog *catalog, char mode, char VERBOSE) {
+
+ int Nitems, nitems;
+ int i, Nmeas, Nmiss, size, ssize;
+ int NewMeasure, NewAverage, NewSecFilt;
+ FILE *f;
+ OldMeasure *tmpmeasure;
+ OldAverage *tmpaverage;
+ OldSecFilt *tmpsecfilt;
+ struct stat filestatus;
+
+ /* possible exit status:
+ 0 - serious failure
+ 1 - success
+ 2 - file non-existent */
+
+ /* read catalog header */
+ if (!fits_read_header (filename, &catalog[0].header)) {
+ if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", filename);
+ return (2);
+ }
+
+ f = fopen (filename, "r");
+ if (f == (FILE *) NULL) {
+ if (VERBOSE) fprintf (stderr, "ERROR: can't read data from catalog file: %s\n", filename);
+ return (FALSE);
+ }
+ fseek (f, catalog[0].header.size, SEEK_SET);
+
+ 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;
+ if (!fits_scan (&catalog[0].header, "NEWMEAS", "%t", 1, &NewMeasure)) NewMeasure = FALSE;
+ if (!fits_scan (&catalog[0].header, "NEWAVES", "%t", 1, &NewAverage)) NewAverage = FALSE;
+ if (!fits_scan (&catalog[0].header, "NEWSECF", "%t", 1, &NewSecFilt)) NewSecFilt = FALSE;
+
+ /* predicted file size */
+
+ size = 0;
+ ssize = NewAverage ? sizeof (Average) : sizeof (OldAverage);
+ size += ssize * catalog[0].Naverage;
+
+ ssize = NewMeasure ? sizeof (Measure) : sizeof (OldMeasure);
+ size += ssize * catalog[0].Nmeasure;
+
+ ssize = sizeof (Missing);
+ size += ssize * catalog[0].Nmissing;
+
+ ssize = NewSecFilt ? sizeof (SecFilt) : sizeof (OldSecFilt);
+ size += ssize * catalog[0].Nsecfilt * catalog[0].Naverage;
+
+ size += catalog[0].header.size;
+
+ /* check that file size makes sense */
+ if (stat (filename, &filestatus) == -1) {
+ fprintf (stderr, "ERROR: failed to get status of star catalog\n");
+ exit (1);
+ }
+ if (size != filestatus.st_size) {
+ fprintf (stderr, "ERROR: star catalog has inconsistent size\n");
+ exit (1);
+ }
+
+ if (catalog[0].Naverage == 0) {
+ if (VERBOSE) fprintf (stderr, "no stars yet in catalog %s\n", filename);
+ fclose (f);
+ return (TRUE);
+ }
+
+ if (mode & LOAD_AVES) {
+ /* read average values */
+ if (NewAverage) {
+ ALLOCATE (catalog[0].average, Average, MAX (catalog[0].Naverage, 1));
+ Nitems = catalog[0].Naverage;
+ nitems = Fread (catalog[0].average, sizeof(Average), Nitems, f, "average");
+ } else {
+ ALLOCATE (tmpaverage, Average, MAX (catalog[0].Naverage, 1));
+ Nitems = catalog[0].Naverage;
+ nitems = Fread (tmpaverage, sizeof(OldAverage), Nitems, f, "oldaverage");
+ catalog[0].average = FixOldAverage (tmpaverage, Nitems);
+ }
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to read averages from catalog file %s\n", filename);
+ fclose (f);
+ return (FALSE);
+ }
+ } else {
+ /* skip over averages */
+ fseek (f, catalog[0].Naverage * sizeof(Average), SEEK_CUR);
+ }
+
+ if (mode & LOAD_MEAS) {
+ /* read measurements */
+ if (NewMeasure) {
+ ALLOCATE (catalog[0].measure, Measure, MAX (catalog[0].Nmeasure, 1));
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fread (catalog[0].measure, sizeof(Measure), Nitems, f, "measure");
+ } else {
+ ALLOCATE (tmpmeasure, OldMeasure, MAX (catalog[0].Nmeasure, 1));
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fread (tmpmeasure, sizeof(OldMeasure), Nitems, f, "oldmeasure");
+ catalog[0].measure = FixOldMeasure (tmpmeasure, Nitems);
+ }
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to read measures from catalog file %s\n", filename);
+ fclose (f);
+ return (FALSE);
+ }
+ } else {
+ /* skip over measures */
+ fseek (f, catalog[0].Nmeasure * sizeof(Measure), SEEK_CUR);
+ }
+
+ if (mode & LOAD_MISS) {
+ /* read missing */
+ ALLOCATE (catalog[0].missing, Missing, MAX (catalog[0].Nmissing, 1));
+ Nitems = catalog[0].Nmissing;
+ nitems = Fread (catalog[0].missing, sizeof(Missing), Nitems, f, "missing");
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to read missing from catalog file %s\n", filename);
+ fclose (f);
+ return (FALSE);
+ }
+ } else {
+ /* skip over missings */
+ fseek (f, catalog[0].Nmissing * sizeof(Missing), SEEK_CUR);
+ }
+
+ if (mode & LOAD_SECF) {
+ /* read missing */
+ if (NewSecFilt) {
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ ALLOCATE (catalog[0].secfilt, SecFilt, MAX (Nitems, 1));
+ nitems = Fread (catalog[0].secfilt, sizeof(SecFilt), Nitems, f, "secfilt");
+ } else {
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ ALLOCATE (tmpsecfilt, OldSecFilt, MAX (Nitems, 1));
+ nitems = Fread (tmpsecfilt, sizeof(OldSecFilt), Nitems, f, "oldsecfilt");
+ catalog[0].secfilt = FixOldSecFilt (tmpsecfilt, Nitems);
+ }
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to read secondary filters from catalog file %s\n", filename);
+ fclose (f);
+ return (FALSE);
+ }
+ } else {
+ /* skip over secfilts */
+ fseek (f, catalog[0].Naverage * catalog[0].Nsecfilt * sizeof(SecFilt), SEEK_CUR);
+ }
+
+ fclose (f);
+
+ if (VERBOSE) fprintf (stderr, "read %d stars from catalog file %s (%d measurements, %d missing, %d secondary filters)\n",
+ catalog[0].Naverage, filename, catalog[0].Nmeasure, catalog[0].Nmissing, catalog[0].Nsecfilt);
+
+ 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)) {
+ fprintf (stderr, "ERROR: data in catalog %s is corrupt, sums don't check\n", filename);
+ fprintf (stderr, "ERROR: Nmeas: %d, %d\n", Nmeas, catalog[0].Nmeasure);
+ fprintf (stderr, "ERROR: Nmiss: %d, %d\n", Nmiss, catalog[0].Nmissing);
+ return (FALSE);
+ }
+ }
+ return (TRUE);
+}
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/notes.txt
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/notes.txt (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/notes.txt (revision 21775)
@@ -0,0 +1,45 @@
+
+string / argv functions:
+
+int get_argument (int argc, char **argv, char *arg);
+int remove_argument (int Narg, int &argc, char **argv);
+int stripwhite (char *string);
+int strnumcmp (char *str1, char *str2);
+char *strcreate (char *string);
+char *strncreate (char *string, int Nbytes);
+int scan_line (FILE *f, char *line);
+int dparse (double *value, int Nentry, char *line);
+
+char *_parse_nextword (char *string);
+
+char *pathname (char *file);
+char *filerootname (char *file);
+char *fileextname (char *file);
+char *filebasename (char *file);
+char *findexec (int argc, char **argv);
+int mkdirhier (char *path);
+char *getcwd_cfht (char *path, int size);
+
+
+locking functions / catalog
+
+int setlockfile2 (char *filename, double timeout, int type, int *state);
+int clearlockfile2 (char *filename, double timeout, int type, int *state);
+
+int lock_catalog (Catalog *catalog, int lockmode);
+int unlock_catalog (Catalog *catalog);
+int load_catalog (Catalog *catalog, char mode, int VERBOSE);
+void save_catalog (Catalog *catalog, char VERBOSE);
+
+locks:
+
+LCK_SOFT - block writing
+LCK_XCLD - block reading and writing
+LCK_HARD - block reading and writing, persistent
+
+setlockfile (LCK_SOFT):
+
+errors:
+ HARD,XCLD: if file missing, create (666?)
+ SOFT: if file missing: error LCK_ACCESS
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ohana.txt
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ohana.txt (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/ohana.txt (revision 21775)
@@ -0,0 +1,57 @@
+
+these routines make up the ohana library. to compile a program with
+them, include ohana.h and link with -libohana.
+
+the ohana library includes the following macros:
+
+ALLOCATE (void *ptr, type, int Nelement)
+REALLOCATE (void *ptr, type, int Nelement)
+
+SIGN(X)
+ROUND(X)
+SQR(X)
+SQ(X)
+MIN(X,Y)
+MAX(X,Y)
+SWAP(X,Y)
+
+whitespace (char c)
+
+the ohana library includes the following functions:
+
+int stripwhite (char *string)
+char *strcreate (char *string)
+char *strncreate (char *string, int n)
+int scan_line (FILE *f, char *line)
+char *_parse_nextword (char *string)
+int dparse (double *X, int NX, char *line)
+
+int get_argument (int argc, char **argv, char *arg)
+int remove_argument (int N, int *argc, char **argv)
+
+int Fread (char ptr, int size, int nitems, FILE *f, char *type)
+int Fwrite (char ptr, int size, int nitems, FILE *f, char *type)
+int ByteSwap (char *ptr, int size, int nitems, char *type)
+
+int _check_permissions (char *filename)
+char *pathname (char *name)
+char *filebasename (char *name)
+char *findexec (int argc, char **argv)
+
+char *SelectConfigFile (int *argc, char **argv, char *progname)
+char *LoadConfigFile (char *filename)
+char *ScanConfig (char *config, char *field, char *mode, void *ptr)
+char *expandline (char *line, char *config)
+
+XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords)
+RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords)
+fXY_to_RD (float *ra, float *dec, float x, float y, Coords *coords)
+fRD_to_XY (float *x, float *y, float ra, float dec, Coords *coords)
+int GetCoords (Coords *coords, Header *header)
+
+save_catalog (char *filename, Catalog *catalog, char VERBOSE)
+load_catalog (char *filename, Catalog *catalog, char mode, char VERBOSE)
+
+LoadPhotcodes (char *filename, PhotCodeData *photcodes)
+GetPhotcodebyName (PhotCodeData *photcodes, char *name)
+double applyphot (PhotCode photcode, int Nphotcode, Average *ave, Measure *meas)
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/save_catalog.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/save_catalog.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/doc/save_catalog.c (revision 21775)
@@ -0,0 +1,89 @@
+# include
+# include
+
+/* the function in this file uses the old definition of Catalog,
+ and is deprecated by the new locking scheme */
+
+typedef struct {
+ Average *average;
+ Measure *measure;
+ Missing *missing;
+ SecFilt *secfilt;
+ Header header;
+ int Naverage, Nmeasure, Nmissing, Nsecfilt;
+ int *found;
+ int *image;
+ float *X, *Y;
+} Catalog;
+
+save_catalog (char *filename, Catalog *catalog, char VERBOSE) {
+
+ int i, Nitems, nitems, mode;
+ FILE *f;
+
+ if (catalog[0].Naverage == 0) {
+ if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n");
+ return (2);
+ }
+
+ f = fopen (filename, "w");
+ if (f == NULL) {
+ fprintf (stderr, "ERROR: can't create new catalog file: %s\n", filename);
+ return (0);
+ }
+
+ mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
+ chmod (filename, mode);
+
+ 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, "NEWMEAS", "%t", 1, TRUE);
+ fits_modify (&catalog[0].header, "NEWAVES", "%t", 1, TRUE);
+ fits_modify (&catalog[0].header, "NEWSECF", "%t", 1, TRUE);
+
+ nitems = Fwrite (catalog[0].header.buffer, 1, catalog[0].header.size, f, "char");
+ if (nitems != catalog[0].header.size) {
+ fprintf (stderr, "ERROR: failed to write header\n");
+ return (0);
+ }
+
+ Nitems = catalog[0].Naverage;
+ nitems = Fwrite (catalog[0].average, sizeof(Average), Nitems, f, "average");
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to write catalog file aves %s\n", filename);
+ return (0);
+ }
+
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fwrite (catalog[0].measure, sizeof(Measure), Nitems, f, "measure");
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to write catalog file meas %s\n", filename);
+ return (0);
+ }
+
+ Nitems = catalog[0].Nmissing;
+ nitems = Fwrite (catalog[0].missing, sizeof(Missing), Nitems, f, "missing");
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to write catalog file miss %s\n", filename);
+ return (0);
+ }
+
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ nitems = Fwrite (catalog[0].secfilt, sizeof(SecFilt), Nitems, f, "secfilt");
+ if (nitems != Nitems) {
+ fprintf (stderr, "ERROR: failed to write catalog file secondary filters %s\n", filename);
+ return (0);
+ }
+
+ fclose (f);
+ return (1);
+
+}
+
+ /* possible exit status:
+ 0 - serious failure
+ 1 - success
+ 2 - data non-existent */
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/include/Xohana.h
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/include/Xohana.h (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/include/Xohana.h (revision 21775)
@@ -0,0 +1,8 @@
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/include/loneos.h
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/include/loneos.h (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/include/loneos.h (revision 21775)
@@ -0,0 +1,467 @@
+
+/*** this file uses data types which must have fixed sizes regardless
+ of the platform. It originally used the basic C primitives:
+ float, double, int, short int, unsigned long int, etc.
+ this breaks under 64 bit (and probably on other systems).
+ I should define internal data types which should be set by the
+ use of # define statements if needed. I will cheat for now and use
+ the time_t to replace unsigned long int in this file
+***/
+
+# ifndef LONEOS_H
+# define LONEOS_H
+
+# define e_time unsigned int
+
+/* image data modes in RegImage */
+enum {T_UNDEF = -1, T_NONE, T_OBJECT, T_DARK, T_BIAS, T_FLAT, T_MASK, T_FRINGE, T_SCATTER, T_MODES, T_FRINGEPTS, T_ANY, N_TYPE};
+enum {M_UNDEF = -1, M_NONE, M_MEF, M_SPLIT, M_SINGLE, M_CUBE, M_SLICE, M_MODES, N_MODE};
+
+# if (0)
+ # define MODE_NONE 0 /* not defined */
+ # define MODE_MEF 1 /* chip of mef image */
+ # define MODE_SPLIT 2 /* chip of split image */
+ # define MODE_SINGLE 3 /* single chip, not mosaic */
+ # define MODE_CUBE 4 /* cube entry */
+ # define MODE_SLICE 5 /* slice of a cube */
+ # define MODE_MODES 6 /* not used? */
+
+ # define TYPE_NONE 0
+ # define TYPE_OBJECT 1
+ # define TYPE_DARK 2
+ # define TYPE_BIAS 3
+ # define TYPE_FLAT 4
+/** these are no longer used by imregister-3.0
+ status:imrough uses the imregister-3.0 values
+ status:imsearch should migrate to the imregister.h version
+ or be removed!
+**/
+# endif
+
+
+# define IMREG_DIST 0x01 /* image distributed, only imregister-3.0 */
+
+/* catalog values to be loaded */
+# define LOAD_AVES 0x01
+# define LOAD_MEAS 0x02
+# define LOAD_MISS 0x04
+# define LOAD_SECF 0x08
+
+/* invalid mag value */
+# define NO_MAG 0x7fff
+# define NO_ERR 0xff
+
+/* photometry code types */
+# define PHOT_PRI 0x01
+# define PHOT_SEC 0x02
+# define PHOT_DEP 0x03
+# define PHOT_REF 0x04
+# define PHOT_ALT 0x05 /* never stored, only for look-ups */
+
+/* Image.code values. these are codes to note bad images */
+# define ID_IMAGE_NEW 0x0000 /* no nrphot attempted */
+# define ID_IMAGE_NOCAL 0x0001 /* used within nrphot to mean "don't apply fit" */
+# define ID_IMAGE_POOR 0x0002 /* relphot says image is bad */
+# define ID_IMAGE_SKIP 0x0004 /* external information image is bad */
+# define ID_IMAGE_FEW 0x0008 /* currently too few measurements for good value */
+
+/* Measure.flags values */
+# define ID_MEAS_NOCAL 0x0001
+# define ID_MEAS_POOR 0x0002
+# define ID_MEAS_SKIP 0x0004
+# define ID_MEAS_AREA 0x0008
+# define BLEND_IMAGE 0x0100
+# define BLEND_CATALOG 0x0200
+# define BLEND_IMAGE_NEIGHBOR 0x1000
+# define ID_MEAS_TRAIL 0x2000
+# define ID_MEAS_GHOST 0x4000
+
+/* Average.code values */
+# define ID_STAR_FEW 0x0001 /* used within relphot: skip star */
+# define ID_STAR_POOR 0x0002 /* used within relphot: skip star */
+# define ID_PROPER 0x0400 /* star with large proper motion */
+# define ID_TRANSIENT 0x1000 /* is this mutually exclusive with USNO? */
+# define ID_VARIABLE 0x2000 /* not currently set? */
+# define ID_ASTEROID 0x2000 /* identified with an asteroid */
+# define ID_BAD_OBJECT 0x4000 /* if all measurements are bad, set this bit */
+# define ID_MOVING 0x8000
+# define ID_ROCK 0xa000 /* 0x8000 + 0x2000 */
+# define ID_GHOST 0xc001 /* 0x8000 + 0x4000 + 0x0001 */
+# define ID_TRAIL 0xc002 /* 0x8000 + 0x4000 + 0x0002 */
+# define ID_BLEED 0xc003 /* 0x8000 + 0x4000 + 0x0003 */
+# define ID_COSMIC 0xc004 /* 0x8000 + 0x4000 + 0x0004 */
+
+/*
+ BLEND_IMAGE: the star on an image is matched with more
+ than one star in the catalog (image has worse seeing than catalog)
+ BLEND_CATALOG: the star in the catalog is matched with more
+ than one star on the image (image has better seeing than catalog)
+ CALIBRATED: relative photometry has been performed on this measurement
+ BLEND_IMAGE_NEIGHBOR: the star on an image is matched with more
+ than one star in the catalog, but not in the same catalog file.
+*/
+
+/* structure for data on a catalog region */
+typedef struct {
+ char filename[256];
+ double DEC[2], RA[2];
+} GSCRegion;
+
+/* structure for Image Registration Database */
+typedef struct {
+ char filename[64];
+ char pathname[128];
+ char filter[32];
+ char instrument[32];
+ char ccd; /* an identifier for CCD in mosaic (0 for single ccd) */
+ char mode; /* MEF, SINGLE, SPLIT */
+ char type; /* image type */
+ char flag; /* image status flags */
+ char junk[24];
+
+ float exptime;
+ float airmass;
+ float sky;
+ float bias;
+ float fwhm;
+
+ float telfocus;
+ float xprobe, yprobe, zprobe;
+ float dettemp;
+ float teltemp[4];
+ float rotangle;
+ float ra, dec;
+
+ e_time obstime;
+ e_time regtime;
+} RegImage; /* 360 bytes / image */
+
+/* structure for Image Registration Database */
+typedef struct {
+ float sky;
+ float bias;
+ float fwhm;
+ float exptime;
+ float airmass;
+ e_time obstime;
+ short int ccdnum;
+ char junk1, junk2;
+ char filter[32];
+} RufImage; /* 60 bytes */
+
+/* Old Detrend Database structure */
+typedef struct {
+ e_time tstart;
+ e_time tstop;
+ e_time treg;
+ float sigma, clipsigma;
+ int type;
+ int filter;
+ int ccd;
+ char filename[254];
+ char flag1, flag2;
+} DetRegOld; /* 288 bytes */
+
+/* Detrend Database structure */
+typedef struct {
+ e_time tstart;
+ e_time tstop;
+ e_time treg;
+ float exptime;
+ int type;
+ int filter;
+ int ccd;
+ int Nentry;
+ int Norder;
+ char mode;
+ char altpath; /* true: data on alt db paths */
+ char dummy[58]; /* for future expansion */
+ char label[64];
+ char filename[256];
+} DetReg; /* 416 bytes */
+
+/* Instant Zeropoint structure */
+typedef struct {
+ float ZP; /* measured zero point */
+ float ZPo; /* expected zero point */
+ float dZP; /* error on zero point */
+ float K; /* airmass term used */
+ float X; /* color term used */
+ e_time tstart; /* start time of observation */
+ e_time tstop; /* end time of observation */
+ short int c1; /* photcode for color term */
+ short int c2; /* photcode for color term */
+ short int photcode; /* appropriate photcode */
+ char label[66];
+} OldPhotPars; /* 100 bytes */
+
+/* Instant Zeropoint structure */
+typedef struct {
+ float ZP; /* measured zero point */
+ float ZPo; /* expected zero point */
+ float dZP; /* error on zero point */
+ float K; /* airmass term used */
+ float X; /* color term used */
+ e_time tstart; /* start time of observation */
+ e_time tstop; /* end time of observation */
+ short int c1; /* photcode for color term */
+ short int c2; /* photcode for color term */
+ short int photcode; /* appropriate photcode */
+ char label[64];
+ short int refcode; /* reference photometry system */
+ int Ntime; /* number of measurement times */
+ int Nmeas; /* number of data points */
+} PhotPars; /* 108 bytes */
+
+/* average data as stored in the LONEOS database */
+typedef struct {
+ float R; /* RA in decimal degrees */
+ float D; /* DEC in decimal degrees */
+ short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */
+ unsigned short int Nm, Nn; /* number of measurements, missing */
+ short int Xp; /* scatter in 1/100 arcsec (-327.67 to +327.67 valid range) */
+ short int Xm; /* chisq for primary mag [1000*value] */
+ unsigned short int code; /* an ID code (ie, star, ghost, satelite, etc) */
+ signed int offset; /* offset to first measurement */
+ signed int missing; /* offset to first missing obs */
+ short int dM; /* formal error on pri mag [1000*log(value)] */
+ short int Xg; /* 'best' chisq value */
+} Average; /* 32 bytes / Average */
+
+/* average data as stored in the LONEOS database -- must be translated on load */
+typedef struct {
+ float R; /* RA in decimal degrees */
+ float D; /* DEC in decimal degrees */
+ short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */
+ unsigned short int Nm, Nn; /* number of measurements, missing */
+ short int Xp; /* scatter in 1/100 arcsec (-327.67 to +327.67 valid range) */
+ short int Xm; /* 1000*log(chisq) for magnitude measurement */
+ unsigned short int code; /* an ID code (ie, star, ghost, satelite, etc) */
+ signed int offset; /* offset to first measurement */
+ signed int missing; /* offset to first missing obs */
+} OldAverage; /* 28 bytes / Average */
+
+/* data for individual measurement for each star */
+/* OLD LONEOS VERSION. NOT BYTE COMPATIBLE WITH NEW VERSION
+ (can be loaded with translation, losing the exptime,
+ and limited the number of Naverage to 0xffffff - in load_catalog) */
+typedef struct {
+ short int dR, dD; /* 1/100 of arcsec (-327.67 to +327.67 valid range) */
+ short int M; /* thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Mcal; /* image cal mag, thousandths of mag (-32.767 to 32.767 valid range) */
+ unsigned char dM; /* thousandths of mag (0.000 to 0.255 valid range) */
+ char dophot; /* dophot type (1-9) */
+ unsigned short int source; /* code to identify photometry source */
+ unsigned int t; /* time in seconds (0 - 143 years valid range) */
+ unsigned int average; /* reference to corresponding average entry, upper byte of Measure.average stores flags */
+} OldMeasure; /* 20 bytes / Measure */
+
+/* data for individual measurement for each star */
+typedef struct {
+ short int dR, dD; /* 1/100 of arcsec (-327.67 to +327.67 valid range) */
+ short int M; /* thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Mcal; /* image cal mag, thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Mgal; /* 'galaxy' mag, thousandths of mag (-32.767 to 32.767 valid range) */
+ short int airmass; /* (airmass - 1), thousandths of mag (-32.767 to 32.767 valid range) */
+ short int FWx; /* object fwhm major axis - 1/100 of arcsec */
+ unsigned char dM; /* thousandths of mag (0.000 to 0.255 valid range) */
+ unsigned char fwy, theta; /* object fwhm minor/major axes ratio, angle wrt ccd X dir (0xff/360 deg precision) */
+ char dophot; /* dophot type (1-9) */
+ unsigned short int source; /* code to identify photometry source */
+ unsigned int t; /* time in seconds (0 - 143 years valid range) */
+ unsigned int averef; /* reference to corresponding average entry, upper byte of Measure.average stores flags */
+ short int dt; /* exposure time in units of 2500*log(exptime) (millimags!) */
+ unsigned short int flags; /* 16 bit-flags for various uses */
+} Measure; /* 32 bytes / Measure */
+
+# if (0)
+/* data for individual measurement for each star */
+typedef struct {
+ short int dR, dD; /* 1/100 of arcsec (-327.67 to +327.67 valid range) */
+ short int M; /* thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Mcal; /* image cal mag, thousandths of mag (-32.767 to 32.767 valid range) */
+ unsigned char dM; /* thousandths of mag (0.000 to 0.255 valid range) */
+ char dophot; /* dophot type (1-9) */
+ unsigned short int source; /* code to identify photometry source */
+ unsigned int t; /* time in seconds (0 - 143 years valid range) */
+ unsigned int averef; /* reference to corresponding average entry, upper byte of Measure.average stores flags */
+ short int dt; /* exposure time in units of 2500*log(exptime) (millimags!) */
+ unsigned short int flags; /* 16 bit-flags for various uses */
+} Measure; /* 24 bytes / Measure */
+# endif
+
+/* data on missing stars */
+typedef struct {
+ unsigned int t; /* time in seconds (0 - 143 years valid range) */
+} Missing; /* 4 byte / Missing */
+
+/* data on missing stars */
+typedef struct {
+ short int M; /* other mags - thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Xm; /* chisq on alt mag */
+ short int dM; /* scatter on alt mag */
+} SecFilt; /* 6 byte / SecFilt */
+
+/* data on missing stars */
+typedef struct {
+ short int M; /* other mags - thousandths of mag (-32.767 to 32.767 valid range) */
+ short int Xm; /* chisq on alt mag */
+} OldSecFilt; /* 4 byte / SecFilt */
+
+/* a catalog contains this data */
+typedef struct {
+ /* data in the file: */
+ Average *average;
+ Measure *measure;
+ Missing *missing;
+ SecFilt *secfilt;
+ Header header;
+ int Naverage, Nmeasure, Nmissing, Nsecfilt;
+ /* description of the catalog file: */
+ char *filename;
+ int lockmode;
+ FILE *f;
+ /* pointers for data manipulation */
+ int *found;
+ int *image;
+ int *mosaic;
+ float *X;
+ float *Y;
+} Catalog;
+
+/* the standard astrometric coordinates structure */
+typedef struct {
+ double crval1, crval2;
+ float crpix1, crpix2;
+ float cdelt1, cdelt2;
+ float pc1_1, pc1_2;
+ float pc2_1, pc2_2;
+ float polyterms[7][2];
+ char ctype[15];
+ char Npolyterms;
+} Coords; /* 120 bytes / Coords */
+
+typedef struct {
+ unsigned short int code; /* code number (stored in Measure.source) */
+ char name[32]; /* name for filter combination */
+ char type; /* primary, secondary, dependent, reference */
+ short int C, dC, dX; /* primary phot calibration terms (millimags) */
+ float K; /* secondary phot calibration terms (millimags) */
+ int c1, c2; /* color is average.M[c1] - average.M[c2] */
+ int equiv; /* this dependent filter is equivalent to equiv primary/secondary */
+ int Nc; /* number of color terms */
+ float X[4]; /* color terms X[0]*mc + X[1]*mc^2 + X[2]*mc^3, etc */
+} PhotCode;
+
+typedef struct {
+ int Ncode;
+ int Nsecfilt;
+ int hashcode[0x10000];
+ int hashNsec[0x10000];
+ PhotCode *code;
+} PhotCodeData;
+
+/* image structure for binary data storage of image data */
+typedef struct {
+ Coords coords; /* 120 bytes */
+ unsigned int tzero; /* readout time row 0 in sec (0 - 142 years valid range) */
+ unsigned int nstar; /* number of stars on image */
+ short int secz; /* thousanths of airmass (valid range -32.000 -- 32.000) */
+ short int NX, NY; /* dimensions of image */
+ short int apmifit, dapmifit; /* aperture correction and error in thousandths of mag */
+ short int source; /* identifier for CCD (each ever used will have a unique letter) */
+ short int Mcal; /* thousandths of mag (-32.000 -- 32.000 valid range) */
+ short int dMcal; /* thousandths of mag (-32.000 -- 32.000 valid range) */
+ short int Xm; /* 10*log(image chi-square) */
+ char name[32]; /* name of original image */
+ unsigned char detection_limit; /* tenths of mag (0.0 - 25.6 valid range) */
+ unsigned char saturation_limit; /* tenths of mag (0.0 - 25.6 valid range) */
+ unsigned char cerror; /* astrometric error: 1/50 of arcsec (0 -- 5.12 valid range) */
+ unsigned char fwhm_x, fwhm_y; /* PSF terms in 25*arcsec (valid range 0.0 -- 10.2" ") */
+ unsigned char trate; /* 10000 * scan rate in sec/pix (0 -- 0.0256 valid range,
+ typically 0.0146. this is used only to determine the time
+ of the observation, not to find the coordinates. 1 byte
+ gives 0.11 sec accuracy */
+ float exptime; /* exposure time, seconds */
+ char code; /* flag to mark an image as bad or whatever */
+ unsigned char ccdnum; /* mosaic CCD ID number (IMAGEID) */
+ char dummy[20]; /* extra space for the future */
+ short int order; /* number of terms used for Mrel */
+ short int Mx, My;
+ short int Mxx, Mxy, Myy;
+ short int Mxxx, Mxxy, Mxyy, Myyy;
+ short int Mxxxx, Mxxxy, Mxxyy, Mxyyy, Myyyy;
+} Image; /* 240 bytes / Image */
+/* Image needs to be aligned with 8-byte boundaries on Suns */
+
+/**** prototypes ****/
+
+int lock_catalog (Catalog *catalog, int lockmode);
+int unlock_catalog (Catalog *catalog);
+int load_catalog (Catalog *catalog, char mode, int VERBOSE);
+int save_catalog (Catalog *catalog, char VERBOSE);
+Measure *FixOldMeasure (OldMeasure *in, int Nvalues);
+Average *FixOldAverage (OldAverage *in, int Nvalues);
+SecFilt *FixOldSecFilt (OldSecFilt *in, int Nvalues);
+
+int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords);
+int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords);
+int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords);
+int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords);
+int GetCoords (Coords *coords, Header *header);
+short int putMi (double value);
+double getMi (short int value);
+void returnMcal (Image *image, double *c);
+void assignMcal (Image *image, double *c, int order);
+double applyMcal (Image *image, double x, double y);
+double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order);
+
+PhotCode *GetPhotcodebyName (char *name);
+PhotCode *GetPhotcodeEquivbyName (char *name);
+PhotCode *GetPhotcodebyCode (int code);
+PhotCode *GetPhotcodeEquivbyCode (int code);
+char *GetPhotcodeNamebyCode (int code);
+
+double PhotInst (Measure *measure);
+double PhotAbs (Measure *measure);
+double PhotCat (Measure *measure);
+double PhotSys (Measure *measure, Average *average, SecFilt *secfilt);
+double PhotRel (Measure *measure, Average *average, SecFilt *secfilt);
+double PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code);
+double PhotAve (PhotCode *code, Average *average, SecFilt *secfilt);
+double PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure);
+double PhotXm (PhotCode *code, Average *average, SecFilt *secfilt);
+double PhotdM (PhotCode *code, Average *average, SecFilt *secfilt);
+
+short iPhotInst (Measure *measure);
+short iPhotAbs (Measure *measure);
+short iPhotCat (Measure *measure);
+short iPhotSys (Measure *measure, Average *average, SecFilt *secfilt);
+short iPhotRel (Measure *measure, Average *average, SecFilt *secfilt);
+short iPhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code);
+short iPhotAve (PhotCode *code, Average *average, SecFilt *secfilt);
+short iPhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure);
+short iPhotXm (PhotCode *code, Average *average, SecFilt *secfilt);
+short iPhotdM (PhotCode *code, Average *average, SecFilt *secfilt);
+
+short iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code);
+int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color);
+
+int LoadPhotcodes (char *filename);
+int GetPhotcodeCodebyName (char *name);
+int GetPhotcodeEquivCodebyName (char *name);
+int GetPhotcodeEquivCodebyCode (int code);
+int GetPhotcodeNsec (int code);
+int GetPhotcodeNsecfilt ();
+void SetZeroPoint (double ZP);
+int *GetPhotcodeEquivList (int code, int *nlist);
+void ParseColorTerms (char *terms, float *X, int *N);
+
+void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch);
+
+int get_image_type (char *name);
+char *get_type_name (int type);
+int get_image_mode (char *name);
+char *get_mode_name (int mode);
+
+# endif
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/include/ohana.h
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/include/ohana.h (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/include/ohana.h (revision 21775)
@@ -0,0 +1,151 @@
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+
+/* OHANA included stuff */
+# include
+
+# ifndef OHANA_H
+# define OHANA_H
+
+extern double hypot();
+
+# define TRUE (1)
+# define FALSE (0)
+# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# define SQR(X) (double) (((double)(X))*((double)(X)))
+# define SQ(X) (double) (((double)(X))*((double)(X)))
+# define MIN(X,Y) ((X) < (Y) ? (X) : (Y))
+# define MAX(X,Y) ((X) > (Y) ? (X) : (Y))
+# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;}
+
+# define LCK_UNLOCK 0 /* file is unlocked */
+# define LCK_ACCESS 1 /* can't get access to file */
+# define LCK_TIMEOUT 2 /* timeout setting lock */
+# define LCK_HARDLCK 3 /* error setting hard lockfile */
+# define LCK_HUNLOCK 4 /* error clearing hard lockfile */
+
+# define LCK_EMPTY 5 /* locked file is empty */
+# define LCK_FULL 6 /* locked file is not empty */
+# define LCK_UNKNOWN 7 /* can't stat file to get size */
+
+# define LCK_SOFT 10 /* soft lock */
+# define LCK_HARD 11 /* hard lock */
+# define LCK_XCLD 12 /* exclusive soft lock */
+
+# ifndef M_PI
+# define M_PI 3.14159265358979323846264
+# endif
+
+int setlockfile2 (char *filename, double timeout, int type, int *state);
+int clearlockfile2 (char *filename, int fd, int type, int *state);
+
+FILE *fsetlockfile (char *filename, double timeout, int type, int *state);
+int fclearlockfile (char *filename, FILE *f, int type, int *state);
+
+# ifndef PROTO
+# define PROTO(A) A
+# endif
+
+/*
+# else
+# ifndef PROTO
+# define PROTO(A) ()
+# endif
+# include
+# endif
+*/
+
+# ifndef ALLOCATE
+# define ALLOCATE(X,T,S) \
+ X = (T *) malloc ((unsigned) (MAX(((S)*sizeof(T)),1)));\
+ if (X == NULL) { \
+ fprintf(stderr,"failed malloc at __LINE__ in __FILE__\n");\
+ exit (10); }
+# define REALLOCATE(X,T,S) \
+ X = (T *) realloc(X,(unsigned) (MAX(((S)*sizeof(T)),1))); \
+ if (X == NULL) { \
+ fprintf(stderr,"failed realloc at __LINE__ in __FILE__\n"); \
+ exit (10); }
+# define CHECK_REALLOCATE(X,T,S,N,D) \
+ if ((N) >= (S)) { \
+ S += D; \
+ X = (T *) realloc (X,(unsigned) (MAX(((S)*sizeof(T)),1))); \
+ if (X == NULL) { \
+ fprintf(stderr,"failed to realloc X __LINE__ in __FILE__\n"); \
+ exit (10); } }
+# endif /* ALLOCATE */
+
+# ifndef FOPEN
+# define FOPEN(F,NAME) \
+ F = fopen (NAME, "r"); \
+ if (F == NULL) { \
+ fprintf (stderr, "failed to open %s\n", NAME); \
+ exit (0); \
+ }
+# endif /* FOPEN */
+
+# define DEG_RAD 57.295779513082322
+# define RAD_DEG 0.017453292519943
+
+# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec))
+
+/* in string.c */
+int stripwhite (char *string);
+int strnumcmp (char *str1, char *str2);
+char *strcreate (char *string);
+char *strncreate (char *string, int n);
+int scan_line (FILE *f, char *line);
+int dparse (double *X, int NX, char *line);
+int fparse (float *X, int NX, char *line);
+int get_argument (int argc, char **argv, char *arg);
+int remove_argument (int N, int *argc, char **argv);
+
+int Fread (void *ptr, int size, int nitems, FILE *f, char *type);
+int Fwrite (void *ptr, int size, int nitems, FILE *f, char *type);
+int ByteSwap (char *ptr, int size, int nitems, char *type);
+int ConvertStruct (char *buffer, int size, int Nbytes, char *type);
+
+char *pathname (char *name);
+char *filebasename (char *name);
+char *filerootname (char *name);
+char *findexec (int argc, char **argv);
+int mkdirhier (char *path);
+
+FILE *fsetlockfile (char *filename, double timeout, int type, int *state);
+int fclearlockfile (char *filename, FILE *f, int type, int *state);
+int setlockfile2 (char *filename, double timeout, int type, int *state);
+int clearlockfile2 (char *filename, int fd, int type, int *state);
+
+char *SelectConfigFile (int *argc, char **argv, char *progname);
+char *LoadConfigFile (char *filename);
+char *ScanConfig (char *config, char *field, char *mode, int N,...);
+char *expandline (char *line, char *config);
+char *fileextname (char *file);
+char *LoadRawConfigFile (char *, int);
+
+/*
+# define F_SETFL 4
+# define O_NONBLOCK 0200000
+# define AF_UNIX 1
+# define SOCK_STREAM 1
+*/
+
+# endif
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/lib/.cvsignore
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/lib/.cvsignore (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/lib/.cvsignore (revision 21775)
@@ -0,0 +1,1 @@
+*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/Fread.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/Fread.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/Fread.c (revision 21775)
@@ -0,0 +1,466 @@
+# 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;
+# define SWAP_BYTE(X) \
+ tmp = byte[X+0]; byte[X+0] = byte[X+1]; byte[X+1] = tmp;
+# define SWAP_WORD(X) \
+ tmp = byte[X+0]; byte[X+0] = byte[X+3]; byte[X+3] = tmp; \
+ tmp = byte[X+1]; byte[X+1] = byte[X+2]; byte[X+2] = tmp;
+# define SWAP_DBLE(X) \
+ tmp = byte[X+0]; byte[X+0] = byte[X+7]; byte[X+7] = tmp; \
+ tmp = byte[X+1]; byte[X+1] = byte[X+6]; byte[X+6] = tmp; \
+ tmp = byte[X+2]; byte[X+2] = byte[X+5]; byte[X+5] = tmp; \
+ tmp = byte[X+3]; byte[X+3] = byte[X+4]; byte[X+4] = tmp;
+# define IMAGE_SIZE 240
+# define OLDAVERAGE_SIZE 28
+# define AVERAGE_SIZE 32
+# define OLDMEASURE_SIZE 20
+# define MEASURE_SIZE 32
+# define MISSING_SIZE 4
+# define OLDSECFILT_SIZE 4
+# define SECFILT_SIZE 6
+# define REGIMAGE_SIZE 360
+# define SPECTRUM_SIZE 216
+# define RUFIMAGE_SIZE 60
+# define OLDDETREG_SIZE 288
+# define DETREG_SIZE 416
+# define PHOTPARS_SIZE 108
+# define OLDPHOTPARS_SIZE 100
+# define TRUE 1
+# define FALSE 0
+
+# ifdef linux
+# define BYTE_SWAP
+# endif
+
+# ifdef sid
+# define BYTE_SWAP
+# endif
+
+/* add other architectures here, ie dec, alpha, etc */
+# ifdef dec
+# define BYTE_SWAP
+# endif
+
+int Fread (char *ptr, int size, int nitems, FILE *f, char *type) {
+
+ int valid, status;
+
+ status = fread (ptr, size, nitems, f);
+
+ valid = ByteSwap (ptr, size, nitems, type);
+
+ if (!valid) return (FALSE);
+ return (status);
+}
+
+int Fwrite (char *ptr, int size, int nitems, FILE *f, char *type) {
+
+ int valid, status;
+
+ valid = ByteSwap (ptr, size, nitems, type);
+
+ if (!valid) return (FALSE);
+
+ status = fwrite (ptr, size, nitems, f);
+ return (status);
+}
+
+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, valid;
+ unsigned char *byte, *byte0, *byte1, *byte2, *byte3, *byte4, *byte5, *byte6, *byte7, tmp;
+
+ valid = FALSE;
+ if (!strcmp (type, "char")) {
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "short")) {
+ byte0 = (char *)ptr;
+ byte1 = (char *)ptr + 1;
+ for (i = 0; i < nitems; i++, byte0 += 2, byte1 += 2) {
+ tmp = *byte0;
+ *byte0 = *byte1;
+ *byte1 = tmp;
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "int") || !strcmp (type, "float")) {
+ byte0 = (char *)ptr;
+ byte1 = (char *)ptr + 1;
+ byte2 = (char *)ptr + 2;
+ byte3 = (char *)ptr + 3;
+ for (i = 0; i < nitems; i++, byte0 += 4, byte1 += 4, byte2 += 4, byte3 += 4) {
+ tmp = *byte0;
+ *byte0 = *byte3;
+ *byte3 = tmp;
+ tmp = *byte1;
+ *byte1 = *byte2;
+ *byte2 = tmp;
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "double")) {
+ byte0 = (char *)ptr;
+ byte1 = (char *)ptr + 1;
+ byte2 = (char *)ptr + 2;
+ byte3 = (char *)ptr + 3;
+ byte4 = (char *)ptr + 4;
+ byte5 = (char *)ptr + 5;
+ byte6 = (char *)ptr + 6;
+ byte7 = (char *)ptr + 7;
+ for (i = 0; i < nitems; i++, byte0 += 8, byte1 += 8, byte2 += 8, byte3 += 8, byte4 += 8, byte5 += 8, byte6 += 8, byte7 += 8) {
+ tmp = *byte0;
+ *byte0 = *byte7;
+ *byte7 = tmp;
+ tmp = *byte1;
+ *byte1 = *byte6;
+ *byte6 = tmp;
+ tmp = *byte1;
+ *byte2 = *byte5;
+ *byte5 = tmp;
+ tmp = *byte1;
+ *byte3 = *byte4;
+ *byte4 = tmp;
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "average")) {
+ if (size != AVERAGE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (Average) %d vs %d\n", size, AVERAGE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += AVERAGE_SIZE) {
+ SWAP_WORD (0); /* R */
+ SWAP_WORD (4); /* D */
+ SWAP_BYTE (8); /* M */
+ SWAP_BYTE (10); /* Nm */
+ SWAP_BYTE (12); /* Nn */
+ SWAP_BYTE (14); /* Xp */
+ SWAP_BYTE (16); /* Xm */
+ SWAP_BYTE (18); /* code */
+ SWAP_WORD (20); /* offset */
+ SWAP_WORD (24); /* missing */
+ SWAP_BYTE (28); /* dM */
+ SWAP_BYTE (30); /* Xg */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "oldaverage")) {
+ if (size != OLDAVERAGE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (OldAverage) %d vs %d\n", size, OLDAVERAGE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += OLDAVERAGE_SIZE) {
+ SWAP_WORD (0); /* R */
+ SWAP_WORD (4); /* D */
+ SWAP_BYTE (8); /* M */
+ SWAP_BYTE (10); /* Nm */
+ SWAP_BYTE (12); /* Nn */
+ SWAP_BYTE (14); /* Xp */
+ SWAP_BYTE (16); /* Xm */
+ SWAP_BYTE (18); /* code */
+ SWAP_WORD (20); /* offset */
+ SWAP_WORD (24); /* missing */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "regimage")) {
+ if (size != REGIMAGE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (RegImage) %d vs %d\n", size, REGIMAGE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += REGIMAGE_SIZE) {
+ SWAP_WORD (284); /* exptime */
+ SWAP_WORD (288); /* airmass */
+ SWAP_WORD (292); /* sky */
+ SWAP_WORD (296); /* bias */
+ SWAP_WORD (300); /* fwhm */
+ SWAP_WORD (304); /* telfocus */
+ SWAP_WORD (308); /* xprobe */
+ SWAP_WORD (312); /* yprobe */
+ SWAP_WORD (316); /* zprobe */
+ SWAP_WORD (320); /* dettemp */
+ SWAP_WORD (324); /* teltemp[0] */
+ SWAP_WORD (328); /* teltemp[1] */
+ SWAP_WORD (332); /* teltemp[2] */
+ SWAP_WORD (336); /* teltemp[3] */
+ SWAP_WORD (340); /* rotangle */
+ SWAP_WORD (344); /* ra */
+ SWAP_WORD (348); /* dec */
+ SWAP_WORD (352); /* obstime */
+ SWAP_WORD (356); /* regtime */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "spectrum")) {
+ if (size != SPECTRUM_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (Spectrum) %d vs %d\n", size, SPECTRUM_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += SPECTRUM_SIZE) {
+ SWAP_WORD (0); /* ra */
+ SWAP_WORD (4); /* dec */
+ SWAP_WORD (8); /* exptime */
+ SWAP_WORD (12); /* airmass */
+ SWAP_WORD (16); /* Ws */
+ SWAP_WORD (20); /* We */
+ SWAP_WORD (24); /* dW */
+ SWAP_WORD (28); /* Nspec */
+ SWAP_WORD (32); /* obstime */
+ SWAP_WORD (36); /* regtime */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "rufimage")) {
+ if (size != RUFIMAGE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (RufImage) %d vs %d\n", size, RUFIMAGE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += RUFIMAGE_SIZE) {
+ SWAP_WORD (0); /* sky */
+ SWAP_WORD (4); /* bias */
+ SWAP_WORD (8); /* fwhm */
+ SWAP_WORD (12); /* exptime */
+ SWAP_WORD (16); /* airmass */
+ SWAP_WORD (20); /* obstime */
+ SWAP_BYTE (24); /* ccdnum */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "olddetreg")) {
+ if (size != OLDDETREG_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (DetReg) %d vs %d\n", size, OLDDETREG_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += OLDDETREG_SIZE) {
+ SWAP_WORD (0); /* tstart */
+ SWAP_WORD (4); /* tstop */
+ SWAP_WORD (8); /* treg */
+ SWAP_WORD (12); /* sigma */
+ SWAP_WORD (16); /* clipsigma */
+ SWAP_WORD (20); /* type */
+ SWAP_WORD (24); /* filter */
+ SWAP_WORD (28); /* ccd */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "detreg")) {
+ if (size != DETREG_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (DetReg) %d vs %d\n", size, DETREG_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += DETREG_SIZE) {
+ SWAP_WORD (0); /* tstart */
+ SWAP_WORD (4); /* tstop */
+ SWAP_WORD (8); /* treg */
+ SWAP_WORD (12); /* exptime */
+ SWAP_WORD (16); /* type */
+ SWAP_WORD (20); /* filter */
+ SWAP_WORD (24); /* ccd */
+ SWAP_WORD (28); /* Nentry */
+ SWAP_WORD (32); /* Norder */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "oldphotpars")) {
+ if (size != OLDPHOTPARS_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (OldPhotPars) %d vs %d\n", size, OLDPHOTPARS_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += OLDPHOTPARS_SIZE) {
+ SWAP_WORD (0); /* ZP */
+ SWAP_WORD (4); /* ZPo */
+ SWAP_WORD (8); /* dZP */
+ SWAP_WORD (12); /* K */
+ SWAP_WORD (16); /* A */
+ SWAP_WORD (20); /* tstart */
+ SWAP_WORD (24); /* tstop */
+ SWAP_BYTE (28); /* c1 */
+ SWAP_BYTE (30); /* c2 */
+ SWAP_BYTE (32); /* photcode */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "photpars")) {
+ if (size != PHOTPARS_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (PhotPars) %d vs %d\n", size, PHOTPARS_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += PHOTPARS_SIZE) {
+ SWAP_WORD (0); /* ZP */
+ SWAP_WORD (4); /* ZPo */
+ SWAP_WORD (8); /* dZP */
+ SWAP_WORD (12); /* K */
+ SWAP_WORD (16); /* A */
+ SWAP_WORD (20); /* tstart */
+ SWAP_WORD (24); /* tstop */
+ SWAP_BYTE (28); /* c1 */
+ SWAP_BYTE (30); /* c2 */
+ SWAP_BYTE (32); /* photcode */
+ SWAP_WORD (100); /* Ntime */
+ SWAP_WORD (104); /* Nmeas */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "oldmeasure")) {
+ if (size != OLDMEASURE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (OldMeasure) %d vs %d\n", size, OLDMEASURE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += OLDMEASURE_SIZE) {
+ SWAP_BYTE (0); /* dR */
+ SWAP_BYTE (2); /* dD */
+ SWAP_BYTE (4); /* M */
+ SWAP_BYTE (6); /* Mcal */
+ SWAP_BYTE (10); /* source */
+ SWAP_WORD (12); /* t */
+ SWAP_WORD (16); /* average */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "measure")) {
+ if (size != MEASURE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (Measure) %d vs %d\n", size, MEASURE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += MEASURE_SIZE) {
+ SWAP_BYTE (0); /* dR */
+ SWAP_BYTE (2); /* dD */
+ SWAP_BYTE (4); /* M */
+ SWAP_BYTE (6); /* Mcal */
+ SWAP_BYTE (8); /* Mgal */
+ SWAP_BYTE (10); /* Map */
+ SWAP_BYTE (12); /* FWx */
+ /* 14, 15, 16, 17 - char */
+ SWAP_BYTE (18); /* source */
+ SWAP_WORD (20); /* t */
+ SWAP_WORD (24); /* averef */
+ SWAP_BYTE (28); /* dt */
+ SWAP_BYTE (30); /* flags */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "secfilt")) {
+ if (size != SECFILT_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (SecFilt) %d vs %d\n", size, SECFILT_SIZE);
+ return (0);
+ }
+ byte = (char *) ptr;
+ for (i = 0; i < nitems; i++, byte += SECFILT_SIZE) {
+ SWAP_BYTE (0); /* M */
+ SWAP_BYTE (2); /* Xm */
+ SWAP_BYTE (4); /* dM */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "oldsecfilt")) {
+ if (size != OLDSECFILT_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (OldSecFilt) %d vs %d\n", size, OLDSECFILT_SIZE);
+ return (0);
+ }
+ byte = (char *) ptr;
+ for (i = 0; i < nitems; i++, byte += OLDSECFILT_SIZE) {
+ SWAP_BYTE (0); /* M */
+ SWAP_BYTE (2); /* Xm */
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "missing")) {
+ if (size != MISSING_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (Missing) %d vs %d\n", size, MISSING_SIZE);
+ return (0);
+ }
+ byte0 = (char *)ptr;
+ byte1 = (char *)ptr + 1;
+ byte2 = (char *)ptr + 2;
+ byte3 = (char *)ptr + 3;
+ for (i = 0; i < nitems; i++, byte0 += 4, byte1 += 4, byte2 += 4, byte3 += 4) {
+ tmp = *byte0;
+ *byte0 = *byte3;
+ *byte3 = tmp;
+ tmp = *byte1;
+ *byte1 = *byte2;
+ *byte2 = tmp;
+ }
+ valid = TRUE;
+ }
+
+ if (!strcmp (type, "image")) {
+ if (size != IMAGE_SIZE) {
+ fprintf (stderr, "mismatch in type sizes (Image) %d vs %d\n", size, IMAGE_SIZE);
+ return (0);
+ }
+ byte = (char *)ptr;
+ for (i = 0; i < nitems; i++, byte += IMAGE_SIZE) {
+ SWAP_DBLE (0); /* coords.crval1 */
+ SWAP_DBLE (8); /* coords.crval2 */
+ for (j = 16; j < 104; j+=4) { /* coords.crpix, delt, pc_ij, polyterms[7][2] */
+ SWAP_WORD (j);
+ }
+ SWAP_WORD (120); /* tzero */
+ SWAP_WORD (124); /* nstar */
+ for (j = 128; j < 146; j+=2) { /* sec z, NX, NY, apmifit, dapmifit, source, Mcal, dMcal, Xm */
+ SWAP_BYTE (j);
+ }
+ SWAP_WORD (184); /* exptime */
+ SWAP_WORD (210); /* order */
+ for (j = 212; j < 240; j+=2) {
+ SWAP_BYTE (j); /* Mx - Myyy */
+ }
+ }
+ valid = TRUE;
+ }
+
+ if (!valid)
+ fprintf (stderr, "unknown type %s\n", type);
+
+ return (valid);
+
+# else
+
+ return (TRUE);
+
+# endif
+
+}
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/LoadPhotcodes.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/LoadPhotcodes.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/LoadPhotcodes.c (revision 21775)
@@ -0,0 +1,793 @@
+# include
+# include
+
+# define NCTERMS 4
+
+static PhotCodeData *photcodes = NULL;
+static double ZERO_POINT;
+static short int iZERO_POINT;
+static short int Nseclist[0x10000];
+
+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);
+}
+
+/******** 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);
+}
+
+/******** internal photometry conversions (keeps values in short int millimags) *********/
+short iPhotInst (Measure *measure) {
+
+ int Np;
+ short M;
+
+ Np = photcodes[0].hashcode[measure[0].source];
+ if (Np == -1) return (NO_MAG);
+
+ if (photcodes[0].code[Np].type == PHOT_REF) {
+ M = measure[0].M;
+ return (M);
+ }
+
+ M = measure[0].M - measure[0].dt - iZERO_POINT;
+
+ return (M);
+
+}
+
+short iPhotCat (Measure *measure) {
+
+ int Np;
+ short Mcat;
+ PhotCode *code;
+
+ Np = photcodes[0].hashcode[measure[0].source];
+ if (Np == -1) return (NO_MAG);
+
+ if (photcodes[0].code[Np].type == PHOT_REF) {
+ Mcat = measure[0].M;
+ return (Mcat);
+ }
+ code = &photcodes[0].code[Np];
+ Mcat = measure[0].M - iZERO_POINT + code[0].K*(measure[0].airmass - 1000) + code[0].C;
+
+ return (Mcat);
+}
+
+short iPhotSys (Measure *measure, Average *average, SecFilt *secfilt) {
+
+ int i, Np;
+ short Mcat, Mcol, Msys, mc;
+ double Mc;
+ PhotCode *code;
+
+ Np = photcodes[0].hashcode[measure[0].source];
+ if (Np == -1) return (NO_MAG);
+
+ if (photcodes[0].code[Np].type == PHOT_REF) {
+ Msys = measure[0].M;
+ return (Msys);
+ }
+ code = &photcodes[0].code[Np];
+ Mcat = measure[0].M - iZERO_POINT + code[0].K*(measure[0].airmass - 1000.0) + code[0].C;
+
+ /* for DEP, color must be made of PRI/SEC */
+ mc = iPhotColor (average, secfilt, NULL, code);
+ if (mc == NO_MAG) return (Mcat);
+ mc = mc - code[0].dX;
+
+ Mc = mc;
+ Mcol = 0;
+ for (i = 0; i < code[0].Nc; i++) {
+ Mcol += code[0].X[i]*Mc;
+ Mc *= mc * 0.001; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
+ }
+ Msys = Mcat + Mcol;
+ return (Msys);
+}
+
+short iPhotRel (Measure *measure, Average *average, SecFilt *secfilt) {
+
+ int i, Np;
+ short Mcat, Mcol, Mrel, mc;
+ double Mc;
+ PhotCode *code;
+
+ Np = photcodes[0].hashcode[measure[0].source];
+ if (Np == -1) return (NO_MAG);
+
+ if (photcodes[0].code[Np].type == PHOT_REF) {
+ Mcat = measure[0].M;
+ return (Mcat);
+ }
+ code = &photcodes[0].code[Np];
+ Mrel = measure[0].M - iZERO_POINT + code[0].K*(measure[0].airmass - 1000.0) + code[0].C - measure[0].Mcal;
+
+ /* for DEP, color must be made of PRI/SEC */
+ mc = iPhotColor (average, secfilt, NULL, code);
+ if (mc == NO_MAG) return (Mrel);
+ mc = mc - code[0].dX;
+
+ Mc = mc;
+ Mcol = 0;
+ for (i = 0; i < code[0].Nc; i++) {
+ Mcol += code[0].X[i]*Mc;
+ Mc *= mc * 0.001; /* 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 */
+short iPhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {
+
+ int i, Np;
+ short Mcal, Mrel, Mcol, mc;
+ double 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);
+ }
+
+ if (photcodes[0].code[Np].type == PHOT_REF) {
+ Mrel = thisone[0].M;
+ return (Mrel);
+ }
+ if (code[0].code != photcodes[0].code[Np].equiv) {
+ return (NO_MAG);
+ }
+
+ Mcal = iPhotRel (thisone, average, secfilt) + code[0].C;
+
+ mc = iPhotColor (average, secfilt, measure, code);
+ if (mc == NO_MAG) return (Mcal);
+ mc = mc - code[0].dX;
+
+ Mc = mc;
+ Mcol = 0;
+ for (i = 0; i < code[0].Nc; i++) {
+ Mcol += code[0].X[i]*Mc;
+ Mc *= mc * 0.001; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
+ }
+ Mcal += Mcol;
+ return (Mcal);
+}
+
+/* color term may not use DEP magnitude */
+short iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {
+
+ int i, Ns1, Ns2, Ns;
+ short m1, m2, mc;
+ PhotCode *color;
+
+ m1 = m2 = NO_MAG;
+
+ if (measure == NULL) {
+ Ns1 = photcodes[0].hashNsec[code[0].c1];
+ Ns2 = photcodes[0].hashNsec[code[0].c2];
+
+ m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M;
+ m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M;
+ mc = ((m1 == NO_MAG) || (m2 == NO_MAG)) ? NO_MAG : (m1 - m2);
+ return (mc);
+ }
+
+ /* find magnitude matching first color term */
+ color = GetPhotcodebyCode (code[0].c1);
+ if (color == NULL) return (NO_MAG);
+ if (color[0].type == PHOT_REF) {
+ for (i = 0; (i < average[0].Nm) && (m1 == NO_MAG); i++) {
+ if (measure[i].source == color[0].code) {
+ m1 = measure[i].M;
+ }
+ }
+ } else {
+ Ns = photcodes[0].hashNsec[color[0].code];
+ m1 = (Ns == -1) ? average[0].M : secfilt[Ns].M;
+ }
+
+ /* find magnitude matching second color term */
+ color = GetPhotcodebyCode (code[0].c2);
+ if (color == NULL) return (NO_MAG);
+ if (color[0].type == PHOT_REF) {
+ for (i = 0; (i < average[0].Nm) && (m2 == NO_MAG); i++) {
+ if (measure[i].source == color[0].code) {
+ m2 = measure[i].M;
+ }
+ }
+ } else {
+ Ns = photcodes[0].hashNsec[color[0].code];
+ m2 = (Ns == -1) ? average[0].M : secfilt[Ns].M;
+ }
+ mc = ((m1 == NO_MAG) || (m2 == NO_MAG)) ? NO_MAG : (m1 - m2);
+ return (mc);
+}
+
+/* return calibrated magnitude from average/secfilt for given photcode */
+short iPhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {
+
+ int i, Ns;
+ short Mave, Mref, Mcol, mc;
+ double Mc;
+
+ Ns = photcodes[0].hashNsec[code[0].code];
+ Mave = (Ns == -1) ? average[0].M : secfilt[Ns].M;
+ Mref = Mave + code[0].C;
+
+ mc = iPhotColor (average, secfilt, measure, code);
+ if (mc == NO_MAG) return (Mref);
+ mc = mc - code[0].dX;
+
+ Mc = mc;
+ Mcol = 0;
+ for (i = 0; i < code[0].Nc; i++) {
+ Mcol += code[0].X[i]*Mc;
+ Mc *= mc * 0.001; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
+ }
+ Mref += Mcol;
+ return (Mref);
+}
+
+/***/
+short iPhotAve (PhotCode *code, Average *average, SecFilt *secfilt) {
+
+ int Ns;
+ short Mave;
+
+ Ns = photcodes[0].hashNsec[code[0].code];
+ Mave = (Ns == -1) ? average[0].M : secfilt[Ns].M;
+ return (Mave);
+}
+
+/*** note that this is NOT a wrapper around iPhotdM ***/
+short iPhotdM (PhotCode *code, Average *average, SecFilt *secfilt) {
+
+ int Ns;
+ short dM;
+
+ Ns = photcodes[0].hashNsec[code[0].code];
+ dM = (Ns == -1) ? average[0].dM : secfilt[Ns].dM;
+ return (dM);
+}
+
+/*** note that this is NOT a wrapper around iPhotXm ***/
+short iPhotXm (PhotCode *code, Average *average, SecFilt *secfilt) {
+
+ int Ns;
+ short Xm;
+
+ Ns = photcodes[0].hashNsec[code[0].code];
+ Xm = (Ns == -1) ? average[0].Xm : secfilt[Ns].Xm;
+ 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 = 0.001*measure[i].M;
+ goto filter1;
+ }
+ }
+ return (FALSE);
+ } else {
+ Ns = photcodes[0].hashNsec[code[0].code];
+ M1 = (Ns == -1) ? 0.001*average[0].M : 0.001*secfilt[Ns].M;
+ }
+
+ 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 = 0.001*measure[i].M;
+ goto filter2;
+ }
+ }
+ return (FALSE);
+ } else {
+ Ns = photcodes[0].hashNsec[code[0].code];
+ M2 = (Ns == -1) ? 0.001*average[0].M : 0.001*secfilt[Ns].M;
+ }
+
+ filter2:
+
+ dM = M1 - M2;
+ *color = dM;
+
+ return (TRUE);
+}
+
+/* photcode table should have the following format:
+
+# code name type zero airmass offset c1 c2 slope primary
+ 1 B pri 24.0 0.15 - - - - - -
+ 2 B pri 24.0 0.15 - - - - - -
+ 3 B1 sec 22.5 0.18 0.15 1 2 0.10 0.50 1
+ 1000 USNO_B ref - - - - - - - -
+
+*/
+
+
+/*
+ Nc1 = photcodes[0].code[Np].c1;
+ Ns1 = photcodes[0].hashNsec[Nc1];
+
+ Nc2 = photcodes[0].code[Np].c2;
+ Ns2 = photcodes[0].hashNsec[Nc2];
+
+ Xlam = photcodes[0].code[Np].X[0];
+ Klam = photcodes[0].code[Np].K;
+
+ m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M;
+ m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M;
+*/
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/config.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/config.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/config.c (revision 21775)
@@ -0,0 +1,457 @@
+# include
+# ifndef whitespace
+# define whitespace(c) (((c) == ' ') || ((c) == '\t'))
+# endif
+
+# define D_NBYTES 4096
+
+static char *ConfigVariable = (char *) NULL;
+static int NDefineVariable;
+static char **DefineVariable;
+static char **DefineValue;
+
+char *SelectConfigFile (int *argc, char **argv, char *progname) {
+
+ /*
+ config file selection rules (first ones override later ones):
+
+ 1) -c Filename : use Filename
+ 2) PROGNAME : use environment variable as config file
+ 3) progname.rc : use alternate name in local dir as config file
+ 4) .prognamerc : use rc file in local dir as config file
+ 5) ~/.prognamerc : use rc file in homedir as config file
+
+ special variable definitions:
+ 1) -C WORD : set CONFIG variable to WORD
+ 2) -D NAME WORD : set NAME variable to WORD
+ -D overrides variables in param file
+
+ these command-line options are removed and a complete arg list left behind
+ */
+
+ char *filename, *find, *home;
+ struct stat filestat;
+ uid_t uid;
+ gid_t gid;
+ int i, N, NDEF, status;
+
+ /* first look for -C CONFIG variable */
+ if ((N = get_argument (*argc, argv, "-C"))) {
+ remove_argument (N, argc, argv);
+ if (ConfigVariable != (char *) NULL) free (ConfigVariable);
+ ConfigVariable = strcreate (argv[N]);
+ remove_argument (N, argc, argv);
+ }
+
+ /* next look for -D NAME WORD variables */
+ NDEF = 10;
+ ALLOCATE (DefineVariable, char *, NDEF);
+ ALLOCATE (DefineValue, char *, NDEF);
+ for (i = 0; (N = get_argument (*argc, argv, "-D")); i++) {
+ remove_argument (N, argc, argv);
+ DefineVariable[i] = strcreate (argv[N]);
+ remove_argument (N, argc, argv);
+ DefineValue[i] = strcreate (argv[N]);
+ remove_argument (N, argc, argv);
+ if (i == NDEF - 1) {
+ NDEF += 10;
+ REALLOCATE (DefineVariable, char *, NDEF);
+ REALLOCATE (DefineValue, char *, NDEF);
+ }
+ }
+ NDefineVariable = i;
+ REALLOCATE (DefineVariable, char *, MAX (1, i));
+
+ /* look for -c FILENAME for config file */
+ if ((N = get_argument (*argc, argv, "-c"))) {
+ remove_argument (N, argc, argv);
+ filename = strcreate (argv[N]);
+ remove_argument (N, argc, argv);
+ return (filename);
+ }
+
+ /* look for PROGNAME env var */
+ find = strcreate (progname);
+ for (i = 0; i < strlen(find); i++) find[i] = toupper (find[i]);
+ filename = getenv (find);
+ free (find);
+ if (filename != (char *) NULL) {
+ find = strcreate (filename);
+ return (find);
+ }
+
+ uid = getuid();
+ gid = getgid();
+ ALLOCATE (find, char, strlen(progname) + 32);
+
+ /* look for progname.rc */
+ sprintf (find, "%s.rc", progname);
+ status = stat (find, &filestat);
+ if (status == 0) {
+ /* file exists, can we read it? */
+ if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR)) ||
+ ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP)) ||
+ ( (filestat.st_mode & S_IROTH))) {
+ return (find);
+ }
+ }
+
+ /* we eliminate this option: this is just confusing */
+# if (0)
+ /* look for ./.prognamerc */
+ sprintf (find, ".%src", progname);
+ status = stat (find, &filestat);
+ if (status == 0) {
+ /* file exists, can we read it? */
+ if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR)) ||
+ ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP)) ||
+ ( (filestat.st_mode & S_IROTH))) {
+ return (find);
+ }
+ }
+ free (find);
+# endif
+
+ /* look for ~/.prognamerc */
+ home = getenv ("HOME");
+ if (home == (char *) NULL) { return ((char *) NULL); }
+ ALLOCATE (find, char, 1024);
+ sprintf (find, "%s/.%src", home, progname);
+ status = stat (find, &filestat);
+ if (status == 0) {
+ /* file exists, can we read it? */
+ if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR)) ||
+ ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP)) ||
+ ( (filestat.st_mode & S_IROTH))) {
+ return (find);
+ }
+ }
+ free (find);
+ return ((char *) NULL);
+}
+
+char *LoadConfigFile (char *filename) {
+
+ FILE *f;
+ int i, done, Nbytes, NBYTES, nbytes, Nout, Ncpy, INPUT, Nlevel;
+ char *ibuffer, *obuffer, *tbuffer;
+ char *last, *next;
+ char infile[256], line[256];
+
+ /* open file */
+ f = fopen (filename, "r");
+ if (f == NULL) {
+ return ((char *) NULL);
+ }
+
+ /* allocate tmp space, 2 extra bytes for a final return and EOL */
+ Nbytes = 0;
+ NBYTES = D_NBYTES;
+ ALLOCATE (ibuffer, char, NBYTES + 2);
+
+ /* load data from file */
+ done = FALSE;
+ while ((nbytes = fread (&ibuffer[Nbytes], sizeof(char), D_NBYTES, f)) == D_NBYTES) {
+ Nbytes += nbytes;
+ NBYTES += D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES + 2);
+ }
+ Nbytes += nbytes;
+ fclose (f);
+
+ /* add final return & EOL if non-existent */
+ if (ibuffer[Nbytes-1] != '\n') {
+ ibuffer[Nbytes] = '\n';
+ Nbytes ++;
+ }
+ if (ibuffer[Nbytes]) ibuffer[Nbytes] = 0;
+
+ /* add the optional variables to the end of the input buffer */
+ if (ConfigVariable != (char *) NULL) {
+ snprintf (line, 256, "%s %s\n", "CONFIG", ConfigVariable);
+ Ncpy = strlen (line);
+ if (Nbytes + Ncpy >= NBYTES - 1) {
+ NBYTES = Nbytes + Ncpy + D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES);
+ }
+ memcpy (&ibuffer[Nbytes], line, Ncpy);
+ Nbytes += Ncpy;
+ ibuffer[Nbytes] = 0;
+ }
+ for (i = 0; i < NDefineVariable; i++) {
+ snprintf (line, 256, "%s %s\n", DefineVariable[i], DefineValue[i]);
+ Ncpy = strlen (line);
+ if (Nbytes + Ncpy >= NBYTES - 1) {
+ NBYTES = Nbytes + Ncpy + D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES);
+ }
+ memcpy (&ibuffer[Nbytes], line, Ncpy);
+ Nbytes += Ncpy;
+ ibuffer[Nbytes] = 0;
+ }
+
+ /* loop over input buffer, interpolating 'input' lines until none are added */
+ Nlevel = 0;
+ do {
+ INPUT = FALSE;
+
+ /* allocate output buffer & set counter */
+ NBYTES = Nbytes + D_NBYTES;
+ ALLOCATE (obuffer, char, NBYTES);
+ Nout = 0;
+
+ /* copy from ibuffer to obuffer, interpolating 'input' lines */
+ last = next = ibuffer;
+ for (i = 1; (next = ScanConfig (ibuffer, "input", "%s", i, infile)) != (char *) NULL; i++) {
+ /* copy data from last point to before 'input' */
+ Ncpy = next - last - 5;
+ if (Nout + Ncpy >= NBYTES - 1) {
+ NBYTES = Nout + Ncpy + D_NBYTES;
+ REALLOCATE (obuffer, char, NBYTES);
+ }
+ memcpy (&obuffer[Nout], last, Ncpy);
+ Nout += Ncpy;
+ obuffer[Nout] = 0;
+
+ /* insert data from 'input' file */
+ tbuffer = LoadRawConfigFile (infile, FALSE);
+ if (tbuffer != (char *) NULL) {
+ Ncpy = strlen (tbuffer);
+ if (Nout + Ncpy >= NBYTES - 1) {
+ NBYTES = Nout + Ncpy + D_NBYTES;
+ REALLOCATE (obuffer, char, NBYTES);
+ }
+ memcpy (&obuffer[Nout], tbuffer, Ncpy);
+ free (tbuffer);
+ Nout += Ncpy;
+ obuffer[Nout] = 0;
+ }
+
+ /* pointer goes to end of input line */
+ last = strchr (next, '\n');
+ if (last == (char *) NULL) break;
+ last ++;
+ INPUT = TRUE;
+ }
+ /* last set of bytes after last input */
+ Ncpy = strlen (last);
+ if (Nout + Ncpy >= NBYTES - 1) {
+ NBYTES = Nout + Ncpy + D_NBYTES;
+ REALLOCATE (obuffer, char, NBYTES);
+ }
+ memcpy (&obuffer[Nout], last, Ncpy);
+ Nout += Ncpy;
+ obuffer[Nout] = 0;
+ free (ibuffer);
+ ibuffer = obuffer;
+ Nlevel ++;
+ } while (INPUT && (Nlevel < 20));
+ if (Nlevel == 20) {
+ fprintf (stderr, "warning: config reached max depth of 20\n");
+ }
+
+ /* 'obuffer' now has complete set of interpolated lines from 'filename' */
+ return (obuffer);
+}
+
+char *ScanConfig (char *config, char *field, char *mode, int Nentry, ...) {
+
+ int i;
+ char *p, *p2, *tmp, *tfield, *start, *expandline();
+ va_list argp;
+ double value;
+
+ if (config == (char *) NULL) return ((char *) NULL);
+ va_start (argp, Nentry);
+
+ ALLOCATE (tfield, char, strlen (field) + 3);
+ sprintf (tfield, "\n%s", field);
+
+ /* we search for Nentry matching fields,
+ or until the end if Nentry == 0 */
+ p = (char *) NULL;
+ p2 = config;
+ for (i = 0; (i < Nentry) || !Nentry;) {
+ tmp = strstr (p2, tfield);
+ if (tmp == (char *) NULL) {
+ break;
+ }
+ p2 = tmp + strlen (tfield);
+ if (whitespace (*p2)) {
+ p = p2;
+ i++;
+ }
+ }
+ free (tfield);
+ if (Nentry && (i != Nentry)) {
+ p = va_arg (argp, char *);
+ p[0] = 0;
+ return ((char *) NULL);
+ }
+ if (p == (char *) NULL) {
+ p = va_arg (argp, char *);
+ p[0] = 0;
+ return ((char *) NULL);
+ }
+
+ start = p;
+ if (!strcmp (mode, "%s")) {
+ p2 = strchr (p, '\n');
+ if (p2 == (char *) NULL) p2 = config + strlen(config);
+ ALLOCATE (tmp, char, p2-p + 2);
+ memcpy (tmp, p, (p2-p));
+ tmp[(p2-p)] = 0;
+ stripwhite (tmp);
+ tmp = expandline (tmp, config);
+ p2 = va_arg (argp, char *);
+ strcpy (p2, tmp);
+ free (tmp);
+ } else {
+
+ /* try to get a numerical value from the field */
+ value = strtod (p, &p2);
+ if ((*p2 == 'd') || (*p2 == 'D'))
+ value *= pow (10.0, atof (p2 + 1));
+
+ if (!strcmp (mode, "%d")) *va_arg (argp, int *) = value;
+ if (!strcmp (mode, "%u")) *va_arg (argp, unsigned *) = value;
+ if (!strcmp (mode, "%ld")) *va_arg (argp, long *) = value;
+ if (!strcmp (mode, "%hd")) *va_arg (argp, short *) = value;
+ if (!strcmp (mode, "%f")) *va_arg (argp, float *) = value;
+ if (!strcmp (mode, "%lf")) *va_arg (argp, double *) = value;
+
+ }
+
+ va_end (argp);
+ return (start);
+
+}
+
+
+char *expandline (char *line, char *config) {
+
+ int Nin, Nout, Ncpy, NBYTES;
+ char *p1, *p2, word[256], value[256];
+ char *outline;
+
+ NBYTES = 256;
+ ALLOCATE (outline, char, NBYTES);
+ Nout = 0;
+ Nin = 0;
+ while ((p1 = strchr (&line[Nin], '$')) != (char *) NULL) {
+ Ncpy = p1 - line - Nin;
+ if (Nout + Ncpy >= NBYTES) {
+ NBYTES = Nout + Ncpy + 256;
+ REALLOCATE (outline, char, NBYTES - 2);
+ }
+ memcpy (&outline[Nout], &line[Nin], Ncpy);
+ Nout += Ncpy;
+
+ p1 ++;
+ for (p2 = p1; isalnum (*p2) || (*p2 == '_') || (*p2 == '-'); p2++);
+ memcpy (word, p1, p2 - p1);
+ word[p2-p1] = 0;
+ Nin += Ncpy + 1 + p2 - p1;
+
+ /* search for last entry of word */
+ if (!ScanConfig (config, word, "%s", 0, value)) {
+ fprintf (stderr, "variable %s not found in config file\n", word);
+ return (line);
+ }
+ Ncpy = strlen(value);
+ if (Nout + Ncpy >= NBYTES - 2) {
+ NBYTES = Nout + Ncpy + 256;
+ REALLOCATE (outline, char, NBYTES);
+ }
+ memcpy (&outline[Nout], value, Ncpy);
+ Nout += Ncpy;
+ }
+ Ncpy = strlen(&line[Nin]);
+ if (Nout + Ncpy >= NBYTES - 2) {
+ NBYTES = Nout + Ncpy + 256;
+ REALLOCATE (outline, char, NBYTES);
+ }
+ memcpy (&outline[Nout], &line[Nin], Ncpy);
+ Nout += Ncpy;
+ outline[Nout] = 0;
+ free (line);
+ return (outline);
+
+}
+
+char *LoadRawConfigFile (char *filename, int options) {
+
+ FILE *f;
+ int i, done, Nbytes, NBYTES, nbytes, Ncpy;
+ char *ibuffer;
+ char line[256];
+
+ /* open file */
+ f = fopen (filename, "r");
+ if (f == NULL) {
+ return ((char *) NULL);
+ }
+
+ /* allocate tmp space, 2 extra bytes for a final return and EOL */
+ Nbytes = 0;
+ NBYTES = D_NBYTES;
+ ALLOCATE (ibuffer, char, NBYTES + 2);
+
+ /* load data from file */
+ done = FALSE;
+ while ((nbytes = fread (&ibuffer[Nbytes], sizeof(char), D_NBYTES, f)) == D_NBYTES) {
+ Nbytes += nbytes;
+ NBYTES += D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES + 2);
+ }
+ Nbytes += nbytes;
+ fclose (f);
+
+ /* add final return & EOL, if non-existent */
+ if (ibuffer[Nbytes-1] != '\n') {
+ ibuffer[Nbytes] = '\n';
+ Nbytes ++;
+ }
+ if (ibuffer[Nbytes]) ibuffer[Nbytes] = 0;
+
+ if (options) {
+ /* write optional variables to bottom of buffer, overriding entries in the file */
+ if (ConfigVariable != (char *) NULL) {
+ snprintf (line, 256, "%s %s\n", "CONFIG", ConfigVariable);
+ Ncpy = strlen (line);
+ if (Nbytes + Ncpy >= NBYTES) {
+ NBYTES = Nbytes + Ncpy + D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES + 2);
+ }
+ memcpy (&ibuffer[Nbytes], line, Ncpy);
+ Nbytes += Ncpy;
+ ibuffer[Nbytes] = 0;
+ }
+ for (i = 0; i < NDefineVariable; i++) {
+ snprintf (line, 256, "%s %s\n", DefineVariable[i], DefineValue[i]);
+ Ncpy = strlen (line);
+ if (Nbytes + Ncpy >= NBYTES) {
+ NBYTES = Nbytes + Ncpy + D_NBYTES;
+ REALLOCATE (ibuffer, char, NBYTES + 2);
+ }
+ memcpy (&ibuffer[Nbytes], line, Ncpy);
+ Nbytes += Ncpy;
+ ibuffer[Nbytes] = 0;
+ }
+ }
+
+ return (ibuffer);
+}
+
+/*
+char *LoadConfigFile (char *filename) {
+
+ char *config;
+
+ config = LoadSubConfigFile (filename, TRUE);
+
+ return (config);
+
+}
+
+*/
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/coordops.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/coordops.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/coordops.c (revision 21775)
@@ -0,0 +1,536 @@
+# include
+# include
+
+int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) {
+
+ int Polynomial, Zenithal, Cartesian, PseudoCyl;
+ char *type;
+ double L, M, X, Y, T, Z;
+ 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 and PLY are identical & apply projection; WRP is cartesian */
+ Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
+ Zenithal = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+ Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+ PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+ 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;
+ 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 (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) {
+ 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 (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
+ 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")) {
+ Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
+ alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 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, 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");
+ Zenithal = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+ Cartesian = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+ PseudoCyl = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+ if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
+
+ /**** Locally Cartesian Projections ****/
+ if (Cartesian) {
+ 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 (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) {
+ L = DEG_RAD * sphi / stht;
+ M = -DEG_RAD * cphi / stht;
+ }
+ if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
+ L = DEG_RAD * sphi;
+ M = -DEG_RAD * cphi;
+ }
+ if (!strcmp(type, "-ZEA")) {
+ Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
+ L = Rc * sphi;
+ M = -Rc * cphi;
+ status = TRUE;
+ }
+ }
+
+ /**** Other Standard Projections ****/
+ if (PseudoCyl) {
+ if (!strcmp(type, "-AIT")) {
+ phi = RAD_DEG*(ra - coords[0].crval1);
+ theta = RAD_DEG*(dec - coords[0].crval2);
+ P = 1.0 + cos (theta) * cos (0.5*phi);
+ if (P != 0.0) {
+ A = DEG_RAD * sqrt (2.0 / P);
+ L = 2.0 * A * cos (theta) * sin (0.5*phi);
+ M = A * sin (theta);
+ } else {
+ L = 0.0;
+ M = 0.0;
+ }
+ }
+ if (!strcmp(type, "-GLS")) {
+ phi = ra - coords[0].crval1;
+ theta = dec - coords[0].crval2;
+ L = phi * cos(RAD_DEG * theta);
+ M = theta;
+ }
+ if (!strcmp(type, "-PAR")) {
+ phi = ra - coords[0].crval1;
+ theta = dec - coords[0].crval2;
+ L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
+ M = 180.0 * sin (RAD_DEG*theta/3.0);
+ }
+ }
+
+ /* convert L,M to X,Y */
+ determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+ X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M);
+ Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L);
+
+ /** extra polynomial terms **/
+ if (Polynomial) {
+ for (i = 0; i < 3; i++) {
+ Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
+ Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
+ if (coords[0].Npolyterms > 1) {
+ Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
+ Mo += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
+ }
+ if (coords[0].Npolyterms > 2) {
+ Lo += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
+ Mo += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
+ }
+ dL = (L - Lo);
+ dM = (M - Mo);
+
+ X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM);
+ Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dL);
+ }
+ }
+ /* check for correct size (iterate?) */
+
+ *x = X / coords[0].cdelt1 + coords[0].crpix1;
+ *y = Y / coords[0].cdelt2 + coords[0].crpix2;
+
+ return (status);
+
+}
+
+int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords) {
+
+ int status;
+ double tmpx, tmpy;
+
+ status = RD_to_XY (&tmpx, &tmpy, ra, dec, coords);
+ *x = tmpx;
+ *y = tmpy;
+
+ return (status);
+
+}
+
+int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords) {
+
+ int status;
+ double tmpr, tmpd;
+
+ status = XY_to_RD (&tmpr, &tmpd, x, y, coords);
+ *ra = tmpr;
+ *dec = tmpd;
+
+ return (status);
+
+}
+
+int GetCoords (Coords *coords, Header *header) {
+
+ int status, itmp, Polynomial;
+ double Lambda, rotate, scale;
+ double equinox;
+ char *ctype;
+
+ rotate = 0.0;
+ coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0;
+ coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0;
+ coords[0].pc1_1 = coords[0].pc2_2 = 1.0;
+ coords[0].pc2_1 = coords[0].pc1_2 = 0.0;
+ coords[0].Npolyterms = 1;
+ strcpy (coords[0].ctype, "NONE");
+
+ status = FALSE;
+ if (fits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) {
+ status = fits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1);
+ status &= fits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1);
+ status &= fits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2);
+ status &= fits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2);
+
+ if (fits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) {
+ status &= fits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2);
+ if (fits_scan (header, "CROTA2", "%f", 1, &rotate)) {
+ Lambda = coords[0].cdelt2 / coords[0].cdelt1;
+ coords[0].pc1_1 = cos(rotate*RAD_DEG);
+ coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
+ coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda;
+ coords[0].pc2_2 = cos(rotate*RAD_DEG);
+ }
+ if (fits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) {
+ status &= fits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2);
+ status &= fits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1);
+ status &= fits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2);
+ }
+
+ ctype = &coords[0].ctype[4];
+ Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
+
+ if (Polynomial) {
+ if (fits_scan (header, "NPLYTERM", "%d", 1, &itmp)) {
+ coords[0].Npolyterms = itmp;
+ } else {
+ coords[0].Npolyterms = 3;
+ }
+ switch (coords[0].Npolyterms) {
+ case 3:
+ status &= fits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]);
+ status &= fits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]);
+ status &= fits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]);
+ status &= fits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]);
+ status &= fits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]);
+ status &= fits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]);
+ status &= fits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]);
+ status &= fits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]);
+ case 2:
+ status &= fits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]);
+ status &= fits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]);
+ status &= fits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]);
+ status &= fits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]);
+ status &= fits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]);
+ status &= fits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]);
+ case 0:
+ case 1:
+ break;
+ }
+ }
+ } else {
+ if (fits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) {
+ status &= fits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2);
+ status &= fits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1);
+ status &= fits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2);
+ /* renormalize */
+ scale = hypot (coords[0].pc1_1, coords[0].pc1_2);
+ coords[0].cdelt1 = coords[0].cdelt2 = scale;
+ coords[0].pc1_1 /= scale;
+ coords[0].pc1_2 /= scale;
+ coords[0].pc2_1 /= scale;
+ coords[0].pc2_2 /= scale;
+ } else {
+ status = FALSE;
+ }
+ }
+ } else {
+ /* some of my thesis data uses this simple linear model - convert on read? */
+ if (fits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) {
+ status = fits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1);
+ status &= fits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2);
+ status &= fits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2);
+ status &= fits_scan (header, "DEC_X", "%f", 1, &coords[0].pc2_1);
+ status &= fits_scan (header, "DEC_Y", "%f", 1, &coords[0].pc2_2);
+ coords[0].crpix1 = coords[0].crpix2 = 0.0;
+ coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
+ strcpy (coords[0].ctype, "GENE");
+ }
+ }
+ if (status) {
+ if (!fits_scan (header, "EQUINOX", "%lf", 1, &equinox)) {
+ if (!fits_scan (header, "EPOCH", "%lf", 1, &equinox)) {
+ equinox = 2000.0;
+ }
+ }
+ if (fabs (equinox - 2000.0) > 0.1) {
+ coords_precess (&coords[0].crval1, &coords[0].crval2, equinox, 2000.0);
+ }
+ }
+ if (!status) {
+ fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype);
+ coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0;
+ coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0;
+ coords[0].pc1_1 = coords[0].pc2_2 = 1.0;
+ coords[0].pc2_1 = coords[0].pc1_2 = 0.0;
+ strcpy (coords[0].ctype, "NONE");
+ }
+ return (status);
+}
+
+int PutCoords (Coords *coords, Header *header) {
+
+ int OldAIPS;
+ char csys[16], ctype[16];
+ double rotate, Lambda;
+
+ /* modifications to the ctype? */
+ OldAIPS = FALSE;
+ fits_modify (header, "CTYPE2", "%s", 1, coords[0].ctype);
+ if (!strcmp(coords[0].ctype, "MM")) {
+ fits_modify (header, "CTYPE1", "%s", 1, "LL");
+ OldAIPS = TRUE;
+ } else {
+ strcpy (csys, "NONE");
+ if (!strncmp (coords[0].ctype, "DEC-", 4)) strcpy (csys, "RA--");
+ if (!strncmp (coords[0].ctype, "GLAT", 4)) strcpy (csys, "GLON");
+ if (!strncmp (coords[0].ctype, "ELAT", 4)) strcpy (csys, "ELON");
+ if (!strncmp (coords[0].ctype, "HLAT", 4)) strcpy (csys, "HLON");
+ if (!strncmp (coords[0].ctype, "SLAT", 4)) strcpy (csys, "SLON");
+ if (!strcmp (csys, "NONE")) return (FALSE);
+ sprintf (ctype, "%s-%s", csys, &coords[0].ctype[5]);
+ fits_modify (header, "CTYPE1", "%s", 1, ctype);
+ }
+
+ fits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1);
+ fits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2);
+ fits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1);
+ fits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2);
+ fits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1);
+ fits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2);
+
+ if (OldAIPS) {
+ Lambda = coords[0].cdelt2 / coords[0].cdelt1;
+ rotate = DEG_RAD*atan2 (coords[0].pc2_1*Lambda, coords[0].pc1_1);
+ fits_modify (header, "CROTA1", "%f", 1, rotate);
+ fits_modify (header, "CROTA2", "%f", 1, rotate);
+ return (TRUE);
+ }
+
+ fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
+ fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
+ fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
+ fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
+ fits_modify (header, "NPLYTERM", "%d", 1, coords[0].Npolyterms);
+
+ /* RA Terms */
+ if (coords[0].Npolyterms > 1) {
+ fits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]); /* polyterms[0]); */
+ fits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]); /* polyterms[1]); */
+ fits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]); /* polyterms[2]); */
+ }
+ if (coords[0].Npolyterms > 2) {
+ fits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]); /* polyterms[3]); */
+ fits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]); /* polyterms[4]); */
+ fits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]); /* polyterms[5]); */
+ fits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]); /* polyterms[6]); */
+ }
+
+ /* Dec Terms */
+ if (coords[0].Npolyterms > 1) {
+ fits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]); /* polyterms[7]); */
+ fits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]); /* polyterms[8]); */
+ fits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]); /* polyterms[9]); */
+ }
+ if (coords[0].Npolyterms > 2) {
+ fits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]); /* polyterms[10]); */
+ fits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]); /* polyterms[11]); */
+ fits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]); /* polyterms[12]); */
+ fits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]); /* polyterms[13]); */
+ }
+ return (TRUE);
+}
+
+void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch) {
+
+ double T;
+ double A, D, RA, DEC, zeta, z, theta;
+ double SA, CA, SD, CD;
+
+ T = (out_epoch - in_epoch) / 100.0;
+
+ zeta = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T);
+ theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T);
+ z = 0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T;
+
+ A = *ra;
+ D = *dec;
+ SD = cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D);
+ CD = sqrt (1 - SD*SD);
+ SA = sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD;
+ CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD;
+
+ DEC = DEG_RAD*asin(SD);
+ RA = DEG_RAD*atan2(SA, CA) + z;
+
+ if (RA < 0)
+ RA += 360;
+
+ *ra = RA;
+ *dec = DEC;
+}
+
+/* -PLY projection is an extrapolation of the -TAN projection.
+ In addition to the usual linear terms of CRPIXi, PC00i00j, there are
+ higher order polynomial terms (up to 3rd order):
+ Axis 1 terms:
+ PCA1X2Y0 = coords.polyterm[0][0] = x^2
+ PCA1X1Y1 = coords.polyterm[1][0] = xy
+ PCA1X0Y2 = coords.polyterm[2][0] = y^2
+ PCA1X3Y0 = coords.polyterm[3][0] = x^3
+ PCA1X2Y1 = coords.polyterm[4][0] = x^2 y
+ PCA1X1Y2 = coords.polyterm[5][0] = x y^2
+ PCA1X0Y3 = coords.polyterm[6][0] = y^2
+ Axis 2 terms:
+ PCA2X2Y0 = coords.polyterm[0][1] = x^2
+ PCA2X1Y1 = coords.polyterm[1][1] = xy
+ PCA2X0Y2 = coords.polyterm[2][1] = y^2
+ PCA2X3Y0 = coords.polyterm[3][1] = x^3
+ PCA2X2Y1 = coords.polyterm[4][1] = x^2 y
+ PCA2X1Y2 = coords.polyterm[5][1] = x y^2
+ PCA2X0Y3 = coords.polyterm[6][1] = y^2
+*/
+
+# if (0)
+
+ /** convert pixel coordinates to cartesian system **/
+ X = coords[0].cdelt1*(x - coords[0].crpix1);
+ Y = coords[0].cdelt2*(y - coords[0].crpix2);
+ if (Polynomi) {
+ if (coords[0].Npolyterms > 2) {
+ X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
+ Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]);
+ }
+ if (coords[0].Npolyterms > 2) {
+ X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]);
+ Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]);
+ }
+ }
+
+ L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
+ M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
+/** this code is the old method used for higher order terms. they
+ are essentially 6th order, with weird coupled terms.
+ I don't think any real data used these terms, but they should
+ be re-calculated, I would think
+**/
+
+# endif
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/findexec.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/findexec.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/findexec.c (revision 21775)
@@ -0,0 +1,251 @@
+# include
+
+int _check_permissions (char *filename) {
+
+ struct stat filestat;
+ uid_t uid;
+ gid_t gid;
+ int status;
+
+ uid = getuid();
+ gid = getgid();
+
+ /* check permission to exec file */
+ status = stat (filename, &filestat);
+ if (status == 0) { /* file exists, are permissions OK? */
+ if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IXUSR)) ||
+ ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IXGRP)) ||
+ ( (filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IXOTH))) {
+ return (TRUE);
+ } else {
+ return (FALSE);
+ }
+ }
+ return (FALSE);
+}
+
+/* pathname
+
+ given path/filename, returns path
+ given just filename, returns .
+ given path1/path2/, returns path1
+*/
+
+char *pathname (char *infile) {
+
+ int i;
+ char *c, *file;
+
+ /* make working version */
+ file = strcreate (infile);
+
+ /* strip off trailing / */
+ for (i = strlen (file); (i > 0) && (file[i] == '/'); i--) file[i] = 0;
+
+ c = strrchr (file, '/');
+ if (c == (char *) NULL) {
+ strcpy (file, ".");
+ } else {
+ *c = 0;
+ }
+
+ return (file);
+
+}
+
+/* filerootname
+
+ given /path/file.ext return file
+ given /path/file return file
+ given /path/file/ return file
+*/
+
+char *filerootname (char *infile) {
+
+ int i;
+ char *file, *root, *p1, *p2;
+
+ /* make working version */
+ file = strcreate (infile);
+
+ /* strip off trailing / */
+ for (i = strlen (file); (i > 0) && (file[i] == '/'); i--) file[i] = 0;
+
+ /* find last / */
+ p1 = strrchr (file, '/');
+ if (p1 == (char *) NULL)
+ p1 = file;
+ else
+ p1 ++;
+
+ /* find last . */
+ p2 = strrchr (file, '.');
+ if (p2 == (char *) NULL) p2 = p1 + strlen(p1);
+
+ /* create new string, free working space */
+ root = strncreate (p1, p2-p1);
+ free (file);
+
+ return (root);
+
+}
+
+/* fileextname
+
+ given /path/file.ext return ext
+ given /path/file return NULL
+ given /path/file/ return NULL
+ given file.ext return ext
+*/
+
+char *fileextname (char *file) {
+
+ char *root, *p1, *p2;
+
+ /* find last / */
+ p1 = strrchr (file, '/');
+ if (p1 == (char *) NULL) p1 = file;
+
+ /* find last . after p1 */
+ p2 = strrchr (p1, '.');
+ if (p2 == (char *) NULL) return ((char *) NULL);
+ p2 ++;
+
+ /* create new string, free working space */
+ root = strncreate (p2, strlen(p2));
+ return (root);
+
+}
+
+/* given /path/file.ext return file.ext */
+char *filebasename (char *name) {
+
+ char *c, *file;
+
+ ALLOCATE (file, char, strlen(name) + 1);
+ strcpy (file, name);
+ c = strrchr (file, '/');
+ if (c == (char *) NULL) {
+ strcpy (file, name);
+ } else {
+ strcpy (file, c+1);
+ }
+
+ return (file);
+
+}
+
+char *findexec (int argc, char **argv) {
+
+ int i, N, done, status;
+ char *c, *e, *dir, path[1024], name[1024];
+
+ /* if given an absolute or relative path, use it */
+ if (strchr (argv[0], '/') != (char *) NULL) {
+ status = _check_permissions (argv[0]);
+ if (status) {
+ if (realpath (argv[0], path) == (char *) NULL) return ((char *) NULL);
+ dir = pathname (path);
+ return (dir);
+ } else {
+ return ((char *) NULL);
+ }
+ }
+ N = 0;
+ for (i = argc+1; argv[i] != (char *) NULL; i++) {
+ if (!strncmp (argv[i], "PATH", 4)) {
+ N = i;
+ break;
+ }
+ }
+
+ if (N) {
+ c = &argv[N][5];
+ e = strchr (c, ':');
+ done = FALSE;
+ i = 0;
+ while (!done) {
+ bzero (path, 1024);
+ if (e == (char *) NULL) {
+ done = TRUE;
+ strncpy (path, c, strlen(c));
+ } else {
+ strncpy (path, c, e-c);
+ c = e+1;
+ e = strchr (c, ':');
+ }
+ sprintf (name, "%s/%s", path, argv[0]);
+ status = _check_permissions (name);
+
+ if (status) {
+ if (realpath (name, path) == (char *) NULL) continue;
+ dir = pathname (path);
+ return (dir);
+ }
+ }
+ }
+ return ((char *) NULL);
+}
+
+int mkdirhier (char *path) {
+
+ char *tpath;
+
+ if (mkdir (path, 0777)) {
+ if (errno == ENOENT) {
+ tpath = pathname (path);
+ if (!mkdirhier (tpath)) {
+ if (mkdir (path, 0777)) {
+ return (-1);
+ } else {
+ errno = 0;
+ return (0);
+ }
+ } else {
+ return (-1);
+ }
+ } else {
+ return (-1);
+ }
+ } else {
+ return (0);
+ }
+}
+
+char *getcwd_cfht (char *path, int size) {
+
+ char *hostname, *newpath, *p;
+
+ path = getcwd (path, size);
+
+ if (!strncmp (path, "/local/data", strlen ("/local/data"))) {
+
+ ALLOCATE (hostname, char, size);
+
+ if (gethostname (hostname, size-1)) {
+ fprintf (stderr, "ERROR: can't get hostname\n");
+ free (hostname);
+ return ((char *) NULL);
+ }
+ if ((p = strchr (hostname, '.')) != (char *) NULL) *p = 0;
+
+ ALLOCATE (newpath, char, size);
+ /* path might be just /local/data or /local/data/ */
+
+ p = path + strlen ("/local/data/");
+ if (strlen (path) <= strlen ("/local/data/")) {
+ snprintf (newpath, size, "/data/%s", hostname);
+ } else {
+ snprintf (newpath, size, "/data/%s/%s", hostname, p);
+ }
+
+ free (hostname);
+ strcpy (path, newpath);
+
+ free (newpath);
+
+ }
+
+ return (path);
+
+}
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/glockfile.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/glockfile.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/glockfile.c (revision 21775)
@@ -0,0 +1,298 @@
+# include
+
+/* lockfile concepts:
+
+ set a lockfile on [path]/filename.
+ we allow three types of locks:
+ LCK_HARD - exclusive lock with persistent lock file
+ LCK_XCLD - exclusive lock, no persistent lock file
+ LCK_SOFT - shared lock
+
+ both lock types set a filesystem-level lock on the given file,
+ using exclusive for hard and shared for soft locks
+
+ a hard lock also uses a lockfile to make the lock permanent:
+ the lockfile must be actively removed and will not vanish
+ even if the locking program dies
+
+ a soft lock will vanish if the locking program dies.
+
+ if a hard lock is set, both hard & soft locks will block
+ if a soft lock is set, a hard lock will block, but not a soft lock
+
+ the hard lockfile uses the file [path]/.filename.lck
+
+ setlockfile2 returns a file descriptor to an open READ/WRITE access file,
+ or -1 on failure, and sets state.
+*/
+
+FILE *fsetlockfile (char *filename, double timeout, int type, int *state) {
+
+ int fd;
+ FILE *f;
+
+ fd = setlockfile2 (filename, timeout, type, state);
+ if (fd == -1) return ((FILE *) NULL);
+
+ /* setlockfile2 returns a pointer to a vile open for update */
+ f = fdopen (fd, "r+");
+ if (f == (FILE *) NULL) {
+ clearlockfile2 (filename, fd, type, state);
+ fprintf (stderr, "can't make stream from descriptor\n");
+ return ((FILE *) NULL);
+ }
+ return (f);
+}
+
+/* clears lock. removes hardlock even if file pointer is not supplied */
+int fclearlockfile (char *filename, FILE *f, int type, int *state) {
+
+ int fd, status;
+
+ fd = -1;
+ if (f != NULL) {
+ fflush (f);
+ fd = fileno (f);
+ }
+ status = clearlockfile2 (filename, fd, type, state);
+ if (f != NULL) fclose (f);
+ return (status);
+}
+
+int setlockfile2 (char *filename, double timeout, int type, int *state) {
+
+ int i, done, nbytes, status;
+ char *lockname, buffer[64];
+ char *file, *path;
+ int fd, flock;
+ int filemode;
+ struct stat filestat;
+ struct flock filelock;
+ struct timeval now, then;
+
+ /* initial values at -1, close if the are set to another value */
+ fd = flock = -1;
+ filemode = S_IRUSR | S_IRGRP | S_IROTH;
+ file = path = lockname = (char *) NULL;
+ gettimeofday (&then, (void *) NULL);
+
+ /* define lock type */
+ filelock.l_start = 0;
+ filelock.l_whence = SEEK_SET;
+ filelock.l_len = 0;
+ filelock.l_pid = getpid ();
+ switch (type) {
+ case LCK_HARD:
+ case LCK_XCLD:
+ filelock.l_type = F_WRLCK; /* set an exclusive lock */
+ filemode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
+ break;
+ case LCK_SOFT:
+ filelock.l_type = F_RDLCK; /* set a shared lock */
+ filemode = S_IRUSR | S_IRGRP | S_IROTH;
+ break;
+ default:
+ fprintf (stderr, "invalid lock type %d\n", type);
+ }
+
+ /* if soft & file does not exist, return LCK_ACCESS */
+ if (type == LCK_SOFT) {
+ status = stat (filename, &filestat);
+ if ((status == -1) && (errno == ENOENT)) {
+ *state = LCK_EMPTY;
+ goto failure;
+ }
+ }
+
+ /* try to open file */
+ fd = open (filename, O_RDWR | O_CREAT, filemode);
+ if (fd == -1) {
+ *state = LCK_ACCESS;
+ goto failure;
+ }
+
+ /* we first try to set a FS level lock on the file */
+ for (done = 0; !done; ) {
+
+ /* check for timeout */
+ gettimeofday (&now, (void *) NULL);
+ if (DTIME (now, then) > timeout) {
+ *state = LCK_TIMEOUT;
+ goto failure;
+ }
+ usleep (1000);
+
+ /* try to lock file */
+ if (fcntl (fd, F_SETLK, &filelock) != -1) done = 1;
+ }
+
+ if (type == LCK_HARD) {
+ /* we've locked filename, now check lockname */
+
+ /* locking the lockfile before checking / setting contents will
+ ensure the data is synced across NFS */
+
+ /* set up name to lockfile */
+ path = pathname (filename);
+ file = filebasename (filename);
+ ALLOCATE (lockname, char, strlen (path) + strlen (file) + 10);
+ sprintf (lockname, "%s/.%s.lck", path, file);
+
+ /* try to open lockfile */
+ flock = open (lockname, O_RDWR | O_CREAT, filemode);
+ if (flock == -1) {
+ *state = LCK_HARDLCK;
+ goto failure;
+ }
+
+ /* try a few times to lock lockfile */
+ for (i = 0; (i < 20) && (fcntl (flock, F_SETLK, &filelock) == -1); i++) {
+ usleep (10000);
+ }
+ if (i == 20) {
+ *state = LCK_HARDLCK + 10;
+ goto failure;
+ }
+
+ /* we've locked the file, now read the contents */
+ nbytes = read (flock, buffer, 4);
+ if (nbytes == 4) { /* lock file has a word in it */
+ buffer[4] = 0;
+ if (!strcmp (buffer, "BUSY")) {
+ *state = LCK_HARDLCK + 20;
+ goto failure;
+ }
+ }
+
+ /* we've really got the lock, write BUSY to protect it */
+ lseek (flock, 0, SEEK_SET);
+ write (flock, "BUSY", 4);
+
+ /* now unlock lockfile */
+ filelock.l_type = F_UNLCK;
+ if (fcntl (flock, F_SETLK, &filelock) == -1) {
+ *state = LCK_HARDLCK + 30;
+ goto failure;
+ }
+ if (close (flock)) {
+ fprintf (stderr, "error closing lock\n");
+ goto failure;
+ }
+ flock = -1;
+ }
+
+ /* check if file is empty or not */
+ if (fstat (fd, &filestat)) {
+ *state = LCK_UNKNOWN;
+ } else {
+ if (filestat.st_size == 0) {
+ *state = LCK_EMPTY;
+ } else {
+ *state = LCK_FULL;
+ }
+ }
+
+ if (path != NULL) free (path);
+ if (file != NULL) free (file);
+ if (lockname != NULL) free (lockname);
+ return (fd);
+
+failure:
+ if (fd != -1) close (fd);
+ if (flock != -1) close (flock);
+ if (path != NULL) free (path);
+ if (file != NULL) free (file);
+ if (lockname != NULL) free (lockname);
+ return (-1);
+}
+
+/* if we pass in -1, all we can do is clear the hardlock */
+int clearlockfile2 (char *filename, int fd, int type, int *state) {
+
+ char *lockname;
+ int i, status;
+ char *path, *file;
+ int flock, filemode;
+ struct stat filestat;
+ struct flock filelock;
+
+ file = path = lockname = (char *) NULL;
+
+ /* define lock */
+ filelock.l_type = F_UNLCK;
+ filelock.l_start = 0;
+ filelock.l_whence = SEEK_SET;
+ filelock.l_len = 0;
+ filelock.l_pid = getpid ();
+ filemode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
+
+ /* first clear hard lockfile */
+ if (type == LCK_HARD) {
+
+ /* define lockfile */
+ path = pathname (filename);
+ file = filebasename (filename);
+ ALLOCATE (lockname, char, strlen (path) + strlen (file) + 10);
+ sprintf (lockname, "%s/.%s.lck", path, file);
+
+ /* check for lockfile existence */
+ status = stat (lockname, &filestat);
+ if (status == -1) {
+ *state = LCK_HUNLOCK;
+ goto failure;
+ }
+
+ /* open lockfile */
+ flock = open (lockname, O_RDWR | O_CREAT, filemode);
+ if (flock == -1) {
+ *state = LCK_HUNLOCK;
+ goto failure;
+ }
+
+ /* try a few times to lock lockfile */
+ filelock.l_type = F_WRLCK;
+ for (i = 0; (i < 20) && (fcntl (flock, F_SETLK, &filelock) == -1); i++) {
+ usleep (10000);
+ }
+ if (i == 20) {
+ *state = LCK_HUNLOCK;
+ goto failure;
+ }
+
+ /* set value to IDLE */
+ lseek (flock, 0, SEEK_SET);
+ write (flock, "IDLE", 4);
+
+ /* unlock lockfile */
+ filelock.l_type = F_UNLCK;
+ if (fcntl (flock, F_SETLK, &filelock) == -1) {
+ *state = LCK_HUNLOCK;
+ goto failure;
+ }
+ if (close (flock)) {
+ fprintf (stderr, "error closing lock\n");
+ goto failure;
+ }
+ }
+
+ /* now unlock the file */
+ if (fd != -1) {
+ if (fcntl (fd, F_SETLK, &filelock) == -1) {
+ *state = LCK_HUNLOCK;
+ goto failure;
+ }
+ }
+
+ *state = LCK_UNLOCK;
+ if (path != NULL) free (path);
+ if (file != NULL) free (file);
+ if (lockname != NULL) free (lockname);
+ return (1);
+
+failure:
+ if (path != NULL) free (path);
+ if (file != NULL) free (file);
+ if (lockname != NULL) free (lockname);
+ return (0);
+
+}
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/imreg_datatypes.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/imreg_datatypes.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/imreg_datatypes.c (revision 21775)
@@ -0,0 +1,51 @@
+# 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: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/phot_catalog.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/phot_catalog.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/phot_catalog.c (revision 21775)
@@ -0,0 +1,357 @@
+# include
+# include
+
+SecFilt *FixOldSecFilt (OldSecFilt *in, int Nvalues);
+Average *FixOldAverage (OldAverage *in, int Nvalues);
+Measure *FixOldMeasure (OldMeasure *in, int Nvalues);
+
+/* 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;
+
+ /* set lock on database, create stream f */
+ catalog[0].f = (FILE *) NULL;
+ catalog[0].f = fsetlockfile (catalog[0].filename, 3600.0, lockmode, &dbstate);
+ if (dbstate == LCK_EMPTY) return (2);
+ if (catalog[0].f == (FILE *) NULL) return (0);
+
+ fseek (catalog[0].f, 0, SEEK_SET);
+ catalog[0].lockmode = lockmode;
+
+ 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);
+}
+
+/* this is (a bit) dangerous: the load function will automatically
+ convert from the old (Loneos) format to the new (CFHT) format.
+ But, if you load an old-format catalog, then save the new-format
+ version, the new file will be corrupted unless you load all four
+ data components. Actually, 'save_catalog' assumes you've loaded
+ all four anyway, so it is not such a problem... */
+
+int load_catalog (Catalog *catalog, char mode, int VERBOSE) {
+
+ int Nitems, nitems;
+ int i, Nmeas, Nmiss, size, ssize;
+ int NewMeasure, NewAverage, NewSecFilt;
+ FILE *f;
+ OldMeasure *tmpmeasure;
+ OldAverage *tmpaverage;
+ OldSecFilt *tmpsecfilt;
+ struct stat filestatus;
+
+ f = catalog[0].f;
+ fseek (f, 0, SEEK_SET);
+
+ /* read catalog header */
+ if (!fits_load_header (f, &catalog[0].header)) {
+ if (VERBOSE) fprintf (stderr, "catalog file does not exist: %s\n", catalog[0].filename);
+ return (FALSE);
+ }
+
+ /* this is done by the load anyway... */
+ 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;
+ if (!fits_scan (&catalog[0].header, "NEWMEAS", "%t", 1, &NewMeasure)) NewMeasure = FALSE;
+ if (!fits_scan (&catalog[0].header, "NEWAVES", "%t", 1, &NewAverage)) NewAverage = FALSE;
+ if (!fits_scan (&catalog[0].header, "NEWSECF", "%t", 1, &NewSecFilt)) NewSecFilt = FALSE;
+
+ /* predicted file size - for double checking data validity */
+ size = catalog[0].header.size;
+ ssize = NewAverage ? sizeof (Average) : sizeof (OldAverage);
+ size += ssize * catalog[0].Naverage;
+ ssize = NewMeasure ? sizeof (Measure) : sizeof (OldMeasure);
+ size += ssize * catalog[0].Nmeasure;
+ ssize = sizeof (Missing);
+ size += ssize * catalog[0].Nmissing;
+ ssize = NewSecFilt ? sizeof (SecFilt) : sizeof (OldSecFilt);
+ size += ssize * 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");
+ ssize = NewAverage ? sizeof (Average) : sizeof (OldAverage);
+ fprintf (stderr, "average: %d = %d bytes\n", catalog[0].Naverage, catalog[0].Naverage*ssize);
+ ssize = NewMeasure ? sizeof (Measure) : sizeof (OldMeasure);
+ fprintf (stderr, "measure: %d = %d bytes\n", catalog[0].Nmeasure, catalog[0].Nmeasure*ssize);
+ ssize = sizeof (Missing);
+ fprintf (stderr, "missing: %d = %d bytes\n", catalog[0].Nmissing, catalog[0].Nmissing*ssize);
+ ssize = NewSecFilt ? sizeof (SecFilt) : sizeof (OldSecFilt);
+ fprintf (stderr, "secfilt: %d = %d bytes\n", catalog[0].Nsecfilt, catalog[0].Nsecfilt*catalog[0].Naverage*ssize);
+ 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);
+ }
+
+ if (mode & LOAD_AVES) {
+ /* read average values */
+ if (NewAverage) {
+ ALLOCATE (catalog[0].average, Average, MAX (catalog[0].Naverage, 1));
+ Nitems = catalog[0].Naverage;
+ nitems = Fread (catalog[0].average, sizeof(Average), Nitems, f, "average");
+ } else {
+ ALLOCATE (tmpaverage, OldAverage, MAX (catalog[0].Naverage, 1));
+ Nitems = catalog[0].Naverage;
+ nitems = Fread (tmpaverage, sizeof(OldAverage), Nitems, f, "oldaverage");
+ catalog[0].average = FixOldAverage (tmpaverage, Nitems);
+ }
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to read averages from catalog file %s (%d vs %d)\n", catalog[0].filename, nitems, Nitems);
+ return (FALSE);
+ }
+ } else {
+ /* skip over averages */
+ fseek (f, catalog[0].Naverage * sizeof(Average), SEEK_CUR);
+ }
+
+ if (mode & LOAD_MEAS) {
+ /* read measurements */
+ if (NewMeasure) {
+ ALLOCATE (catalog[0].measure, Measure, MAX (catalog[0].Nmeasure, 1));
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fread (catalog[0].measure, sizeof(Measure), Nitems, f, "measure");
+ } else {
+ ALLOCATE (tmpmeasure, OldMeasure, MAX (catalog[0].Nmeasure, 1));
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fread (tmpmeasure, sizeof(OldMeasure), Nitems, f, "oldmeasure");
+ catalog[0].measure = FixOldMeasure (tmpmeasure, Nitems);
+ }
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to read measures from catalog file %s (%d vs %d)\n", catalog[0].filename, nitems, Nitems);
+ return (FALSE);
+ }
+ } else {
+ /* skip over measures */
+ fseek (f, catalog[0].Nmeasure * sizeof(Measure), SEEK_CUR);
+ }
+
+ if (mode & LOAD_MISS) {
+ /* read missing */
+ ALLOCATE (catalog[0].missing, Missing, MAX (catalog[0].Nmissing, 1));
+ Nitems = catalog[0].Nmissing;
+ nitems = Fread (catalog[0].missing, sizeof(Missing), Nitems, f, "missing");
+ 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);
+ }
+ } else {
+ /* skip over missings */
+ fseek (f, catalog[0].Nmissing * sizeof(Missing), SEEK_CUR);
+ }
+
+ if (mode & LOAD_SECF) {
+ /* read missing */
+ if (NewSecFilt) {
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ ALLOCATE (catalog[0].secfilt, SecFilt, MAX (Nitems, 1));
+ nitems = Fread (catalog[0].secfilt, sizeof(SecFilt), Nitems, f, "secfilt");
+ } else {
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ ALLOCATE (tmpsecfilt, OldSecFilt, MAX (Nitems, 1));
+ nitems = Fread (tmpsecfilt, sizeof(OldSecFilt), Nitems, f, "oldsecfilt");
+ catalog[0].secfilt = FixOldSecFilt (tmpsecfilt, Nitems);
+ }
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to read secondary filters from catalog file %s (%d vs %d)\n", catalog[0].filename, nitems, Nitems);
+ return (FALSE);
+ }
+ } else {
+ /* skip over secfilts */
+ fseek (f, catalog[0].Naverage * catalog[0].Nsecfilt * sizeof(SecFilt), 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);
+}
+
+/* this is not robust: an error will not return us to the original
+ state. But, it is faster: only limited data written to disk? *
+*/
+
+int save_catalog (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);
+ fits_modify (&catalog[0].header, "NEWMEAS", "%t", 1, TRUE);
+ fits_modify (&catalog[0].header, "NEWAVES", "%t", 1, TRUE);
+ fits_modify (&catalog[0].header, "NEWSECF", "%t", 1, TRUE);
+
+ f = catalog[0].f;
+ /* rewind file pointers and truncate */
+ fseek (f, 0, SEEK_SET);
+ ftruncate (fileno (catalog[0].f), 0);
+
+ nitems = Fwrite (catalog[0].header.buffer, 1, catalog[0].header.size, f, "char");
+ if (nitems != catalog[0].header.size) {
+ if (VERBOSE) fprintf (stderr, "failed to write header\n");
+ return (FALSE);
+ }
+
+ Nitems = catalog[0].Naverage;
+ nitems = Fwrite (catalog[0].average, sizeof(Average), Nitems, f, "average");
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to write catalog file aves %s\n", catalog[0].filename);
+ return (FALSE);
+ }
+
+ Nitems = catalog[0].Nmeasure;
+ nitems = Fwrite (catalog[0].measure, sizeof(Measure), Nitems, f, "measure");
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to write catalog file meas %s\n", catalog[0].filename);
+ return (FALSE);
+ }
+
+ Nitems = catalog[0].Nmissing;
+ nitems = Fwrite (catalog[0].missing, sizeof(Missing), Nitems, f, "missing");
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to write catalog file miss %s\n", catalog[0].filename);
+ return (FALSE);
+ }
+
+ Nitems = catalog[0].Naverage * catalog[0].Nsecfilt;
+ nitems = Fwrite (catalog[0].secfilt, sizeof(SecFilt), Nitems, f, "secfilt");
+ if (nitems != Nitems) {
+ if (VERBOSE) fprintf (stderr, "failed to write catalog file secondary filters %s\n", catalog[0].filename);
+ return (FALSE);
+ }
+
+ return (TRUE);
+
+}
+
+/**** tools to convert between formats - this happens automatically in load_catalog */
+
+Measure *FixOldMeasure (OldMeasure *in, int Nvalues) {
+
+ int i;
+ Measure *out;
+
+ ALLOCATE (out, Measure, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].dR = in[i].dR;
+ out[i].dD = in[i].dD;
+ out[i].M = in[i].M;
+ out[i].dM = in[i].dM;
+ out[i].Mcal = in[i].Mcal;
+ out[i].dophot = in[i].dophot;
+ out[i].source = in[i].source;
+ out[i].t = in[i].t;
+
+ out[i].dt = 0xffff;
+ out[i].averef = in[i].average & 0x00ffffff;
+ out[i].flags = (in[i].average & 0xff000000) >> 24;
+ }
+
+ free (in);
+
+ return (out);
+}
+
+Average *FixOldAverage (OldAverage *in, int Nvalues) {
+
+ int i;
+ Average *out;
+
+ ALLOCATE (out, Average, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].R = in[i].R;
+ out[i].D = in[i].D;
+ out[i].M = in[i].M;
+ 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].dM = 0xffff;
+ }
+
+ free (in);
+
+ return (out);
+}
+
+SecFilt *FixOldSecFilt (OldSecFilt *in, int Nvalues) {
+
+ int i;
+ SecFilt *out;
+
+ ALLOCATE (out, SecFilt, Nvalues);
+
+ for (i = 0; i < Nvalues; i++) {
+ out[i].M = in[i].M;
+ out[i].Xm = in[i].Xm;
+
+ out[i].dM = 0xffff;
+ }
+
+ free (in);
+
+ return (out);
+}
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/photfits.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/photfits.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/photfits.c (revision 21775)
@@ -0,0 +1,254 @@
+# include
+# include
+# define MAX_ORDER 3
+
+/* 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 */
+/* all terms of order > 0 are stored as 16bit floats
+ the zero-order term is stored as a short int (-32k,+32k) */
+
+void returnMcal (Image *image, double *c) {
+
+ switch (image[0].order) {
+ case 0:
+ c[0] = image[0].Mcal;
+ return;
+ case 1:
+ c[0] = image[0].Mcal;
+ c[1] = getMi (image[0].Mx);
+ c[2] = getMi (image[0].My);
+ return;
+ case 2:
+ c[0] = image[0].Mcal;
+ 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;
+ 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;
+ 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 = c[0];
+ return;
+ case 1:
+ image[0].Mcal = c[0];
+ image[0].Mx = putMi(c[1]);
+ image[0].My = putMi(c[2]);
+ return;
+ case 2:
+ image[0].Mcal = 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 = 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 = 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 = 0.0;
+ image[0].order = 0;
+ return;
+ }
+}
+
+/* return value of image fit at x,y */
+double applyMcal (Image *image, double x, double y) {
+
+ double Mcal, c[15];
+
+ returnMcal (image, c); /* convert image parameters to c[i] coeffs */
+ Mcal = 0.0;
+ switch (image[0].order) {
+ case 0:
+ Mcal = c[0];
+ break;
+ case 1:
+ Mcal = c[0] + x*c[1] + y*c[2];
+ break;
+ case 2:
+ Mcal = c[0] + x*(c[1] + x*c[2]) + y*(c[3] + x*c[4] + y*c[5]);
+ break;
+ case 3:
+ Mcal = c[0] + x*(c[1] + x*(c[2] + x*c[3])) + y*(c[4] + x*(c[5] + x*c[6]) + y*(c[7] + x*c[8] + y*c[9]));
+ break;
+ case 4:
+ Mcal = c[0] + x*(c[1] + x*(c[2] + x*(c[3] + x*c[4]))) + y*(c[5] + x*(c[6] + x*(c[7] + x*c[8])) + y*(c[9] + x*(c[10] + x*c[11]) + y*(c[12] + x*c[13] + y*c[14])));
+ break;
+ }
+ return (Mcal);
+}
+
+double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order) {
+
+ int i;
+ double *x, *y, *m, *dm;
+ double S, S2, s, dS;
+
+ S2 = S = 0;
+ if (N < 2) {
+ return (0.0);
+ }
+
+ x = X; y = Y; m = M; dm = dM;
+ switch (order) {
+ case 0:
+ for (i = 0; i < N; i++, m++, dm++) {
+ dS = *m - c[0];
+ *dm = dS;
+ S += dS;
+ S2 += dS*dS;
+ }
+ break;
+ case 1:
+ for (i = 0; i < N; i++, m++, x++, y++, dm++) {
+ s = c[0] + *x*c[1] + *y*c[2];
+ dS = *m - s;
+ *dm = dS;
+ S += dS;
+ S2 += dS*dS;
+ }
+ break;
+ case 2:
+ for (i = 0; i < N; i++, m++, x++, y++, dm++) {
+ s = c[0] + *x*(c[1] + *x*c[2]) + *y*(c[3] + *x*c[4] + *y*c[5]);
+ dS = *m - s;
+ *dm = dS;
+ S += dS;
+ S2 += dS*dS;
+ }
+ break;
+ case 3:
+ for (i = 0; i < N; i++, m++, x++, y++, dm++) {
+ s = c[0] + *x*(c[1] + *x*(c[2] + *x*c[3])) + *y*(c[4] + *x*(c[5] + *x*c[6]) + *y*(c[7] + *x*c[8] + *y*c[9]));
+ dS = *m - s;
+ *dm = dS;
+ S += dS;
+ S2 += dS*dS;
+ }
+ break;
+ case 4:
+ for (i = 0; i < N; i++, m++, x++, y++, dm++) {
+ s = c[0] + *x*(c[1] + *x*(c[2] + *x*(c[3] + *x*c[4]))) + *y*(c[5] + *x*(c[6] + *x*(c[7] + *x*c[8])) + *y*(c[9] + *x*(c[10] + *x*c[11]) + *y*(c[12] + *x*c[13] + *y*c[14])));
+ dS = *m - s;
+ *dm = dS;
+ S += dS;
+ S2 += dS*dS;
+ }
+ break;
+ }
+ S = S / N;
+ S2 = sqrt (S2 / N - S*S);
+ return (S2);
+}
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/skydb.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/skydb.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/skydb.c (revision 21775)
@@ -0,0 +1,400 @@
+
+typedef struct {
+ float r, d;
+} SkyCoord;
+
+/* each region is bounded on the sky by lines of constant RA & DEC.
+ the region represents a specified depth, and it may or may not have
+ children. The start and end of the children in the array are childS and childE
+ If the region at the given depth is populated with an object table, then 'object' is TRUE.
+ If the region at the given depth is populated with an image table, then 'image' is TRUE.
+*/
+
+typedef struct {
+ char name[16];
+ float Rmin, Rmax;
+ float Dmin, Dmax;
+ int childS, childE;
+ char depth, child;
+ char object, image;
+} SkyRegion; /* 48 bytes */
+
+/* I have defined no accelerators other than the table hierarchy */
+
+/* find region which overlaps c at given depth (-1 : max depth) */
+SkyRegion *SkyFindPoint (SkyRegion *db, SkyCoord c, int depth) {
+
+ Ns = 0;
+ Ne = 1;
+
+ while (1) {
+ No = -1;
+ for (i = Ns; (No == -1) && (i < Ne); i++) {
+ if (c.r < db[i].Rmin) continue;
+ if (c.r > db[i].Rmax) continue;
+ if (c.d < db[i].Dmin) continue;
+ if (c.d > db[i].Dmax) continue;
+ // got the needed region at this depth
+ No = i;
+ }
+ if (No == -1) return ((SkyRegion *) NULL);
+ if (depth == db[No].depth) return (&db[No]);
+ if ((depth > db[No].depth) && !db[No].child) return ((SkyRegion *) NULL);
+
+ /* need to check Ns, Ne, or guarantee valid range */
+ Ns = db[No].childS;
+ Ne = db[No].childE;
+
+ if ((depth > db[No].depth) && db[No].child) continue;
+ return (&db[No]);
+ }
+}
+
+/* find regions at all levels which overlap c */
+SkyRegion **SkyFindLevels (SkyRegion *db, SkyCoord c, int *Nregion) {
+
+ SkyRegion **region;
+
+ Ns = 0;
+ Ne = 1;
+
+ N = 0;
+ NREGION = 10;
+ ALLOCATE (region, SkyRegion *, NREGION);
+
+ while (1) {
+ No = -1;
+ for (i = Ns; (No == -1) && (i < Ne); i++) {
+ if (c.r < db[i].Rmin) continue;
+ if (c.r > db[i].Rmax) continue;
+ if (c.d < db[i].Dmin) continue;
+ if (c.d > db[i].Dmax) continue;
+ // got the needed region at this depth
+ No = i;
+ }
+ if (No == -1) return ((SkyRegion *) NULL);
+
+ region[N] = &db[No];
+ N++;
+ if (N == NREGION) {
+ NREGION += 10;
+ REALLOCATE (region, SkyRegion *, NREGION);
+ }
+
+ /* need to check Ns, Ne, or guarantee valid range */
+ if (db[No].child) {
+ Ns = db[No].childS;
+ Ne = db[No].childE;
+ } else {
+ *Nregion = N;
+ return (region);
+ }
+ }
+}
+
+/* find regions contained within rectangular region c1 - c2 */
+SkyRegion **SkyFindArea (SkyRegion *db, SkyCoord c1, SkyCoord c2, *nlist) {
+
+ int Nlist;
+ SkyRegion *list, *new, *sub;
+
+ while (c1.r > 360.0) c1.r -= 360.0;
+ while (c1.r < 0.0) c1.r += 360.0;
+ while (c2.r > 360.0) c2.r -= 360.0;
+ while (c2.r < 0.0) c2.r += 360.0;
+
+ /* check on c1.r > c2.r : cross boundary */
+ if (c1.r > c2.r) {
+ c0 = c2; c0.r = 360.0;
+ list1 = SkyFindArea (db, c1, c0, &Nlist1);
+ c0 = c1; c0.r = 0.0;
+ list2 = SkyFindArea (db, c0, c2, &Nlist2);
+ Nlist = Nlist1 + Nlist2;
+ ALLOCATE (list, SkyRegion *, Nlist);
+ memcpy (&list[0], list1, Nlist1*sizeof());
+ memcpy (&list[Nlist1], list2, Nlist2*sizeof());
+ free (list1);
+ free (list2);
+ *nlist = Nlist;
+ return (list);
+ }
+
+ Nlist = 1;
+ ALLOCATE (list, SkyRegion *, 1);
+ list[0] = &db[0];
+ getchild = db[0].child;
+
+ while (getchild) {
+
+ getchild = FALSE;
+ Nnew = 0;
+ NNEW = Nlist + 100;
+ ALLOCATE (new, SkyRegion *, NNEW);
+
+ for (i = 0; i < Nlist; i++) {
+ if (list[i][0].child) {
+ sub = SkyFindAreaDB(db, list[i], c1, c2, &Nsub);
+ if (Nnew + Nsub == NNEW) {
+ NNEW += 100;
+ REALLOCATE (new, SkyRegion *, NNEW);
+ }
+ for (i = 0; i < Nsub; i++) {
+ getchild |= sub[i][0].child;
+ new[Nnew] = sub[i];
+ Nnew++;
+ }
+ free (sub);
+ } else {
+ new[Nnew] = list[i];
+ Nnew ++;
+ if (Nnew == NNEW) {
+ NNEW += 100;
+ REALLOCATE (new, SkyRegion *, NNEW);
+ }
+ }
+ }
+
+ free (list);
+ list = new;
+ Nlist = Nnew;
+ }
+
+ return (list);
+ *nlist = Nlist;
+
+}
+
+/* ?? */
+SkyRegion *SkyFindAreaDB (SkyRegion *db, SkyRegion *ref, SkyCoord c1, SkyCoord c2, int *Nregion) {
+
+ SkyRegion **region;
+
+ Ns = ref[0].childS;
+ Ne = ref[0].childE;
+
+ N = 0;
+ NREGION = 100;
+ ALLOCATE (region, SkyRegion *, NREGION);
+
+ /* c1 min, c2 max */
+
+ for (i = Ns; (No == -1) && (i < Ne); i++) {
+ if (c2.r < db[i].Rmin) continue;
+ if (c1.r > db[i].Rmax) continue;
+ if (c2.d < db[i].Dmin) continue;
+ if (c1.d > db[i].Dmax) continue;
+
+ region[N] = &db[i];
+ N++;
+ if (N == NREGION) {
+ NREGION += 100;
+ REALLOCATE (region, SkyRegion *, NREGION);
+ }
+ }
+ *Nregion = N;
+ return (region);
+}
+
+SkyRegion *SkyBuildTable (SkyRegion *seed, int Nseed, int level, int depth) {
+
+ /* given a table of Rmin, Dmin, Rmax, Dmax, name, at level 3, generate
+ all higher levels and 'depth' extra levels */
+
+ /* levels:
+ 0 - fullsky.cpt
+ 1 - n????.cpt / s????.cpt
+ 2 - r????.cpt
+ 3 - ????.cpt
+ 4 - ????.??.cpt
+ */
+
+ SkyRegion *db;
+
+ /* allocate at least 30 for levels 0 & 1 */
+ NREGION = 100;
+ ALLOCATE (db, SkyRegion, NREGION);
+
+ /* full sky */
+ strcpy (db[0].name, "fullsky.cpt");
+ db[0].Rmin = 0; db[0].Rmax = 360;
+ db[0].Dmin = -90; db[0].Dmax = 90;
+ db[0].depth = 0;
+ db[0].child = FALSE;
+ N = 1;
+
+ db[0].childS = N;
+ /* north dec bands */
+ for (dec = 0; dec < 90; dec += 7.5) {
+ sprintf (db[N].name, "n%04d.cpt", (int) 100*dec);
+ db[N].Rmin = 0; db[N].Rmax = 360;
+ db[N].Dmin = dec; db[N].Dmax = dec + 7.5;
+ db[N].depth = 1;
+ db[N].child = FALSE;
+ N++;
+ }
+ /* south dec bands */
+ for (dec = 0; dec > -90; dec -= 7.5) {
+ sprintf (db[N].name, "s%04d.cpt", (int) 100*dec);
+ db[N].Rmin = 0; db[N].Rmax = 360;
+ db[N].Dmin = dec - 7.5; db[N].Dmax = dec;
+ db[N].depth = 1;
+ db[N].child = FALSE;
+ N++;
+ }
+ db[0].childE = N;
+
+ /* subdivide dec bands based on seed */
+ Rnumber = 0;
+ childS = db[0].childS;
+ childE = db[0].childE;
+ for (i = childS; i < childE; i++) {
+
+ Nnew = 0;
+ NNEW = 100;
+ ALLOCATE (new, SkyRegion, NNEW);
+
+ for (j = 0; j < Nseed; j++) {
+ if (seed[j].Rmin > db[i].Rmax) continue;
+ if (seed[j].Rmax < db[i].Rmin) continue;
+ if (seed[j].Dmin > db[i].Dmax) continue;
+ if (seed[j].Dmax < db[i].Dmin) continue;
+
+ for (k = 0; k < Nnew; k++) {
+ if ((seed[j].Rmin == new[k].Rmin) &&
+ (seed[j].Rmax == new[k].Rmax)) {
+ goto next_seed;
+ }
+ if ((seed[j].Rmin == new[k].Rmin) ^^
+ (seed[j].Rmax != new[k].Rmax)) {
+ fprintf (stderr, "inconsistent blocks in seed file\n");
+ return (NULL);
+ }
+ }
+ new[Nnew].Rmin = seed[j].Rmin;
+ new[Nnew].Rmax = seed[j].Rmax;
+ new[Nnew].Dmin = db[i].Dmin;
+ new[Nnew].Dmax = db[i].Dmax;
+ new[Nnew].depth = 2;
+ new[Nnew].child = FALSE;
+ strncpy (root, db[i].name, 5); root[5] = 0;
+ sprintf (new[Nnew].name, "%s/r%04d.cpt", root, Rnumber);
+ Rnumber ++;
+ Nnew ++;
+ if (Nnew == NNEW) {
+ NNEW += 100;
+ REALLOCATE (new, SkyRegion, NNEW);
+ }
+ next_seed:
+ }
+ /* update db list */
+ db[i].childS = N;
+ db[i].childE = N + Nnew;
+ NREGION = N + Nnew + 100;
+ REALLOCATE (db, SkyRegion, NREGION);
+ memcpy (&db[N], new, Nnew*sizeof(SkyRegion));
+ free (new);
+ N += Nnew;
+ }
+
+ /* subdivide ra strips based on seed */
+ childS = db[db[0].childS].childS;
+ childE = db[db[0].childE - 1].childE;
+ nextS = N;
+ for (i = childS; i < childE; i++) {
+
+ Nnew = 0;
+ NNEW = 100;
+ ALLOCATE (new, SkyRegion, NNEW);
+
+ for (j = 0; j < Nseed; j++) {
+ if (seed[j].Rmin > db[i].Rmax) continue;
+ if (seed[j].Rmax < db[i].Rmin) continue;
+ if (seed[j].Dmin > db[i].Dmax) continue;
+ if (seed[j].Dmax < db[i].Dmin) continue;
+
+ for (k = 0; k < Nnew; k++) {
+ if ((seed[j].Dmin == new[k].Dmin) &&
+ (seed[j].Dmax == new[k].Dmax)) {
+ goto next_seed;
+ }
+ if ((seed[j].Dmin == new[k].Dmin) ^^
+ (seed[j].Dmax != new[k].Dmax)) {
+ fprintf (stderr, "inconsistent blocks in seed file\n");
+ return (NULL);
+ }
+ }
+ new[Nnew].Dmin = seed[j].Dmin;
+ new[Nnew].Dmax = seed[j].Dmax;
+ new[Nnew].Rmin = db[i].Rmin;
+ new[Nnew].Rmax = db[i].Rmax;
+ new[Nnew].depth = 3;
+ new[Nnew].child = FALSE;
+ strcpy (new[Nnew].name, seed[j].name);
+ Nnew ++;
+ if (Nnew == NNEW) {
+ NNEW += 100;
+ REALLOCATE (new, SkyRegion, NNEW);
+ }
+ next_seed:
+ }
+
+ /* update db list */
+ db[i].childS = N;
+ db[i].childE = N + Nnew;
+ NREGION = N + Nnew + 100;
+ REALLOCATE (db, SkyRegion, NREGION);
+ memcpy (&db[N], new, Nnew*sizeof(SkyRegion));
+ free (new);
+ N += Nnew;
+ }
+ nextE = N;
+
+ /* subdivide entries once more */
+ childS = nextS;
+ childE = nextE;
+ for (i = childS; i < childE; i++) {
+ Rnumber = 0;
+ db[i].childS = N;
+ dDec = (db[i].Dmax - db[i].Dmin) / 5.0;
+ dRa = (db[i].Rmax - db[i].Rmin) / 5.0;
+ for (nx = 0; nx < 5; nx++) {
+ for (ny = 0; ny < 5; ny++) {
+ db[N].Dmin = db[i].Dmin + dDec * (nx + 0);
+ db[N].Dmax = db[i].Dmax + dDec * (nx + 1);
+ db[N].Rmin = db[i].Rmin + dDec * (nx + 0);
+ db[N].Rmax = db[i].Rmax + dDec * (nx + 1);
+ db[N].depth = 3;
+ db[N].child = FALSE;
+ strncpy (root, db[i].name, 10); root[10] = 0;
+ sprintf (db[N].name, "%s.%02d.cpt", root, Rnumber);
+ Rnumber ++;
+ N ++;
+ if (N == NREGION) {
+ NREGION += 100;
+ REALLOCATE (new, SkyRegion, NREGION);
+ }
+ }
+ }
+ }
+
+}
+
+/* region is pointer to entry in db */
+SkyRegion *SkyExtend (SkyRegion *db, SkyRegion *region) {
+
+}
+
+/* region is pointer to entry in db */
+SkyRegion *SkyContract (SkyRegion *db, SkyRegion *region) {
+
+}
+
+SkyRegion *SkyFindGCircle (SkyRegion *db, SkyCoord c1, SkyCoord c2) {
+
+}
+
+
+/* is entire db in memory, or do we load from disk for each operation?
+ 10,000 at level 3, 100,000 @ 4 : 0.5MB, 5MB */
+
+
Index: /tags/ohana/libohana-1-2/Ohana/src/libohana/src/string.c
===================================================================
--- /tags/ohana/libohana-1-2/Ohana/src/libohana/src/string.c (revision 21775)
+++ /tags/ohana/libohana-1-2/Ohana/src/libohana/src/string.c (revision 21775)
@@ -0,0 +1,157 @@
+# include
+
+#ifndef whitespace
+#define whitespace(c) (((c) == ' ') || ((c) == '\t'))
+#endif
+
+/* Strip whitespace from the start and end of STRING. */
+int stripwhite (char *string) {
+
+ int i;
+
+ if (string == (char *) NULL) return (FALSE);
+
+ for (i = 0; whitespace (string[i]); i++);
+ if (i) memmove (string, string + i, strlen(string+i)+1);
+ for (i = strlen (string) - 1; (i > 0) && whitespace (string[i]); i--);
+ string[++i] = 0;
+ return (i);
+
+}
+
+/* compare two strings either as strings, or as numbers if both are
+ pure numeric strings (base 10) */
+int strnumcmp (char *str1, char *str2) {
+
+ char *end1;
+ char *end2;
+ int num1, num2;
+ int value;
+
+ value = FALSE;
+ num1 = strtol (str1, &end1, 10);
+ num2 = strtol (str2, &end2, 10);
+
+ if (!end1[0] && !end2[0] && (num1 == num2)) {
+ value = TRUE;
+ }
+ if (!strcmp (str1, str2)) {
+ value = TRUE;
+ }
+
+ return (value);
+}
+
+/* create a new string from this string */
+char *strcreate (char *string) {
+
+ char *line;
+
+ if (string == (char *) NULL) return ((char *) NULL);
+
+ ALLOCATE (line, char, MAX (1, strlen(string)) + 1);
+ line = strcpy (line, string);
+
+ return (line);
+}
+
+/* create a new string of length n from this string */
+char *strncreate (char *string, int n) {
+
+ char *line;
+
+ if (string == (char *) NULL) return ((char *) NULL);
+
+ ALLOCATE (line, char, n + 1);
+ memcpy (line, string, n);
+ line[n] = 0;
+ return (line);
+}
+
+int scan_line (FILE *f, char *line) {
+
+ int i, status;
+ char c;
+
+ status = EOF + 1;
+
+ for (i = 0, c = 0; (c != '\n') && (status != EOF); i++) {
+ status = fscanf (f, "%c", &c);
+ line[i] = c;
+ }
+ line[i - 1] = 0; /* this could make things crash! */
+
+ if (i > 1) {
+ status = EOF + 1;
+ }
+
+ return (status);
+}
+
+char *_parse_nextword (char *string) {
+
+ if (string == (char *) NULL) return ((char *) NULL);
+
+ for (; isspace (*string); string++);
+ for (; (*string != 0) && !isspace (*string); string++);
+ for (; isspace (*string); string++);
+ return (string);
+}
+
+int dparse (double *X, int NX, char *line) {
+
+ int i;
+ char *word;
+ char *ptr;
+
+ word = line;
+ for (i = 0; i < NX - 1; i++)
+ word = _parse_nextword (word);
+
+ *X = strtod (word, &ptr);
+ if (ptr == word) return (FALSE);
+ if (word[0] == '-') return (-1);
+ return (1);
+}
+
+int fparse (float *X, int NX, char *line) {
+
+ int i;
+ char *word;
+ char *ptr;
+
+ word = line;
+ for (i = 0; i < NX - 1; i++)
+ word = _parse_nextword (word);
+
+ *X = strtod (word, &ptr);
+ if (ptr == word) return (FALSE);
+ if (word[0] == '-') return (-1);
+ return (1);
+}
+
+int get_argument (int argc, char **argv, char *arg) {
+
+ int i;
+
+ for (i = 0; i < argc; i++) {
+ if (!strcmp(argv[i], arg))
+ return (i);
+ }
+
+ return ((int)NULL);
+}
+
+int remove_argument (int N, int *argc, char **argv) {
+
+ int i;
+
+ if ((N != (int)NULL) && (N != 0)) {
+ (*argc)--;
+ for (i = N; i < *argc; i++) {
+ argv[i] = argv[i+1];
+ }
+ }
+
+ return (N);
+}