Index: trunk/Ohana/src/libdvo/Makefile =================================================================== --- trunk/Ohana/src/libdvo/Makefile (revision 36528) +++ trunk/Ohana/src/libdvo/Makefile (revision 36630) @@ -103,4 +103,5 @@ $(SRC)/ImageMetadata.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/RegionHostTable.$(ARCH).o \ $(SRC)/match_image.$(ARCH).o \ $(SRC)/db_utils.$(ARCH).o \ Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 36528) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 36630) @@ -285,4 +285,46 @@ } HostTable; +// A RegionHost processes data for some region in parallel with other regions +typedef struct RegionHostInfo { + double Rmin; // (Rmin,Rmax),(Dmin,Dmax) arehard RA,DEC boundaries of the + double Rmax; // region for which each host is responsible. A given host + double Dmin; // calibrates the images for which the fiducial point (center) + double Dmax; // lands in the region, and all objects in the region + + double RminCat; // (RminCat,RmaxCat),(DminCat,DmaxCat) are the region for which + double RmaxCat; // the catalogs need to be loaded : this is the outer bounds + double DminCat; // of the region containing all images completely + double DmaxCat; + + char *hostname; + + int hostID; // remove machine ID in SkyTable + int stdio[3]; // fd's for communication with the remote host + int pid; // remote process ID + int status; + IOBuffer stdout; + IOBuffer stderr; + + off_t Nimage; + off_t NIMAGE; + Image *image; + off_t *imseq; + + int *neighbors; // list of neighbor index values + int Nneighbors; // number of neighbors + char isNeighbor; // TRUE if I am a neighbor to the current region host +} RegionHostInfo; + +typedef struct { + double Rmin; + double Rmax; + double Dmin; + double Dmax; + + int Nhosts; + RegionHostInfo *hosts; + short *index; +} RegionHostTable; + // special-case function: CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, off_t *Ndata, char *swapped); @@ -385,4 +427,6 @@ uint32_t flags; int catID; + int objID; + int nOwn; } AverageTiny; @@ -410,4 +454,5 @@ short dYccd; short dRsys; + char myDet; } MeasureTiny; @@ -555,4 +600,5 @@ float *X; float *Y; + int *nOwn; // relastro uses this to count owned detections per object } Catalog; @@ -616,6 +662,4 @@ float PhotInst (Measure *measure); float PhotCat (Measure *measure); -float PhotAper (Measure *measure); -float PhotKron (Measure *measure); float PhotSys (Measure *measure, Average *average, SecFilt *secfilt); float PhotRel (Measure *measure, Average *average, SecFilt *secfilt); @@ -623,12 +667,24 @@ float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt); float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotAveErr (PhotCode *code, Average *average, SecFilt *secfilt); + +float PhotAperInst (Measure *measure); +float PhotAperCat (Measure *measure); +float PhotAperSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotAperRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotAperCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotAperRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); + +float PhotKronInst (Measure *measure); +float PhotKronCat (Measure *measure); +float PhotKronSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotKronRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotKronCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotKronRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt); + float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt); - -float PhotAperInst (Measure *measure); -float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotKronInst (Measure *measure); -float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt); float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt); @@ -855,3 +911,11 @@ void dvo_measureT_init (MeasureTiny *measure); +void InitRegionHosts (RegionHostInfo *hosts, int Nhosts, int NHOSTS); +void FreeRegionHosts (RegionHostInfo *hosts, int Nhosts); +void FreeRegionHostTable (RegionHostTable *table); +RegionHostTable *RegionHostTableLoad (char *catdir, char *rootname); +int RegionHostTableWaitJobs (RegionHostTable *regionHosts, char *file, int lineno); +int RegionHostTableWaitJobsGetIO (RegionHostTable *regionHosts, char *file, int lineno, int VERBOSE); +int RegionHostFindNeighbors (RegionHostTable *table, int Nhost); + # endif // DVO_H Index: trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvodb.h (revision 36528) +++ trunk/Ohana/src/libdvo/include/dvodb.h (revision 36630) @@ -23,31 +23,60 @@ /* magnitude types */ enum {MAG_NONE, - MAG_INST, - MAG_CAT, - MAG_APER, - MAG_APER_INST, - MAG_KRON, - MAG_KRON_INST, - MAG_KRON_ERR, - MAG_SYS, - MAG_REL, - MAG_CAL, + + // these magnitude types are from the measurement table (or derived therein) + // the following imply PSF magnitudes: + MAG_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_CAT, // MAG_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_SYS, // MAG_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_REL, // MAG_SYS - Mcal [specific zero point for image] + MAG_CAL, // MAG_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_ERR, // error on PSF MAG + + // the following imply APERTURE magnitudes: + MAG_APER_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_APER_CAT, // MAG_APER_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_APER_SYS, // MAG_APER_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_APER_REL, // MAG_APER_SYS - Mcal [specific zero point for image] + MAG_APER_CAL, // MAG_APER_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_APER_ERR, // error on APER MAG (if known, else psf error) + MAG_APER, // in mextract -> MAG_APER_REL, in avextract -> MAG_APER_AVE + + // the following imply KRON magnitudes: + MAG_KRON_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_KRON_CAT, // MAG_KRON_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_KRON_SYS, // MAG_KRON_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_KRON_REL, // MAG_KRON_SYS - Mcal [specific zero point for image] + MAG_KRON_CAL, // MAG_KRON_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_KRON_ERR, // error on KRON MAG + MAG_KRON, // in mextract -> MAG_KRON_REL, in avextract -> MAG_KRON_AVE + + // these magnitude types are from the average (secfilt) table: + // the following imply PSF magnitudes (or fluxes) MAG_AVE, MAG_REF, - MAG_ERR, MAG_STDEV, MAG_AVE_ERR, + MAG_CHISQ, + MAG_20, + MAG_80, + MAG_FLUX_PSF, + MAG_FLUX_PSF_ERR, + + // the following imply APERTURE magnitudes: + MAG_APER_AVE, + MAG_APER_REF, + + // the following imply KRON magnitudes or fluxes: + MAG_KRON_AVE, + MAG_KRON_REF, + MAG_FLUX_KRON, + MAG_FLUX_KRON_ERR, + + // other secfilt fields of interest MAG_PHOT_FLAGS, - MAG_CHISQ, MAG_NCODE, MAG_NPHOT, - MAG_20, - MAG_80, MAG_UC_DIST, MAG_STACK_DET_ID, - MAG_FLUX_PSF, - MAG_FLUX_PSF_ERR, - MAG_FLUX_KRON, - MAG_FLUX_KRON_ERR, }; Index: trunk/Ohana/src/libdvo/src/RegionHostTable.c =================================================================== --- trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 36630) +++ trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 36630) @@ -0,0 +1,497 @@ +# include "dvo.h" +# include +# include + +void InitRegionHosts (RegionHostInfo *hosts, int Nhosts, int NHOSTS) { + + int i; + for (i = Nhosts; i < NHOSTS; i++) { + hosts[i].Rmin = NAN; + hosts[i].Rmax = NAN; + hosts[i].Dmin = NAN; + hosts[i].Dmax = NAN; + + hosts[i].RminCat = NAN; + hosts[i].RmaxCat = NAN; + hosts[i].DminCat = NAN; + hosts[i].DmaxCat = NAN; + + hosts[i].hostname = NULL; + hosts[i].stdio[HOST_STDIN] = -1; + hosts[i].stdio[HOST_STDOUT] = -1; + hosts[i].stdio[HOST_STDERR] = -1; + hosts[i].pid = 0; + + hosts[i].image = NULL; + hosts[i].imseq = NULL; + hosts[i].Nimage = 0; + hosts[i].NIMAGE = 0; + + hosts[i].neighbors = NULL; + hosts[i].Nneighbors = 0; + hosts[i].isNeighbor = FALSE; + } + return; +} + +void FreeRegionHosts (RegionHostInfo *hosts, int Nhosts) { + + int i; + for (i = 0; i < Nhosts; i++) { + free (hosts[i].hostname); + } + free (hosts); + return; +} + +void FreeRegionHostTable (RegionHostTable *table) { + + if (!table) return; + if (table->hosts) { + FreeRegionHosts (table->hosts, table->Nhosts); + } + if (table->index) free (table->index); + free (table); + return; +} + +RegionHostTable *RegionHostTableLoad (char *catdir, char *rootname) { + + int i, Nline; + + char *filename = NULL; + + ALLOCATE (filename, char, strlen(catdir) + strlen(rootname) + 2); // one slash and one EOL + sprintf (filename, "%s/%s", catdir, rootname); + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "failed to open host table %s\n", filename); + free (filename); + return NULL; + } + + // simple format: ID hostname pathname + + int NHOSTS = 16; + int Nhosts = 0; + RegionHostInfo *hosts = NULL; + ALLOCATE (hosts, RegionHostInfo, NHOSTS); + InitRegionHosts (hosts, Nhosts, NHOSTS); + + int maxID = 0; + + RegionHostTable *table = NULL; + ALLOCATE (table, RegionHostTable, 1); + table->Rmin = 360; + table->Rmax = 0; + table->Dmin = +90; + table->Dmax = -90; + + for (Nline = 0; TRUE; Nline ++) { + int ID; + char tmphost[1024]; + char line[1024]; + + // XXXX use this for safety: int status = scan_line_maxlen (f, line, 1024); + int status = scan_line (f, line); + if (status == EOF) break; + + // find first non-whitespace char & skip commented lines + for (i = 0; OHANA_WHITESPACE (line[i]); i++); + if (line[i] == '#') continue; + if (line[i] == 0) continue; + + double Rmin, Rmax, Dmin, Dmax; + status = sscanf (line, "%d %1023s %lf %lf %lf %lf", &ID, tmphost, &Rmin, &Rmax, &Dmin, &Dmax); + if (status != 6) { + fprintf (stderr, "error reading line %d of region host table %s\n", Nline, filename); + FreeRegionHosts (hosts, Nhosts); + free (filename); + fclose (f); + return NULL; + } + + // check the validity of ID (0 < ID < MAX_SHORT) + + if (ID < 1) { + fprintf (stderr, "invalid host ID %d\n", ID); + exit (1); + } + if (ID > 255) { + fprintf (stderr, "invalid host ID %d\n", ID); + exit (1); + } + maxID = MAX(maxID, ID); + + hosts[Nhosts].hostID = ID; + hosts[Nhosts].hostname = strcreate(tmphost); + + InitIOBuffer (&hosts[Nhosts].stdout, 1000); + InitIOBuffer (&hosts[Nhosts].stderr, 1000); + + hosts[Nhosts].Rmin = Rmin; + hosts[Nhosts].Rmax = Rmax; + hosts[Nhosts].Dmin = Dmin; + hosts[Nhosts].Dmax = Dmax; + + hosts[Nhosts].RminCat = Rmin; + hosts[Nhosts].RmaxCat = Rmax; + hosts[Nhosts].DminCat = Dmin; + hosts[Nhosts].DmaxCat = Dmax; + + table->Rmin = MIN(Rmin, table->Rmin); + table->Rmax = MAX(Rmax, table->Rmax); + table->Dmin = MIN(Dmin, table->Dmin); + table->Dmax = MAX(Dmax, table->Dmax); + + // InitIOBuffer (&hosts[Nhosts].stdout, 1000); + // InitIOBuffer (&hosts[Nhosts].stderr, 1000); + + Nhosts ++; + if (Nhosts >= NHOSTS) { + NHOSTS += 16; + REALLOCATE (hosts, RegionHostInfo, NHOSTS); + InitRegionHosts (hosts, Nhosts, NHOSTS); + } + } + + table->Nhosts = Nhosts; + table->hosts = hosts; + + ALLOCATE (table->index, short, maxID + 1); + for (i = 0; i <= maxID; i++) table->index[i] = -1; + + for (i = 0; i < table->Nhosts; i++) { + if (table->index[table->hosts[i].hostID] != -1) { + fprintf (stderr, "error: duplicate hostID %d\n", table->hosts[i].hostID); + exit (1); + } + table->index[table->hosts[i].hostID] = i; + } + + free (filename); + fclose (f); + return table; +} + +int RegionHostFindNeighbors (RegionHostTable *table, int Nhost) { + + int i; + + // given a specific host (by table sequence), find all of its neighbors + // (eg, (host->Rmin == althost->Rmax && (host->Dmin <= althost->Dmax) && (host->Dmax >= althost->Dmin) + + RegionHostInfo *myhost = &table->hosts[Nhost]; + myAssert (!myhost->neighbors, "myhost neighbors already allocated"); + myAssert (!myhost->Nneighbors, "myhost Nneighbors not zero?"); + ALLOCATE (myhost->neighbors, int, 1); // always allocate 1 extra + + for (i = 0; i < table->Nhosts; i++) { + + if (i == Nhost) continue; + RegionHostInfo *altHost = &table->hosts[i]; + + int onBorder; + + onBorder = + (myhost->Rmin == altHost->Rmax) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Rmax == altHost->Rmin) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + // handle the 0,360 boundary (Rmin,Rmax are in the range 0,360, but 0 == 360) + onBorder = + (myhost->Rmin == 0.0) && + (altHost->Rmax == 360.0) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Rmax == 360.0) && + (altHost->Rmin == 0.0) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Dmin == altHost->Dmax) && + (myhost->Rmin <= altHost->Rmax) && + (myhost->Rmax >= altHost->Rmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Dmax == altHost->Dmin) && + (myhost->Rmin <= altHost->Rmax) && + (myhost->Rmax >= altHost->Rmin); + if (onBorder) goto add_neighbor; + + continue; + + add_neighbor: + myhost->neighbors[myhost->Nneighbors] = i; + myhost->Nneighbors ++; + REALLOCATE (myhost->neighbors, int, myhost->Nneighbors + 1); + altHost->isNeighbor = TRUE; + } + return TRUE; +} + +// wait for all children to complete, report output to stdout +int RegionHostTableWaitJobs (RegionHostTable *table, char *file, int lineno) { + + int i; + + // we have launched table->Nhosts jobs; wait for all of them to complete... + // if one (N) failed to launch, we will get an ECHILD error from the last (N) calls + int done = FALSE; + for (i = 0; !done && (i < table->Nhosts); i++) { + int status = 0; + + // XXX we'll need to pass in WNOHANG and keep retrying if we want to have a timeout... + int pid = waitpid (-1, &status, 0); + if (!pid) { + // this should only occur if we called waitpid with the WNOHANG option + fprintf (stderr, "programming error (1)? %s %d", file, lineno); + exit (2); + } + if (pid == -1) { + switch (errno) { + case ECHILD: + done = TRUE; + break; + default: + fprintf (stderr, "programming error (2)? %s %d", file, lineno); + exit (2); + } + } + + // when the host has finished, close the open sockets + + // find the host which has finished + int Nout, j; + int found = FALSE; + for (j = 0; j < table->Nhosts; j++) { + if (table->hosts[j].pid != pid) continue; + found = TRUE; + // check on the status of this and report any output? + fprintf (stderr, "job finished for %s (%d)\n", table->hosts[j].hostname, pid); + // read the stderr and stdout + IOBuffer buffer; + InitIOBuffer (&buffer, 100); + EmptyIOBuffer (&buffer, 100, table->hosts[j].stdio[HOST_STDOUT]); + fprintf (stderr, "--- stdout from %s ---\n", table->hosts[j].hostname); + Nout = write (STDOUT_FILENO, buffer.buffer, buffer.Nbuffer); + if (Nout != buffer.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + fprintf (stderr, "\n"); + + InitIOBuffer (&buffer, 100); + EmptyIOBuffer (&buffer, 100, table->hosts[j].stdio[HOST_STDERR]); + fprintf (stderr, "--- stderr from %s ---\n", table->hosts[j].hostname); + Nout = write (STDOUT_FILENO, buffer.buffer, buffer.Nbuffer); + if (Nout != buffer.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + fprintf (stderr, "\n"); + if (WIFEXITED(status)) { + fprintf (stderr, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + table->hosts[j].status = WEXITSTATUS(status); + if (table->hosts[j].status) { + fprintf (stderr, "job failed on %s\n", table->hosts[j].hostname); + continue; + } + } else { + table->hosts[j].status = -1; + fprintf (stderr, "job exited abnormally on %s\n", table->hosts[j].hostname); + continue; + } + } + if (!found) { + fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); + exit (2); + } + } + return TRUE; +} + + +// wait for all children to complete, report output to stdout +int RegionHostTableWaitJobsGetIO (RegionHostTable *table, char *file, int lineno, int VERBOSE) { + + // we have launched table->Nhosts jobs; wait for all of them to complete... + // if one (N) failed to launch, we will get an ECHILD error from the last (N) calls + + // we need to read any data waiting on stderr or stdout from these jobs, or the overfull + // buffers can cause a problem. we alternate between 'select' and 'waitpid' calls with + // timeouts for both + + // add all hosts' sockets to the fd_sets + fd_set rdSet, wtSet; + FD_ZERO (&rdSet); + FD_ZERO (&wtSet); + + // XXX can I set the fd_sets once, since I am not actually closing the fd's? + + int globalStatus = TRUE; + + int i; + int Nmax = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + FD_SET (table->hosts[i].stdio[HOST_STDIN], &wtSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDIN]); + FD_SET (table->hosts[i].stdio[HOST_STDOUT], &rdSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDOUT]); + FD_SET (table->hosts[i].stdio[HOST_STDERR], &rdSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDERR]); + } + Nmax ++; + + // need the list of connected hosts for exit test below + int Nrunning = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + Nrunning ++; + } + + int Nfound = 0; + + // this loop has 2 chunks: (a) check for I/O + (b) check for jobs done + while (1) { + + // Wait up to 0.5 second for host to provide I/O + // timeout gets mucked: need to reset before each select + struct timeval timeout; + timeout.tv_sec = 10; + timeout.tv_usec = 500000; + + int status = select (Nmax, NULL, &wtSet, NULL, &timeout); + if (status == -1) { + perror("select()"); + exit (2); + } + + // we have some sockets to check, check sockets for all hosts + for (i = 0; (status > 0) && (i < table->Nhosts); i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + + if (FALSE && FD_ISSET (table->hosts[i].stdio[HOST_STDIN], &wtSet)) { + // this host is waiting for input : this is an error, so exit + fprintf (stderr, "host %s is waiting for input\n", table->hosts[i].hostname); + abort(); + } + + if ((table->hosts[i].stdio[HOST_STDOUT] > 0) && FD_ISSET (table->hosts[i].stdio[HOST_STDOUT], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i].stdout, table->hosts[i].stdio[HOST_STDOUT]); + // if (table->hosts[i].stdout.Nbuffer > 0x10000) { + if (table->hosts[i].stdout.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i].stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i].stdout.buffer, table->hosts[i].stdout.Nbuffer); + if (Nout != table->hosts[i].stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i].stdout); + if (printHead) fprintf (stdout, "\n"); + } + } + + if ((table->hosts[i].stdio[HOST_STDERR] > 0) && FD_ISSET (table->hosts[i].stdio[HOST_STDERR], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i].stderr, table->hosts[i].stdio[HOST_STDERR]); + // if (table->hosts[i].stderr.Nbuffer > 0x10000) { + if (table->hosts[i].stderr.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i].stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stderr.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i].stderr.buffer, table->hosts[i].stderr.Nbuffer); + if (Nout != table->hosts[i].stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i].stderr); + if (printHead) fprintf (stdout, "\n"); + } + } + } + + // now check if any children have finished... + while (TRUE) { + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + if (!pid) { + // fprintf (stderr, "no hosts to harvest\n"); + usleep (500000); + break; // no outstanding jobs have finished + } + if ((pid == -1) && (errno == ECHILD)) goto escape; // no more jobs on which to wait + if ((pid == -1) && (errno != ECHILD)) { + fprintf (stderr, "programming error (2)? %s %d", file, lineno); + exit (2); + } + + // find the host which has finished + int found = FALSE; + for (i = 0; (i < table->Nhosts) && !found; i++) { + if (table->hosts[i].pid != pid) continue; + found = TRUE; + + RegionHostInfo *host = &table->hosts[i]; + + // check on the status of this and report any output? + if (VERBOSE) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + // read stdout + int printHead; + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (VERBOSE) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + globalStatus = FALSE; + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + globalStatus = FALSE; + continue; + } + } + if (!found) { + fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); + exit (2); + } + Nfound ++; + if (Nfound == Nrunning) goto escape; // we've harvested all jobs + } + } + +escape: + + // close all opened connections + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + close (table->hosts[i].stdio[HOST_STDIN]); + close (table->hosts[i].stdio[HOST_STDOUT]); + close (table->hosts[i].stdio[HOST_STDERR]); + } + + return globalStatus; +} Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36630) @@ -211,5 +211,5 @@ } } else { - value.Flt = PhotdM (field->photcode, average, secfilt); + value.Flt = PhotAveErr (field->photcode, average, secfilt); } break; @@ -244,9 +244,19 @@ case MAG_APER: + case MAG_APER_AVE: value.Flt = PhotAperAve (field->photcode, average, secfilt); break; + case MAG_APER_REF: + value.Flt = PhotAperRef (field->photcode, average, secfilt, measure); + break; + case MAG_KRON: + case MAG_KRON_AVE: value.Flt = PhotKronAve (field->photcode, average, secfilt); + break; + + case MAG_KRON_REF: + value.Flt = PhotKronRef (field->photcode, average, secfilt, measure); break; @@ -283,5 +293,5 @@ break; case AVE_dMAG: - value.Flt = PhotdM (field->photcode, average, secfilt); + value.Flt = PhotAveErr (field->photcode, average, secfilt); break; case AVE_Xm: Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 36630) @@ -123,10 +123,20 @@ case MAG_REL: case MAG_CAL: - case MAG_APER: + case MAG_ERR: + case MAG_APER_INST: - case MAG_KRON: + case MAG_APER_CAT: + case MAG_APER_SYS: + case MAG_APER_REL: + case MAG_APER_CAL: + case MAG_APER_ERR: + case MAG_KRON_INST: + case MAG_KRON_CAT: + case MAG_KRON_SYS: + case MAG_KRON_REL: + case MAG_KRON_CAL: case MAG_KRON_ERR: - case MAG_ERR: + case MAG_PHOT_FLAGS: equiv = myEquiv; @@ -137,12 +147,23 @@ case MAG_AVE: case MAG_REF: + case MAG_STDEV: + case MAG_AVE_ERR: case MAG_CHISQ: - case MAG_AVE_ERR: + case MAG_FLUX_PSF: + case MAG_FLUX_PSF_ERR: + + case MAG_APER_AVE: + + case MAG_KRON_AVE: + case MAG_FLUX_KRON: + case MAG_FLUX_KRON_ERR: + + case MAG_20: + case MAG_80: + case MAG_UC_DIST: + case MAG_STACK_DET_ID: + case MAG_NCODE: case MAG_NPHOT: - case MAG_FLUX_PSF: - case MAG_FLUX_PSF_ERR: - case MAG_FLUX_KRON: - case MAG_FLUX_KRON_ERR: equiv = field->photcode; goto valid_photcode; @@ -178,40 +199,71 @@ value.Flt = PhotRef (equiv, average, secfilt, measure); break; - case MAG_APER: - value.Flt = PhotAper (measure); - break; + case MAG_ERR: + value.Flt = measure[0].dM; + break; + case MAG_AVE_ERR: + value.Flt = PhotAveErr (equiv, average, secfilt); + break; + case MAG_APER_INST: value.Flt = PhotAperInst (measure); break; - case MAG_KRON: - value.Flt = PhotKron (measure); - break; + case MAG_APER_CAT: + value.Flt = PhotAperCat (measure); + break; + case MAG_APER_SYS: + value.Flt = PhotAperSys (measure, average, secfilt); + break; + case MAG_APER: + case MAG_APER_REL: + value.Flt = PhotAperRel (measure, average, secfilt); + break; + case MAG_APER_CAL: + value.Flt = PhotAperCal (measure, average, secfilt, measure, equiv); + break; + case MAG_APER_AVE: + value.Flt = PhotAperAve (equiv, average, secfilt); + break; + case MAG_APER_REF: + value.Flt = PhotAperRef (equiv, average, secfilt, measure); + break; + case MAG_APER_ERR: + value.Flt = measure[0].dM; + break; + case MAG_KRON_INST: value.Flt = PhotKronInst (measure); break; + case MAG_KRON_CAT: + value.Flt = PhotKronCat (measure); + break; + case MAG_KRON_SYS: + value.Flt = PhotKronSys (measure, average, secfilt); + break; + case MAG_KRON: + case MAG_KRON_REL: + value.Flt = PhotKronRel (measure, average, secfilt); + break; + case MAG_KRON_CAL: + value.Flt = PhotKronCal (measure, average, secfilt, measure, equiv); + break; + case MAG_KRON_AVE: + value.Flt = PhotKronAve (equiv, average, secfilt); + break; + case MAG_KRON_REF: + value.Flt = PhotKronRef (equiv, average, secfilt, measure); + break; case MAG_KRON_ERR: - value.Flt = measure[0].dMkron; - break; - case MAG_ERR: - value.Flt = measure[0].dM; - break; - case MAG_AVE_ERR: - value.Flt = PhotdM (equiv, average, secfilt); - break; + value.Flt = measure[0].dMkron; + break; + case MAG_PHOT_FLAGS: if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) { value.Int = measure[0].photFlags; - } else { - value.Int = 0; } break; - case MAG_CHISQ: - value.Flt = PhotXm (equiv, average, secfilt); - break; + case MAG_NCODE: - if (equiv == NULL) { - value.Int = 0; - break; - } + if (equiv == NULL) break; Nsec = GetPhotcodeNsec (equiv->code); if (Nsec == -1) break; @@ -219,12 +271,26 @@ break; case MAG_NPHOT: - if (equiv == NULL) { - value.Int = 0; - break; - } + if (equiv == NULL) break; Nsec = GetPhotcodeNsec (equiv->code); if (Nsec == -1) break; value.Int = secfilt[Nsec].Nused; break; + + case MAG_CHISQ: + value.Flt = PhotXm (equiv, average, secfilt); + break; + case MAG_STDEV: + value.Flt = PhotMstdev (equiv, average, secfilt); + break; + case MAG_20: + value.Flt = PhotM20 (equiv, average, secfilt); + break; + case MAG_80: + value.Flt = PhotM80 (equiv, average, secfilt); + break; + case MAG_UC_DIST: + value.Flt = PhotUCdist (equiv, average, secfilt); + break; + case MAG_FLUX_PSF: value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt); @@ -531,7 +597,8 @@ fieldc = MatchFieldMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); // fieldc = MatchField (measure[0].t, measure[0].photcode); + fieldc = MatchFieldMetadata (measure[0].imageID); } if (fieldc == NULL) break; @@ -547,7 +614,8 @@ fieldc = MatchFieldMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); // fieldc = MatchField (measure[0].t, measure[0].photcode); + fieldc = MatchFieldMetadata (measure[0].imageID); } if (fieldc == NULL) break; @@ -564,7 +632,8 @@ mosaic = MatchMosaicMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); + // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaicMetadata (measure[0].imageID); } if (mosaic == NULL) break; @@ -580,7 +649,8 @@ mosaic = MatchMosaicMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); + // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaicMetadata (measure[0].imageID); } if (mosaic == NULL) break; Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 36630) @@ -30,4 +30,5 @@ int GetMagMode (char *string) { + // these all imply PSF mags: if (!strcasecmp (string, "inst")) return (MAG_INST); if (!strcasecmp (string, "cat")) return (MAG_CAT); @@ -37,14 +38,45 @@ if (!strcasecmp (string, "ave")) return (MAG_AVE); if (!strcasecmp (string, "ref")) return (MAG_REF); + if (!strcasecmp (string, "err")) return (MAG_ERR); + if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); + + // these are explicit PSF mags: + if (!strcasecmp (string, "psf_inst")) return (MAG_INST); + if (!strcasecmp (string, "psf_cat")) return (MAG_CAT); + if (!strcasecmp (string, "psf_sys")) return (MAG_SYS); + if (!strcasecmp (string, "psf_rel")) return (MAG_REL); + if (!strcasecmp (string, "psf_cal")) return (MAG_CAL); + if (!strcasecmp (string, "psf_ave")) return (MAG_AVE); + if (!strcasecmp (string, "psf_ref")) return (MAG_REF); + if (!strcasecmp (string, "psf_err")) return (MAG_ERR); + if (!strcasecmp (string, "psf_aveerr")) return (MAG_AVE_ERR); + + // these are aper mags: if (!strcasecmp (string, "ap")) return (MAG_APER); if (!strcasecmp (string, "aper")) return (MAG_APER); if (!strcasecmp (string, "aperinst")) return (MAG_APER_INST); if (!strcasecmp (string, "aper_inst")) return (MAG_APER_INST); + if (!strcasecmp (string, "aper_cat")) return (MAG_APER_CAT); + if (!strcasecmp (string, "aper_sys")) return (MAG_APER_SYS); + if (!strcasecmp (string, "aper_rel")) return (MAG_APER_REL); + if (!strcasecmp (string, "aper_cal")) return (MAG_APER_CAL); + if (!strcasecmp (string, "aper_ave")) return (MAG_APER_AVE); + if (!strcasecmp (string, "aper_ref")) return (MAG_APER_REF); + if (!strcasecmp (string, "aper_err")) return (MAG_APER_ERR); + + // these are kron mags: if (!strcasecmp (string, "kron")) return (MAG_KRON); + if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); + if (!strcasecmp (string, "kron_cat")) return (MAG_KRON_CAT); + if (!strcasecmp (string, "kron_sys")) return (MAG_KRON_SYS); + if (!strcasecmp (string, "kron_rel")) return (MAG_KRON_REL); + if (!strcasecmp (string, "kron_cal")) return (MAG_KRON_CAL); + if (!strcasecmp (string, "kron_ave")) return (MAG_KRON_AVE); + if (!strcasecmp (string, "kron_ref")) return (MAG_KRON_REF); + if (!strcasecmp (string, "kron_err")) return (MAG_KRON_ERR); if (!strcasecmp (string, "kroninst")) return (MAG_KRON_INST); - if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); if (!strcasecmp (string, "kronerr")) return (MAG_KRON_ERR); - if (!strcasecmp (string, "err")) return (MAG_ERR); - if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); + + // other fields: if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); if (!strcasecmp (string, "flags")) return (MAG_PHOT_FLAGS); Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36630) @@ -119,4 +119,5 @@ average->measureOffset = -1; average->catID = 0; + average->nOwn = 0; } @@ -242,4 +243,5 @@ measure->dYccd = 0; measure->dRsys = 0; + measure->myDet = FALSE; } Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36630) @@ -549,4 +549,5 @@ } +// XXX note that there are 2 bad versions of PS1_V4 : the other one does have have PSF_QF_PERFECT, but has PAD instead int gfits_convert_Measure_PS1_V4alt (Measure_PS1_V4alt *data, off_t size, off_t nitems) { @@ -708,5 +709,5 @@ # ifdef BYTE_SWAP byte = (unsigned char *) data; - for (i = 0; i < nitems; i++, byte += 129) { + for (i = 0; i < nitems; i++, byte += 120) { /** BYTE SWAP **/ SWAP_DBLE (0); // RA @@ -754,5 +755,5 @@ Ncols = ftable[0].header[0].Naxis[0]; if (Ncols != 120) { - fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, 176); + fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, 120); return NULL; } Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 36630) @@ -169,4 +169,5 @@ } +// returns Nsec if code is PRI/SEC, else -1 int GetPhotcodeNsec (int code) { @@ -259,76 +260,4 @@ } -float PhotAper (Measure *measure) { - - int Np; - float Mcat; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Mcat = measure[0].Map; - return (Mcat); - } - code = &photcodes[0].code[Np]; - Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (Mcat); -} - -float PhotAperInst (Measure *measure) { - - int Np; - float Minst; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Minst = measure[0].Map; - return (Minst); - } - Minst = measure[0].Map - measure[0].dt - ZERO_POINT; - - return (Minst); -} - -float PhotKron (Measure *measure) { - - int Np; - float Mcat; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Mcat = measure[0].Mkron; - return (Mcat); - } - code = &photcodes[0].code[Np]; - Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (Mcat); -} - -float PhotKronInst (Measure *measure) { - - int Np; - float Minst; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Minst = measure[0].Mkron; - return (Minst); - } - Minst = measure[0].Mkron - measure[0].dt - ZERO_POINT; - - return (Minst); -} - float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) { @@ -362,23 +291,4 @@ } -float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt) { - - int Np; - float ZP; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - ZP = 0.0; - return (ZP); - } - code = &photcodes[0].code[Np]; - ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (ZP); -} - float PhotRel (Measure *measure, Average *average, SecFilt *secfilt) { @@ -444,4 +354,435 @@ Mcal += Mcol; return (Mcal); +} + +/***/ +float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mave; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + return (Mave); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +float PhotAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dM; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dM = (Ns == -1) ? NAN : secfilt[Ns].dM; + return (dM); +} + +/************************************** APERTURE-related Magnitudes ***********************************************/ + +float PhotAperInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Map; + return (Minst); + } + Minst = measure[0].Map - measure[0].dt - ZERO_POINT; + + return (Minst); +} + +float PhotAperCat (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotAperSys (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mcat); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotAperRel (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mrel); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotAperCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + if (code == NULL) return NAN; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].Map; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) return (NAN); + + Mcal = PhotAperRel (thisone, average, secfilt) + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mcal); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Maper; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Maper = (Ns == -1) ? NAN : secfilt[Ns].Map; + return (Maper); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotAperRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].Map; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +/************************************** KRON-related Magnitudes ***********************************************/ + +float PhotKronInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Mkron; + return (Minst); + } + Minst = measure[0].Mkron - measure[0].dt - ZERO_POINT; + + return (Minst); +} + +float PhotKronCat (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Mkron; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotKronSys (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Msys = measure[0].Mkron; + return (Msys); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mcat); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotKronRel (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Mkron; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mrel); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotKronCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + if (code == NULL) return NAN; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].Mkron; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) return (NAN); + + Mcal = PhotKronRel (thisone, average, secfilt) + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mcal); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mkron = (Ns == -1) ? NAN : secfilt[Ns].Mkron; + return (Mkron); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotKronRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].Mkron; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dMkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dMkron = (Ns == -1) ? NAN : secfilt[Ns].dMkron; + return (dMkron); +} + +/********************* other support functions ********************************/ + +float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt) { + + int Np; + float ZP; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + ZP = 0.0; + return (ZP); + } + code = &photcodes[0].code[Np]; + ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (ZP); } @@ -496,56 +837,4 @@ } -/* return calibrated magnitude from average/secfilt for given photcode */ -float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { - - int i, Ns; - float Mave, Mref, Mcol, mc; - double Mc; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mave = (Ns == -1) ? NAN : secfilt[Ns].M; - Mref = Mave + SCALE*code[0].C; - - mc = PhotColorForCode (average, secfilt, measure, code); - if (isnan(mc)) return (Mref); - mc = mc - SCALE*code[0].dX; - - Mc = mc; - Mcol = 0; - for (i = 0; i < code[0].Nc; i++) { - Mcol += code[0].X[i]*Mc; - Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ - } - Mref += Mcol; - return (Mref); -} - -/***/ -float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Mave; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mave = (Ns == -1) ? NAN : secfilt[Ns].M; - return (Mave); -} - -float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float dM; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - dM = (Ns == -1) ? NAN : secfilt[Ns].dM; - return (dM); -} - float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt) { @@ -558,40 +847,4 @@ Mstdev = (Ns == -1) ? NAN : 0.001*secfilt[Ns].Mstdev; return (Mstdev); -} - -float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Maper; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Maper = (Ns == -1) ? NAN : secfilt[Ns].Map; - return (Maper); -} - -float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Mkron; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mkron = (Ns == -1) ? NAN : secfilt[Ns].Mkron; - return (Mkron); -} - -float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float dMkron; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - dMkron = (Ns == -1) ? NAN : secfilt[Ns].dMkron; - return (dMkron); } Index: trunk/Ohana/src/libdvo/src/dvo_tiny_values.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 36528) +++ trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 36630) @@ -8,4 +8,6 @@ averageT[0].measureOffset = average[0].measureOffset; averageT[0].catID = average[0].catID; + averageT[0].objID = average[0].objID; + averageT[0].nOwn = -1; return; } @@ -33,4 +35,5 @@ measureT[0].dYccd = measure[0].dYccd; measureT[0].dRsys = measure[0].dRsys; + measureT[0].myDet = 0; return ; }