diff options
Diffstat (limited to 'noao/digiphot/photcal/mkobsfile/t_obsfile.x')
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/t_obsfile.x | 1131 |
1 files changed, 1131 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/mkobsfile/t_obsfile.x b/noao/digiphot/photcal/mkobsfile/t_obsfile.x new file mode 100644 index 00000000..6c997cea --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/t_obsfile.x @@ -0,0 +1,1131 @@ +include <fset.h> +include <ctotok.h> +include <ctype.h> +include "../lib/obsfile.h" + +# T_OBSFILE -- Make an observations catalog using the data extracted +# from the APPHOT/DAOPHOT catalogs by TXDUMP/TDUMP or supplied in a simple +# text file by the user, the observing parameters file, the shifts and +# aperture corrections file and the image sets file. OBSFILE expects to +# see the following 8 fields in the input data file: image name, xcenter, +# ycenter, exposure time, filter id, airmass, magnitude and magnitude error +# and uses a simple string containing a list of column numbers to determine +# which data is in which column. The observing parameters file is decoded +# in the same way. + +procedure t_obsfile () + +int photfiles # pointer to the list of photometry files +pointer columnstr # pointer to the list of input columns +pointer infilters # pointer to the input list of filter ids +int normtime # normalize to an exposure time of 1 time unit +int allfilters # only output objects with data in all filters +real tolerance # tolerance in pixels for position matching +int verify # verify interactive user input +int verbose # print status, warning and errors messages + +int incat, sets, sh, ap, obs, outcat, fmt, wrap +int csp, cp, nincolumns, nfilters, nimages, npts, sortimid +pointer sp, outfilters, fname, tfname, column, incolumns, obscolumns, imtable +pointer objid, x, y, mag, merr, imid, id +real minmagerr + +bool clgetb() +int clpopnu(), clgfil(), open(), clplen(), ph_filters(), ph_getlabels() +int ctoi(), ph_setimtable(), ph_esetimtable(), ph_mkimtable(), btoi() +int access() +pointer stopen() +real clgetr() + +begin + # Allocate working memory. + call smark (sp) + + call salloc (infilters, SZ_LINE, TY_CHAR) + call salloc (outfilters, SZ_LINE, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (tfname, SZ_FNAME, TY_CHAR) + + call salloc (columnstr, SZ_LINE, TY_CHAR) + call salloc (column, SZ_FNAME, TY_CHAR) + call salloc (incolumns, CAT_NFIELDS, TY_INT) + call salloc (obscolumns, OBS_NFIELDS, TY_INT) + + # Get the input file list. + photfiles = clpopnu ("photfiles") + if (clplen (photfiles) <= 0) + call error (0, "The input file list is empty") + + # Decode the input columns list. There must be at least CAT_NFIELDS + # column numbers in the incolumns list. The number 0 can be used as + # a place holder for missing columns. + + call clgstr ("incolumns", Memc[columnstr], SZ_LINE) + csp = 1 + nincolumns = 0 + while (ph_getlabels (Memc[columnstr], csp, Memc[column], SZ_FNAME) != + EOF) { + cp = 1 + if (ctoi (Memc[column], cp, Memi[incolumns+nincolumns]) <= 0) + Memi[incolumns+nincolumns] = 0 + nincolumns = nincolumns + 1 + } + if (nincolumns < CAT_NFIELDS) + call error (0, + "There are too few columns defined in the incolumns string") + + # Open the idfilters string. + call clgstr ("idfilters", Memc[infilters], SZ_LINE) + nfilters = ph_filters (Memc[infilters], Memc[outfilters], SZ_LINE) + + # Open the image sets file. + call clgstr ("imsets", Memc[fname], SZ_FNAME) + sets = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the output catalog file. + call clgstr ("observations", Memc[fname], SZ_FNAME) + outcat = open (Memc[fname], NEW_FILE, TEXT_FILE) + call sprintf (Memc[tfname], SZ_FNAME, "f%s.dat") + call pargstr (Memc[fname]) + if (access (Memc[tfname], 0, 0) == YES) + call delete (Memc[tfname]) + fmt = open (Memc[tfname], NEW_FILE, TEXT_FILE) + wrap = btoi (clgetb ("wrap")) + + # Get the minimum magnitude error. + minmagerr = clgetr ("minmagerr") + + # Normalize the exposure times? + normtime = btoi (clgetb ("normtime")) + + # Get the position matching parameters. + tolerance = clgetr ("tolerance") + allfilters = btoi (clgetb ("allfilters")) + + # Verify interactive user input. + verify = btoi (clgetb ("verify")) + verbose = btoi (clgetb ("verbose")) + + # Open the shifts file. + call clgstr ("shifts", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + sh = NULL + else + sh = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the aperture corrections file. + call clgstr ("apercors", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + ap = NULL + else + ap = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the observing parameters file and decode the column list. + + call clgstr ("obsparams", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + obs = NULL + else { + obs = open (Memc[fname], READ_ONLY, TEXT_FILE) + call amovki (0, Memi[obscolumns], OBS_NFIELDS) + if (obs != STDIN) { + call clgstr ("obscolumns", Memc[columnstr], SZ_LINE) + csp = 1 + nincolumns = 1 + Memi[obscolumns] = 1 + while (ph_getlabels (Memc[columnstr], csp, Memc[column], + SZ_FNAME) != EOF && (nincolumns < OBS_NFIELDS)) { + cp = 1 + if (ctoi (Memc[column], cp, + Memi[obscolumns+nincolumns]) <= 0) + Memi[obscolumns+nincolumns] = 0 + else if (Memi[obscolumns+nincolumns] == 1) + nincolumns = nincolumns - 1 + nincolumns = nincolumns + 1 + } + } + } + + # Open a symbol table to hold a list of image charactersitics. + + imtable = stopen ("imtable", 2 * LEN_IMTABLE, LEN_IMTABLE, 10 * + LEN_IMTABLE) + + # Read in the image set file and store each unique image name in + # the symbol table. Initialize the records fields. + + if (sets == NULL) + nimages = 0 + else if (sets == STDIN) + nimages = ph_esetimtable (imtable, nfilters, verify) + else + nimages = ph_setimtable (imtable, sets, nfilters, verbose) + + # For each image referenced in the observing parameters file extract + # the image name, the filter id if defined, the exposure time if + # defined and the airmass if defined. Enter the new values of filter + # id and airmass in the symbol table and compute the effective + # exposure time required to correct the instrumental magnitudes to + # the new exposure time. + + if (obs != NULL) { + if (obs == STDIN) + call ph_eobsimtable (imtable, nimages, verify) + else { + if (sets == STDIN && verbose == YES) { + call fstats (obs, F_FILENAME, Memc[fname], SZ_FNAME) + call eprintf ("Warning: Filter ids in %s ") + call pargstr (Memc[fname]) + call eprintf ("will replace those listed in idfilters\n") + } + call ph_obsimtable (imtable, obs, Memi[obscolumns], verbose) + } + } + + # Read in the x and y shifts for the image. Images for which no + # shift is defined are assigned an x-y shift of 0.0. + + if (sh != NULL) { + if (sh == STDIN) + call ph_eshimtable (imtable, nimages, verify) + else + call ph_shimtable (imtable, sh, verbose) + } + + # Read in the aperture corrections. Images for which no aperture + # correction is defined are assigned an aperture correction of 0.0. + + if (ap != NULL) { + if (ap == STDIN) + call ph_eapimtable (imtable, nimages, verify) + else + call ph_apimtable (imtable, ap, verbose) + } + + # Read in the data. Extract each unique image name and match it + # with the corresponding image name in the symbol table. Skip data + # for images which are not in the symbol table. For each unique image + # name which does match a symbol table image name, extract the exposure + # time, filter id, and airmass and enter them in the symbol table. + # Read the x and y and magnitude and magnitude error data into + # separate data arrays. + + objid = NULL + x = NULL + y = NULL + mag = NULL + merr = NULL + imid = NULL + id = NULL + sortimid = NO + + npts = 0 + while (clgfil (photfiles, Memc[fname], SZ_FNAME) != EOF) { + incat = open (Memc[fname], READ_ONLY, TEXT_FILE) + npts = ph_mkimtable (imtable, incat, Memi[incolumns], objid, x, + y, mag, merr, imid, id, npts, normtime, sortimid, verbose) + call close (incat) + } + + if (objid != NULL) + call realloc (objid, npts * (DEF_LENLABEL+1), TY_CHAR) + call realloc (x, npts, TY_REAL) + call realloc (y, npts, TY_REAL) + call realloc (mag, npts, TY_REAL) + call realloc (merr, npts, TY_REAL) + call realloc (imid, npts, TY_INT) + call realloc (id, npts, TY_INT) + + # Sort the data by image id and object id if position matching is + # turned off or by image id and y coordinate if position matching is + # turned on. + + if (objid == NULL) + call ph_sort (imtable, nimages, "", Memr[x], Memr[y], Memr[mag], + Memr[merr], Memi[imid], Memi[id], npts, sortimid, tolerance) + else + call ph_sort (imtable, nimages, Memc[objid], Memr[x], Memr[y], + Memr[mag], Memr[merr], Memi[imid], Memi[id], npts, sortimid, + tolerance) + # Free the image id sorting index arrays. + + if (imid != NULL) + call mfree (imid, TY_INT) + + # Apply the exposure time corrections, x-y shifts, and the aperture + # corrections. + + call ph_correct (imtable, Memr[x], Memr[y], Memr[mag], Memr[merr], + minmagerr) + + if (verbose == YES) { + call fstats (outcat, F_FILENAME, Memc[fname], SZ_FNAME) + call printf ("\nObservations file: %s\n") + call pargstr (Memc[fname]) + } + + # Match by order in the catalog if tolerance is less than or equal to + # zero or by x-y position if tolerance is greater than zero and output + # the data. + + if (objid == NULL) + call ph_match (imtable, nimages, Memc[outfilters], nfilters, outcat, + "", Memr[x], Memr[y], Memr[mag], Memr[merr], Memi[id], npts, + tolerance, allfilters, wrap, verbose) + else + call ph_match (imtable, nimages, Memc[outfilters], nfilters, outcat, + Memc[objid], Memr[x], Memr[y], Memr[mag], Memr[merr], + Memi[id], npts, tolerance, allfilters, wrap, verbose) + + if (verbose == YES) + call printf ("\n") + + # Write the format file for the output catalog. + + call ph_format (fmt, Memc[outfilters], nfilters) + + # Free the data arrays. + + if (objid != NULL) + call mfree (objid, TY_CHAR) + if (x != NULL) + call mfree (x, TY_REAL) + if (y != NULL) + call mfree (y, TY_REAL) + if (mag != NULL) + call mfree (mag, TY_REAL) + if (merr != NULL) + call mfree (merr, TY_REAL) + if (id != NULL) + call mfree (id, TY_INT) + + # Close the symbol table. + call stclose (imtable) + + # Close the files and file lists. + if (sets != NULL) + call close (sets) + if (sh != NULL) + call close (sh) + if (ap != NULL) + call close (ap) + if (obs != NULL) + call close (obs) + call close (outcat) + call close (fmt) + call clpcls (photfiles) + + call sfree (sp) +end + + +# PH_SORT -- Sort the data. + +procedure ph_sort (imtable, nimages, objid, x, y, mag, merr, imid, id, npts, + sortimid, tolerance) + +pointer imtable # pointer to the symbol table +int nimages # the number of images in the symbol table +char objid[ARB] # the array of object ids objid[1] = EOS if none +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real merr[ARB] # array of magnitude errors +int imid[ARB] # array of image ids +int id[ARB] # array of object ids +int npts # the number of points +int sortimid # sort by image and id +real tolerance # the tolerance in pixels + +int i, nptr +pointer sp, sym, symbol +pointer sthead(), stnext() + +begin + # Allocate working space. + call smark (sp) + call salloc (sym, nimages, TY_INT) + + # Read in the symbol list in the last in first out sense. + + symbol = sthead (imtable) + do i = 1, nimages { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Do the image/object id sort if necessary. This is a real sort + # where the actual contents of the arrays are rearranged. + # After the image/object id sort recompute the offsets specifying the + # location of the image data in the symbol table. The first + # pass through the symbol table to picks up the information for the + # images that have only one entry. The second pass picks up + # information for images that have more than one entry. + + if (sortimid == YES) { + if (objid[1] == EOS) + call ph_4r2isort (x, y, mag, merr, imid, id, npts) + else + call ph_1c4r2isort (objid, x, y, mag, merr, imid, id, npts) + nptr = npts + 1 + do i = 1, nimages { + symbol = Memi[sym+i-1] + if (IMT_NENTRIES(symbol) <= 0) + next + nptr = nptr - IMT_NENTRIES(symbol) + IMT_OFFSET(symbol) = nptr + } + } + + # Sort on y if necessary. This is an index sort and does not alter + # the position of the data. + + if (tolerance > 0) { + do i = 1, nimages { + symbol = Memi[sym+i-1] + if (IMT_NENTRIES(symbol) <= 0) + next + nptr = IMT_OFFSET(symbol) + call ph_qsort (y[nptr], id[nptr], IMT_NENTRIES(symbol), nptr) + } + } + + call sfree (sp) +end + + +# PH_CORRECT -- Correct the x, y and mag values using the x-y shifts, aperture +# corrections, and exposure times. + +procedure ph_correct (imtable, x, y, mag, merr, minmagerr) + +pointer imtable # pointer to the symbol table +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real merr[ARB] # array of magnitude errors +real minmagerr # the minimum magnitude error + +int i, nptr, ndata +pointer sp, imname, symbol, osymbol +real magshift +pointer sthead(), stnext(), stfind() + +begin + # Allocate working space. + + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + # Loop through the images adding the x, y and magnitude shifts to + # the data. + + symbol = sthead (imtable) + while (symbol != NULL) { + + # Get the correct symbol. + ndata = IMT_NENTRIES(symbol) + osymbol = stfind (imtable, IMT_IMNAME(symbol)) + + # If data is present make the corrections once for each image. + + if ((ndata > 0) && (osymbol == symbol)) { + nptr = IMT_OFFSET(symbol) + + # Correct the positions. + call aaddkr (x[nptr], IMT_XSHIFT(symbol), x[nptr], ndata) + call aaddkr (y[nptr], IMT_YSHIFT(symbol), y[nptr], ndata) + + # Correct the magnitudes. + magshift = 2.5 * log10 (IMT_ITIME(symbol)) + IMT_APERCOR(symbol) + do i = nptr, nptr + ndata - 1 { + if (IS_INDEFR(mag[i])) + next + mag[i] = mag[i] + magshift + if (IS_INDEFR(merr[i])) + next + if (merr[i] < minmagerr) + merr[i] = minmagerr + } + } + + # Get the next symbol. + symbol = stnext (imtable, symbol) + } + + call sfree (sp) +end + + +# PH_MATCH -- Match up the photometry lists by order in the input catalog +# or x-y position. + +procedure ph_match (imtable, nimages, filters, nfilters, out, objid, x, y, + mag, merr, ysort, npts, tolerance, allfilters, wrap, verbose) + +pointer imtable # pointer to the image name symbol table +int nimages # the number of images in the symbol table +char filters[ARB] # the filters string +int nfilters # the number of filters +int out # the output file descriptor +char objid[ARB] # the object id array +real x[ARB] # the x input array +real y[ARB] # the y input array +real mag[ARB] # the magnitude array +real merr[ARB] # the magnitude error array +int ysort[ARB] # the y sort index +int npts # the number of points +real tolerance # the x-y matching tolerance in pixels +int allfilters # match in all filters +int wrap # format the output file for easy reading +int verbose # print status, warning and error messages + +int i, k, nsets, imptr, filterno, ndata +pointer sp, imsym, filterid, label, sym, match, osymbol, symbol, index +bool streq() +int strdic(), ph_nthlabel() +pointer sthead(), stnext(), stfind() + +begin + # Write out the output file column headers. + call fprintf (out, "\n") + if (wrap == YES) { + call fprintf (out, +"# FIELD%17tFILTER%34tOTIME%40tAIRMASS%49tXCENTER%59tYCENTER%71tMAG%77tMERR\n") + } else { + do i = 1, nfilters { + if (i == 1) + call fprintf (out, +"# FIELD%17tFILTER%34tOTIME%40tAIRMASS%49tXCENTER%59tYCENTER%71tMAG%77tMERR ") + else + call fprintf (out, +"%2tFILTER%19tOTIME%25tAIRMASS%34tXCENTER%44tYCENTER%56tMAG%62tMERR ") + } + call fprintf (out,"\n") + } + call fprintf (out, "\n") + + # Allocate and or initialize working space. + call smark (sp) + call salloc (imsym, nimages, TY_INT) + call salloc (filterid, SZ_FNAME, TY_CHAR) + call salloc (label, SZ_FNAME, TY_CHAR) + call salloc (sym, nfilters, TY_INT) + if (tolerance > 0.0) + call malloc (match, npts, TY_INT) + else + match = NULL + index = NULL + + + # Read in the symbol table pointers and rearrange them so that the + # image name symbols are ordered in the same way as the image set + # numbers. Compute the total number of image sets at the same time. + + nsets = 0 + symbol = sthead (imtable) + do i = nimages, 1, -1 { + nsets = max (nsets, IMT_IMSETNO(symbol)) + Memi[imsym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Loop over the image sets. + imptr = 1 + do i = 1, nsets { + + # Clear the array of symbols and set the label. + ndata = 0 + call amovki (NULL, Memi[sym], nfilters) + call strcpy (IMT_LABEL(Memi[imsym+imptr-1]), Memc[label], SZ_FNAME) + + # Loop over the images in the set. + repeat { + + # Quit if the image pointer exceeds the number of images. + if (imptr > nimages) + break + + # Get the next symbol for a given image set. + symbol = Memi[imsym+imptr-1] + if (IMT_IMSETNO(symbol) != i) + break + + # Check for duplicate image names. + osymbol = stfind (imtable, IMT_IMNAME(symbol)) + if (osymbol != symbol) { + IMT_OFFSET(symbol) = IMT_OFFSET(osymbol) + IMT_NENTRIES(symbol) = IMT_NENTRIES(osymbol) + IMT_XAIRMASS(symbol) = IMT_XAIRMASS(osymbol) + IMT_OTIME(symbol) = IMT_OTIME(osymbol) + call strcpy (IMT_IFILTER(osymbol), IMT_IFILTER(symbol), + SZ_FNAME) + } + + # Search for the filter id, if one is not found then + # use the first unoccupied space. + + filterno = strdic (IMT_IFILTER(symbol), Memc[filterid], + SZ_FNAME, filters) + if (streq (IMT_IFILTER(symbol), "INDEF")) { + do k = 1, nfilters { + if (Memi[sym+k-1] != NULL) + next + if (ph_nthlabel (filters, k, Memc[filterid], + SZ_FNAME) != k) + Memc[filterid] = EOS + filterno = k + break + } + } + + # Load the data pointers. + if (filterno <= 0) { + + if (verbose == YES) { + call eprintf ( + "Warning: Image set: %d label: %s image: %s\n") + call pargi (i) + call pargstr (IMT_LABEL(symbol)) + call pargstr (IMT_IMNAME(symbol)) + call eprintf ( + "\tWarning: Filter %s ") + call pargstr (IMT_IFILTER(symbol)) + call eprintf ("does not belong to filter set %s\n") + call pargstr (filters) + } + + } else if (Memi[sym+filterno-1] != NULL) { + + if (verbose == YES) { + call eprintf ( + "Warning: Image set: %d label: %s image: %s\n") + call pargi (i) + call pargstr (IMT_LABEL(symbol)) + call pargstr (IMT_IMNAME(symbol)) + call eprintf ("\tData for filter %s is redundant\n") + call pargstr (Memc[filterid]) + } + + } else { + + Memi[sym+filterno-1] = symbol + ndata = max (ndata, IMT_NENTRIES(symbol)) + } + + imptr = imptr + 1 + + } + + # Skip matching if there are no points. + if (ndata <= 0) { + if (verbose == YES) { + call eprintf ("Image set: %d Label: %s has no data\n") + call pargi (i) + call pargstr (Memc[label]) + } + next + } + + # Do the matching. + if (tolerance <= 0.0) { + if (objid[1] == EOS) + call ph_join (out, Memc[label], Memi[sym], filters, + nfilters, "", x, y, mag, merr, allfilters, wrap, + verbose) + else + call ph_join (out, Memc[label], Memi[sym], filters, + nfilters, objid, x, y, mag, merr, allfilters, wrap, + verbose) + } else { + if (index == NULL) + call malloc (index, nfilters * ndata, TY_INT) + else + call realloc (index, nfilters * ndata, TY_INT) + if (objid[1] == EOS) + call ph_merge (out, Memc[label], Memi[sym], filters, + nfilters, "", x, y, mag, merr, ysort, Memi[match], + Memi[index], ndata, tolerance, allfilters, wrap, + verbose) + else + call ph_merge (out, Memc[label], Memi[sym], filters, + nfilters, objid, x, y, mag, merr, ysort, Memi[match], + Memi[index], ndata, tolerance, allfilters, wrap, + verbose) + } + } + + # Free working space. + if (match != NULL) + call mfree (match, TY_INT) + if (index != NULL) + call mfree (index, TY_INT) + call sfree (sp) +end + + +# PH_FORMAT -- Write the format file for the output catalog. + +procedure ph_format (fd, filters, nfilters) + +int fd # file descripter for the format file +char filters[ARB] # list of filter ids +int nfilters # number of filters + +int i, col +pointer sp, filterid +int ph_nthlabel() + +begin + call smark (sp) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + call fprintf (fd, "# Declare the observations file variables\n\n") + call fprintf (fd, "observations\n\n") + + col = FIRST_COLUMN + do i = 1, nfilters { + if (ph_nthlabel (filters, i, Memc[filterid], SZ_FNAME) != i) + Memc[filterid] = EOS + call fprintf (fd, + "T%s%15t%d%30t# time of observation in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "X%s%15t%d%30t# airmass in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "x%s%15t%d%30t# x coordinate in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "y%s%15t%d%30t# y coordinate in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, + "m%s%15t%d%30t# instrumental magnitude in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, + "error(m%s)%15t%d%30t# magnitude error in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + call fprintf (fd, "\n") + col = col + DELTA_COLUMN + } + + call sfree (sp) +end + + +# PH_NXTIMAGE -- Find the first line in the input catalog containing the next +# image name, given the current image name. Return the new image name and line. + +int procedure ph_nxtimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # the line buffer size + +int i, op, colno +long stat +pointer sp, str +bool streq() +int getline(), nscan() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + Memc[str] = EOS + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() != colno) + next + + } until (! streq (image, Memc[str])) + call strcpy (Memc[str], image, SZ_FNAME) + + call sfree (sp) + + return (stat) +end + + +# PH_GETIMAGE -- Read the next line in the input catalog and return the image +# name and the line. + +int procedure ph_getimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # the size of the output buffer + +int op, i, colno +int stat +int getline(), nscan() + +begin + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (image, SZ_FNAME) + if (nscan() != colno) + next + + break + } + + return (stat) +end + + +# PH_IMDATA -- Decode the filter id, integration time and airmass from +# a line of the input catalog. + +procedure ph_imdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the name of the filter +int sz_filterid # maximum length of the filter id +real itime # the exposure time of the observation +real airmass # the airmass of the observation +real otime # the time of observations + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[CAT_IFILTER] + icol = columns[CAT_ITIME] + acol = columns[CAT_XAIRMASS] + ocol = columns[CAT_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + call strcpy ("INDEF", filterid, sz_filterid) + itime = INDEFR + airmass = INDEFR + otime = INDEFR + + call sscan (line) + do i = 1, maxcol { + if ( i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_OBSDATA -- Decode the filter id, integration time and airmass from +# a line of the observing parameters file. + +procedure ph_obsdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the name of the filter +int sz_filterid # maximum length of the filter id +real itime # the exposure time of the observation +real airmass # the airmass of the observation +real otime # the time of observation + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[OBS_IFILTER] + icol = columns[OBS_ITIME] + acol = columns[OBS_XAIRMASS] + ocol = columns[OBS_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + call strcpy ("INDEF", filterid, sz_filterid) + itime = INDEFR + airmass = INDEFR + otime = INDEFR + + call sscan (line) + do i = 1, maxcol { + if ( i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_STARDATA -- Decode the image name, x and y coordinates, magnitude +# and magnitude error from the input catalog. + +int procedure ph_stardata (line, columns, objid, x, y, mag, merr) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char objid[ARB] # the object id +real x, y # the output x and y coordinates +real mag, merr # the output magnitude and magnitude error + +int i, icol, xcol, ycol, mcol, ecol, maxcol, stat +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + icol = columns[CAT_ID] + xcol= columns[CAT_XCENTER] + ycol= columns[CAT_YCENTER] + mcol = columns[CAT_MAG] + ecol = columns[CAT_MERR] + maxcol = max (icol, xcol, ycol, mcol, ecol) + + objid[1] = EOS + x = INDEFR + y = INDEFR + mag = INDEFR + merr = INDEFR + stat = OK + + call sscan (line) + do i = 1, maxcol { + + if (i == icol) { + call gargwrd (objid, DEF_LENLABEL) + if (nscan() != icol) + stat = ERR + } else if (i == xcol) { + call gargr (x) + if (nscan() != xcol) + stat = ERR + } else if (i == ycol) { + call gargr (y) + if (nscan() != ycol) + stat = ERR + } else if (i == mcol) { + call gargr (mag) + if (nscan() != mcol) + stat = ERR + } else if (i == ecol) { + call gargr (merr) + if (nscan() != ecol) + stat = ERR + } else { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() != i) + stat = ERR + } + } + + call sfree (sp) + + return (stat) +end + + +# PH_FILTERS -- Reformat the filter string so that it is a suitable string +# dictionary for the STRDIC routine. + +int procedure ph_filters (infilters, outfilters, max_lenfilters) + +char infilters[ARB] # input list of filters +char outfilters[ARB] # output list of filters +int max_lenfilters # maximum length of the filter string + +int ip, nfilters +pointer sp, filter +int ph_getlabels() + +begin + call smark (sp) + call salloc (filter, max_lenfilters, TY_CHAR) + + ip = 1 + nfilters = 0 + outfilters[1] = EOS + while (ph_getlabels (infilters, ip, Memc[filter], max_lenfilters) != + EOF) { + call strcat (",", outfilters, max_lenfilters) + call strcat (Memc[filter], outfilters, max_lenfilters) + nfilters = nfilters + 1 + } + + call sfree (sp) + + return (nfilters) +end + + +# PH_NTHLABEL -- Get the nth label out of a list of labels. + +int procedure ph_nthlabel (list, item, label, maxch) + +char list[ARB] # input list +int item # item to be extracted +char label[ARB] # extracted label +int maxch # maximum length of the extracted label + +int ip, nitems +int ph_getlabels() + +begin + nitems = 0 + + ip = 1 + while (ph_getlabels (list, ip, label, maxch) != EOF) { + nitems = nitems + 1 + if (nitems == item) + break + } + + return (nitems) +end + + +# PH_GETLABELS -- Get the next label from a list of labels. + +int procedure ph_getlabels (list, ip, label, maxch) + +char list[ARB] # list of labels +int ip # pointer in to the list of labels +char label[ARB] # the output label +int maxch # maximum length of a column name + +int op, token +int ctotok(), strlen() + +begin + # Decode the column labels. + op = 1 + while (list[ip] != EOS) { + + token = ctotok (list, ip, label[op], maxch) + if (label[op] == EOS) + next + if ((token == TOK_UNKNOWN) || (token == TOK_CHARCON)) + break + if ((token == TOK_PUNCTUATION) && (label[op] == ',')) { + if (op == 1) + next + else + break + } + + op = op + strlen (label[op]) + if (IS_WHITE(list[ip])) + break + } + + label[op] = EOS + if ((list[ip] == EOS) && (op == 1)) + return (EOF) + else + return (op - 1) +end |