Repatch (from linux) of OSX IRAF

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