Initial commit

8 files changed, 3248 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/nstar/dpggroup.x b/noao/digiphot/daophot/nstar/dpggroup.x
new file mode 100644
index 00000000..fa180c17
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpggroup.x
@@ -0,0 +1,386 @@
+include "../../lib/ptkeysdef.h"
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+
+# DP_GNSTPSF -- Procedure to initialize for reading the group file fields from
+# a photometry text file . The group file fields are ID, GROUP, X, Y, MAG, ERR,
+# and SKY.
+
+procedure dp_gnstpsf (fields, sel_fields, max_nfields)
+
+int	fields[ARB]		# array of selected fields
+char	sel_fields[ARB]		# names of selected containing fields
+int	max_nfields		# maximum number of fields selected
+
+int	i
+int	strlen()
+
+begin
+	# Initialize the fields string.
+	sel_fields[1] = EOS
+
+	# Encode the fields.
+	do i = 1, max_nfields {
+	    switch (fields[i]) {
+	    case DP_PAPID:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	            call pargstr (ID)
+	    case DP_PAPXCEN:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	    	    call pargstr (XCENTER)
+	    case DP_PAPYCEN:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	    	    call pargstr (YCENTER)
+	    case DP_PAPSKY:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	    	    call pargstr (SKY)
+	    case DP_PAPMAG1:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	            call pargstr (MAG)
+	    case DP_PAPGROUP:
+		call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ")
+	    	    call pargstr (GROUP)
+	    }
+	}
+
+	# Backspace over the terminating blank character.
+	if (sel_fields[1] != EOS)
+	    sel_fields[strlen(sel_fields)] = EOS
+end
+
+
+# DP_TNSTINIT -- Procedure to initialize for reading the group file fields from
+# a photometry table. The group file fields are ID, GROUP, X, Y, MAG, ERR,
+# and SKY.
+
+procedure dp_tnstinit (tp, colpoint)
+
+pointer	tp			# the table descriptor
+pointer	colpoint[ARB]		# the column descriptor
+
+begin
+	# Get the id.
+	# First the ID
+	call tbcfnd (tp, ID, colpoint[1], 1)
+	if (colpoint[1] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (ID)
+	}
+
+	# Get the x position.
+	call tbcfnd (tp, XCENTER, colpoint[2], 1)
+	if (colpoint[2] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (XCENTER)
+	}
+
+	# Get the y position.
+	call tbcfnd (tp, YCENTER, colpoint[3], 1)
+	if (colpoint[3] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (YCENTER)
+	}
+
+	# Get the magnitude.
+	call tbcfnd (tp, MAG, colpoint[4], 1)
+	if (colpoint[4] == NULL)		# No column
+	    call tbcfnd (tp, APMAG, colpoint[4], 1)
+	if (colpoint[4] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	    call pargstr (APMAG)
+	}
+
+	# Get the sky.
+	call tbcfnd (tp, SKY, colpoint[5], 1)
+	if (colpoint[5] == NULL)
+	    call tbcfnd (tp, SKY, colpoint[5], 1)
+	if (colpoint[5] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (SKY)
+	}
+
+	# Get the group number.
+	call tbcfnd (tp, GROUP, colpoint[6], 1)
+	if (colpoint[6] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (GROUP)
+	}
+end
+
+
+# DP_GGROUP -- Read in a single group.
+
+int procedure dp_ggroup (dao, tp, key, fields, indices, colpoint, max_row,
+	max_group, in_record, curr_group)
+
+pointer	dao				# pointer to daophot structure
+int	tp				# input file/table descriptor
+pointer	key				# pointer to text database structure
+char	fields[ARB]			# nstar fields to be read
+int	indices[ARB]			# array of text file field pointers
+int	colpoint[ARB]			# array of column pointers
+int	max_row				# number of rows in table
+int	max_group			# maximum group size
+int	in_record			# pointer to current input record
+int	curr_group			# current group number
+
+bool	nullflag
+int	istar, group, buf_size
+pointer	apsel
+int	dp_nstsel()
+
+begin
+	# If the current input record is set to zero we are at EOF. In_record
+	# is initialized to one on entry to this routine.
+
+	if (in_record == 0)
+	    return (0)
+
+	# Get the next group number. Note that the last star read on the
+	# previous call is the first star in the new group.
+
+	apsel = DP_APSEL(dao)
+	if (in_record == 1) {
+	    buf_size = max_group
+	    group = 0
+	    istar = 0
+	} else {
+	    Memi[DP_APID(apsel)] = Memi[DP_APID(apsel)+istar-1]
+	    Memr[DP_APXCEN(apsel)] = Memr[DP_APXCEN(apsel)+istar-1]
+	    Memr[DP_APYCEN(apsel)] = Memr[DP_APYCEN(apsel)+istar-1]
+	    Memr[DP_APMAG(apsel)] = Memr[DP_APMAG(apsel)+istar-1]
+	    Memr[DP_APMSKY(apsel)] = Memr[DP_APMSKY(apsel)+istar-1]
+	    istar = 1
+	}
+	
+	# Loop over the stars in a single group.
+	repeat {
+
+	    # Set the current group.
+	    curr_group = group
+
+	    # Read in the photometry for a single star.
+
+	    # In this case we have a text database file.
+	    if (key != NULL) {
+
+	        if (dp_nstsel (key, tp, fields, indices, Memi[DP_APID(apsel)+
+		    istar], group, Memr[DP_APXCEN(apsel)+istar],
+		    Memr[DP_APYCEN(apsel)+istar], Memr[DP_APMSKY(apsel)+
+		    istar], Memr[DP_APMAG(apsel)+ istar]) == EOF) {
+		    in_record = 0
+		    break
+		}
+
+	    # In this case we have a table.
+	    } else {
+
+		if (in_record > max_row) {
+		    in_record = 0
+		    break
+		} else {
+		    call tbrgti (tp, colpoint[1], Memi[DP_APID(apsel)+istar],
+		        nullflag, 1, in_record)
+		    call tbrgtr (tp, colpoint[2], Memr[DP_APXCEN(apsel)+istar],
+		        nullflag, 1, in_record)
+		    call tbrgtr (tp, colpoint[3], Memr[DP_APYCEN(apsel)+istar],
+		        nullflag, 1, in_record)
+		    call tbrgtr (tp, colpoint[4], Memr[DP_APMAG(apsel)+istar],
+		        nullflag, 1, in_record)
+		    call tbrgtr (tp, colpoint[5], Memr[DP_APMSKY(apsel)+istar],
+		        nullflag, 1, in_record)
+	            call tbrgti (tp, colpoint[6], group, nullflag, 1, in_record)
+		}
+	    }
+
+	    # Increment the record and star counters.
+	    in_record = in_record + 1
+	    istar = istar + 1
+
+	    # Allocate more memory as needed.
+	    if (istar == buf_size) {
+		buf_size = buf_size + max_group
+		call dp_rmemapsel (dao, indices, NAPPAR, buf_size + 1)
+	    }
+
+	} until ((group != curr_group) && (curr_group != 0))
+
+	# Return the number of stars in the group.
+	if (in_record == 0) {
+	    if (curr_group == 0)
+		curr_group = group
+	    return (istar)
+	} else
+	    return (istar - 1)
+end
+
+
+# DP_NSTSEL --  Read in the required photometry records from a text file.
+
+int procedure dp_nstsel (key, fd, fields, indices, id, group, x, y, sky, mag)
+
+pointer key		# pointer to key structure
+int	fd		# text file descriptor
+char	fields[ARB]	# fields to be output
+int	indices[ARB]	# indices array
+int	id		# star id number
+int	group		# group number
+real	x		# x center
+real	y		# y center
+real	sky		# sky value
+real	mag		# magnitude
+
+int	nchars, nunique, uunique, funique, ncontinue, recptr
+int 	first_rec, nselect, record
+pointer	line
+int	getline(), strncmp(), pt_choose()
+data	first_rec /YES/
+
+begin
+	# Initialize the file read.
+	if (first_rec == YES) {
+	    nunique = 0
+	    uunique = 0
+	    funique = 0
+	    nselect = 0
+	    record = 0
+	    call malloc (line, SZ_LINE, TY_CHAR)
+	}
+
+	# Initialize the record read.
+	ncontinue = 0
+	recptr = 1
+
+	# Loop over the text file records.
+	repeat  {
+
+	    # Read in a line of the text file.
+	    nchars = getline (fd, Memc[line])
+	    if (nchars == EOF)
+		break
+
+	    # Determine the type of record.
+	    if (Memc[line] == KY_CHAR_POUND) {
+
+	        if (strncmp (Memc[line], KY_CHAR_KEYWORD, KY_LEN_STR) == 0) {
+		    call pt_kyadd (key, Memc[line], nchars)
+	        } else if (strncmp (Memc[line], KY_CHAR_NAME,
+		    KY_LEN_STR) == 0) {
+		    nunique = nunique + 1
+		    call pt_kname (key, Memc[line], nchars, nunique)
+	        } else if (strncmp (Memc[line], KY_CHAR_UNITS,
+		    KY_LEN_STR) == 0) {
+		    uunique = uunique + 1
+		    call pt_knunits (key, Memc[line], nchars, uunique)
+	        } else if (strncmp (Memc[line], KY_CHAR_FORMAT,
+		    KY_LEN_STR) == 0) {
+		    funique = funique + 1
+		    call pt_knformats (key, Memc[line], nchars, funique)
+	        }
+
+	    } else if (Memc[line] == KY_CHAR_NEWLINE) {
+		# skip blank lines
+
+	    } else {
+
+		# Construct the text file record.
+		call pt_mkrec (key, Memc[line], nchars, first_rec, recptr,
+		    ncontinue) 
+
+	        # Construct output record when there is no continuation char.
+	        if (Memc[line+nchars-2] != KY_CHAR_CONT) {
+
+		    # Select the appropriate records.
+		    if (nselect <= 0) {
+		        nselect = pt_choose (key, fields)
+			if (nselect < NAPGROUP) {
+			    call eprintf (
+			    "The group file does not have the correct format\n")
+			    break
+			}
+		    }
+
+		    # Construct the output record by moving the selected fields
+		    # into the data structures.
+
+		    call dp_gnst (key, indices, id, group, x, y, sky, mag)
+		    record = record + 1
+		    first_rec = NO
+
+		    # Record is complete so exit the loop.
+		    break
+	        }
+	    }
+
+	}
+
+	# Return EOF or the record number.
+	if (nchars == EOF || (nselect < NAPGROUP)) {
+	    first_rec = YES
+	    nunique = 0
+	    uunique = 0
+	    funique = 0
+	    nselect = 0
+	    call mfree (line, TY_CHAR)
+	    return (EOF)
+	} else
+	    return (record)
+end
+
+
+# DP_GNST  -- Decode the standard GROUP text file fields into the appropriate
+# arrays.
+
+procedure dp_gnst (key, fields, id, group, x, y, sky, mag)
+
+pointer	key		# pointer to keys strucuture
+int	fields[ARB]	# fields array
+int	id		# star id
+int	group		# group id
+real	x		# x position
+real	y		# y position
+real	sky		# sky value
+real	mag		# magnitude
+
+int	i, index, elem, maxch, kip, ip
+int	ctoi(), ctor()
+char	buffer[SZ_LINE]
+
+begin
+	do i = 1, KY_NSELECT(key) {
+
+	    # Find the key.
+	    index = Memi[KY_SELECT(key)+i-1]
+	    elem = Memi[KY_ELEM_SELECT(key)+i-1]
+	    maxch = Memi[KY_LEN_SELECT(key)+i-1]
+	    kip = Memi[KY_PTRS(key)+index-1] + (elem - 1) * maxch
+	    call amovc (Memc[kip], buffer, maxch)
+	    buffer[maxch+1] = EOS
+
+	    # Decode the output value.
+	    ip = 1
+	    switch (fields[i]) {
+	    case DP_PAPID: 
+		if (ctoi (buffer, ip, id) <= 0)
+		    call error (0, "ERROR: Error reading ID field.")
+	    case DP_PAPGROUP: 
+		if (ctoi (buffer, ip, group) <= 0)
+		    call error (0, "ERROR: Error reading GROUP field.")
+	    case DP_PAPXCEN:
+		if (ctor (buffer, ip, x) <= 0)
+		    call error (0, "ERROR: Error reading XCENTER field.")
+	    case DP_PAPYCEN:
+		if (ctor (buffer, ip, y) <= 0)
+		    call error (0, "ERROR: Error reading YCENTER field.")
+	    case DP_PAPSKY:
+	        if (ctor (buffer, ip, sky) <= 0)
+		    call error (0, "ERROR: Error reading MSKY field.")
+	    case DP_PAPMAG1:
+		if (ctor (buffer, ip, mag) <= 0)
+		    call error (0, "ERROR: Error reading MAG field.")
+	    default:
+		call printf ("Error reading the photometry file.\n")
+	    }
+
+	}
+end
diff --git a/noao/digiphot/daophot/nstar/dpmemnstar.x b/noao/digiphot/daophot/nstar/dpmemnstar.x
new file mode 100644
index 00000000..04bd25f4
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpmemnstar.x
@@ -0,0 +1,158 @@
+include "../lib/daophotdef.h"
+include "../lib/nstardef.h"
+
+# DP_NSTARSETUP -- Procedure to set up the NSTAR parameters.
+
+procedure dp_nstarsetup (dp)
+
+pointer	dp		# pointer to daophot structure
+
+pointer	nstar
+
+begin
+	# Allocate Memory
+	call malloc (DP_NSTAR(dp), LEN_NSTARSTRUCT, TY_STRUCT)
+	nstar = DP_NSTAR(dp)
+
+	DP_NNPIX(nstar) = NULL
+	DP_NNUMER(nstar) = NULL
+	DP_NDENOM(nstar) = NULL
+	DP_NRPIXSQ(nstar) = NULL
+	DP_NSKIP(nstar) = NULL
+	DP_NXCLAMP(nstar) = NULL
+	DP_NXOLD(nstar) = NULL
+	DP_NX(nstar) = NULL
+	DP_NV(nstar) = NULL
+	DP_NSUMWT(nstar) = NULL
+	DP_NC(nstar) = NULL
+	DP_NIER(nstar) = NULL
+end
+
+
+# DP_MEMNSTAR -- Procedure to allocate sufficient memory for NSTAR.
+
+procedure dp_memnstar (dao, max_star)
+
+pointer	dao			# pointer to daophot structure
+int	max_star		# maximum number of stars
+
+pointer	nstar
+
+begin
+	nstar = DP_NSTAR(dao)
+
+	if (DP_NNPIX (nstar) != NULL)
+	    call mfree (DP_NNPIX (nstar), TY_INT)
+	call malloc (DP_NNPIX (nstar), max_star + 1, TY_INT)
+
+	if (DP_NNUMER (nstar) != NULL)
+	    call mfree (DP_NNUMER (nstar), TY_REAL)
+	call malloc (DP_NNUMER (nstar), max_star + 1, TY_REAL)
+
+	if (DP_NDENOM (nstar) != NULL)
+	    call mfree (DP_NDENOM (nstar), TY_REAL)
+	call malloc (DP_NDENOM (nstar), max_star + 1, TY_REAL)
+
+	if (DP_NRPIXSQ (nstar) != NULL)
+	    call mfree (DP_NRPIXSQ (nstar), TY_REAL)
+	call malloc (DP_NRPIXSQ (nstar), max_star + 1, TY_REAL)
+
+	if (DP_NSKIP (nstar) != NULL)
+	    call mfree (DP_NSKIP (nstar), TY_INT)
+	call malloc (DP_NSKIP (nstar), max_star + 1, TY_INT)
+
+	if (DP_NIER (nstar) != NULL)
+	    call mfree (DP_NIER (nstar), TY_INT)
+	call malloc (DP_NIER(nstar), max_star + 1, TY_INT)
+
+	if (DP_NSUMWT (nstar) != NULL)
+	    call mfree (DP_NSUMWT (nstar), TY_REAL)
+	call malloc (DP_NSUMWT (nstar), max_star + 1, TY_REAL)
+
+	if (DP_RECENTER(dao) == YES) {
+
+	    if (DP_NXCLAMP (nstar) != NULL)
+	        call mfree (DP_NXCLAMP (nstar), TY_REAL)
+	    call malloc (DP_NXCLAMP (nstar), 3 * max_star + 1, TY_REAL)
+
+	    if (DP_NXOLD (nstar) != NULL)
+	        call mfree (DP_NXOLD (nstar), TY_REAL)
+	    call malloc (DP_NXOLD (nstar), 3 * max_star + 1, TY_REAL)
+
+	    if (DP_NX (nstar) != NULL)
+	        call mfree (DP_NX (nstar), TY_REAL)
+	    call malloc (DP_NX (nstar), 3 * max_star + 1, TY_REAL)
+
+	    if (DP_NC (nstar) != NULL)
+	        call mfree (DP_NC (nstar), TY_REAL)
+	    call malloc (DP_NC (nstar), (3 * max_star + 1) * (3 * max_star + 1),
+	        TY_REAL)
+
+	    if (DP_NV (nstar) != NULL)
+	        call mfree (DP_NV (nstar), TY_REAL)
+	    call malloc (DP_NV (nstar), 3 * max_star + 1, TY_REAL)
+
+	} else {
+
+	    if (DP_NXCLAMP (nstar) != NULL)
+	        call mfree (DP_NXCLAMP (nstar), TY_REAL)
+	    call malloc (DP_NXCLAMP (nstar), max_star + 1, TY_REAL)
+
+	    if (DP_NXOLD (nstar) != NULL)
+	        call mfree (DP_NXOLD (nstar), TY_REAL)
+	    call malloc (DP_NXOLD (nstar), max_star + 1, TY_REAL)
+
+	    if (DP_NX (nstar) != NULL)
+	        call mfree (DP_NX (nstar), TY_REAL)
+	    call malloc (DP_NX (nstar), max_star + 1, TY_REAL)
+
+	    if (DP_NC (nstar) != NULL)
+	        call mfree (DP_NC (nstar), TY_REAL)
+	    call malloc (DP_NC (nstar), (max_star + 1) * (max_star + 1),
+	        TY_REAL)
+
+	    if (DP_NV (nstar) != NULL)
+	        call mfree (DP_NV (nstar), TY_REAL)
+	    call malloc (DP_NV (nstar), max_star + 1, TY_REAL)
+	}
+end
+
+
+# DP_NSCLOSE -- Procedure to close up the NSTAR parameters.
+
+procedure dp_nsclose (dp)
+
+pointer	dp		# pointer to daophot structure
+
+pointer	nstar
+
+begin
+	nstar = DP_NSTAR(dp)
+
+	if (DP_NNPIX (nstar) != NULL)
+	    call mfree (DP_NNPIX (nstar), TY_INT)
+	if (DP_NNUMER (nstar) != NULL)
+	    call mfree (DP_NNUMER (nstar), TY_REAL)
+	if (DP_NDENOM (nstar) != NULL)
+	    call mfree (DP_NDENOM (nstar), TY_REAL)
+	if (DP_NRPIXSQ (nstar) != NULL)
+	    call mfree (DP_NRPIXSQ (nstar), TY_REAL)
+	if (DP_NSKIP (nstar) != NULL)
+	    call mfree (DP_NSKIP (nstar), TY_INT)
+	if (DP_NXCLAMP (nstar) != NULL)
+	    call mfree (DP_NXCLAMP (nstar), TY_REAL)
+	if (DP_NXOLD (nstar) != NULL)
+	    call mfree (DP_NXOLD (nstar), TY_REAL)
+	if (DP_NX (nstar) != NULL)
+	    call mfree (DP_NX (nstar), TY_REAL)
+	if (DP_NV (nstar) != NULL)
+	    call mfree (DP_NV (nstar), TY_REAL)
+	if (DP_NSUMWT (nstar) != NULL)
+	    call mfree (DP_NSUMWT (nstar), TY_REAL)
+	if (DP_NC (nstar) != NULL)
+	    call mfree (DP_NC (nstar), TY_REAL)
+	if (DP_NIER (nstar) != NULL)
+	    call mfree (DP_NIER (nstar), TY_INT)
+
+	call mfree (nstar, TY_STRUCT)
+end
diff --git a/noao/digiphot/daophot/nstar/dpnconfirm.x b/noao/digiphot/daophot/nstar/dpnconfirm.x
new file mode 100644
index 00000000..ecbd3056
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpnconfirm.x
@@ -0,0 +1,34 @@
+include "../lib/daophotdef.h"
+
+# DP_NCONFIRM -- Procedure to confirm the critical nstar parameters.
+
+procedure dp_nconfirm (dao)
+
+pointer	dao		# pointer to the group structure
+
+int	dp_stati()
+
+begin
+	call printf ("\n")
+
+	# Confirm recentering and sky fitting.
+	call dp_vrecenter (dao)
+	call dp_vfitsky (dao)
+	if (dp_stati (dao, FITSKY) == NO)
+	    call dp_vgroupsky (dao)
+
+	# Confirm the psf radius.
+	call dp_vpsfrad (dao)
+
+	# Confirm the fitting radius.
+	call dp_vfitrad (dao)
+
+	# Confirm the maximum group size.
+	call dp_vmaxgroup (dao)
+
+	# Confirm the minimum and maximum good data values.
+	call dp_vdatamin (dao)
+	call dp_vdatamax (dao)
+
+	call printf ("\n")
+end
diff --git a/noao/digiphot/daophot/nstar/dpnstar.x b/noao/digiphot/daophot/nstar/dpnstar.x
new file mode 100644
index 00000000..66ba10a4
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpnstar.x
@@ -0,0 +1,355 @@
+include	<imhdr.h>
+include <tbset.h>
+include <mach.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/nstardef.h"
+
+
+# DP_NPHOT -- Read in the groups and fit all the stars in each group
+# simultaneously.
+
+procedure dp_nphot (dao, im, grp, nst, rej, ap_text)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# input image descriptor
+int	grp			# input group file descriptor
+int	nst			# output photometry file descriptor
+int	rej			# rejections file descriptor
+bool	ap_text			# photometry text file
+
+bool	converge 
+int	out_record, rout_record, in_record, nrow_in_table, old_size, niter
+int	cdimen, stat
+pointer	sp, indices, fields, icolpoint, ocolpoint, key, nstar, apsel
+real	radius, mean_sky
+
+int	dp_ggroup(), tbpsta(), dp_nstarfit(), dp_nxycheck()
+real	dp_nmsky()
+
+begin
+	# Store the original value of the fitting radius.
+	radius = DP_FITRAD(dao)
+	DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao))
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (icolpoint, NAPGROUP, TY_INT)
+	call salloc (indices, NAPPAR, TY_INT)
+	call salloc (fields, SZ_LINE, TY_CHAR)
+	call salloc (ocolpoint, NST_NOUTCOL, TY_INT)
+
+	# Get some daophot pointers.
+	apsel = DP_APSEL(dao)
+	nstar = DP_NSTAR (dao)
+
+	# Allocate space for and define the output table.
+	if (DP_TEXT(dao) == YES) { 
+	    call dp_xpnewnstar (dao, nst)
+	    if (rej != NULL)
+	        call dp_xpnewnstar (dao, rej)
+	} else {
+	    call dp_tpnewnstar (dao, nst, Memi[ocolpoint])
+	    if (rej != NULL)
+	        call dp_tpnewnstar (dao, rej, Memi[ocolpoint])
+	}
+
+	# Set up the input file.
+	call dp_gnindices (Memi[indices])
+	if (ap_text) {
+	    call pt_kyinit (key)
+	    call dp_gnstpsf (Memi[indices], Memc[fields], NAPGROUP)
+	    nrow_in_table = 0
+	} else {
+	    key = NULL
+	    call dp_tnstinit (grp, Memi[icolpoint])
+	    nrow_in_table = tbpsta (grp, TBL_NROWS)
+	}
+
+	# Allocate some memory for fitting.
+	call dp_memapsel (dao, Memi[indices], NAPPAR, DP_MAXGROUP(dao) + 1)
+	call dp_memnstar (dao, DP_MAXGROUP(dao) + 1)
+
+	# Initialize the input/output files for reading/writing.
+	in_record = 1
+	out_record = 0
+	rout_record = 0
+
+	repeat {
+
+	    # Read in the next group of stars.
+	    old_size = dp_ggroup (dao, grp, key, Memc[fields], Memi[indices],
+	        Memi[icolpoint], nrow_in_table, DP_MAXGROUP(dao), in_record,
+		DP_NGNUM(nstar))
+	    if (old_size <= 0)
+		break
+	    DP_NNUM(nstar) = old_size
+
+	    # Convert coordinates if necessary.
+	    call dp_win (dao, im, Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], DP_NNUM(nstar))
+
+	    # Print out the group number and number of stars.
+	    if (DP_VERBOSE (dao) == YES) {
+		call printf ("Group: %4d contains %2d stars\n")
+		    call pargi (DP_NGNUM (nstar))
+		    call pargi (DP_NNUM (nstar))
+	    }
+
+	    # If the group center is undefined or off image skip to next 
+	    # group.
+	    if (dp_nxycheck (im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+	        DP_NNUM(nstar), DP_FITRAD(dao), stat) <= 0) {
+
+		if (DP_VERBOSE(dao) == YES) {
+		    if (stat < 0)
+			call printf (
+			    "\tGroup %d center is undefined: skipping\n")
+		    else if (stat == 0)
+			call printf ("\tGroup %d is off the image: skipping\n")
+		    call pargi (DP_NGNUM(nstar))
+		}
+		next
+
+	    # If too many stars skip to the next group.
+	    } else if (DP_NNUM (nstar) > DP_MAXGROUP(dao)) {
+
+		if (DP_VERBOSE(dao) == YES) {
+		    call printf ("\tGroup %d larger than %d stars: skipping\n")
+			call pargi (DP_NGNUM(nstar))
+			call pargi (DP_MAXGROUP(dao))
+		}
+		niter = 0
+		call realloc (DP_NIER(nstar), DP_NNUM(nstar), TY_INT)
+		call dp_nstindef (dao, DP_NNUM(nstar), NSTERR_BIGGROUP)
+
+	    # If the mean sky value of group is undefined skip to next group.
+	    } else if (IS_INDEFR (dp_nmsky (Memr[DP_APMSKY(apsel)],
+	        DP_NNUM(nstar), mean_sky))) {
+
+		if (DP_VERBOSE(dao) == YES) {
+		    call printf (
+		        "\tGroup %d sky value is undefined: skipping\n")
+		        call pargi (DP_NGNUM(nstar))
+		}
+		niter = 0
+		call dp_nstindef (dao, DP_NNUM(nstar), NSTERR_INDEFSKY)
+
+	    # Do the fit.
+	    } else {
+
+		# Estimate values of INDEF magnitudes.
+		call dp_nmaginit (dao, Memr[DP_APMAG(apsel)],
+		    Memr[DP_APERR(apsel)], DP_NNUM(nstar))
+
+		# Set up for the fit.
+		if (DP_RECENTER(dao) == YES)
+	            cdimen = 3 * DP_NNUM(nstar) + 1
+		else
+	            cdimen = DP_NNUM(nstar) + 1
+
+	        # Iterate.
+		for (niter = 1; niter <= DP_MAXITER(dao); niter = niter + 1) {
+	            DP_NNUM(nstar) = dp_nstarfit (dao, im, DP_NNUM(nstar),
+		        mean_sky, cdimen, niter, converge)
+		    if (DP_NNUM(nstar) <= 0)
+			break
+		    if (converge)
+			break
+	        }
+
+		# Solution did not converge.
+		niter = min (niter, DP_MAXITER(dao))
+	    }
+
+	    # Now write out the results.
+	    if (DP_TEXT(dao) == YES)
+	        call dp_xntwrite (dao, im, nst, rej, niter, old_size)
+	    else
+	        call dp_tntwrite (dao, im, nst, rej, niter, old_size,
+		    out_record, rout_record, Memi[ocolpoint])
+	} 
+
+	if (ap_text)
+	    call pt_kyfree (key)
+
+	call sfree (sp)
+
+	# Restore the original value of the fitting radius.
+	DP_FITRAD(dao) = radius
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+end
+
+
+# DP_NSTINDEF -- Set magnitudes and fitting parameters of unfit stars to indef.
+
+procedure dp_nstindef (dao, max_nstars, pier)
+
+pointer	dao		# pointer to the daophot structure
+int	max_nstars	# number of stars
+int	pier		# the photometry error code
+
+int	i
+pointer	apsel, nstar
+
+begin
+	apsel = DP_APSEL(dao)
+	nstar = DP_NSTAR(dao)
+	do i = 1, max_nstars {
+	    Memr[DP_APMAG(apsel)+i-1] = INDEFR
+	    Memr[DP_APERR(apsel)+i-1] = INDEFR
+	    Memr[DP_APCHI(apsel)+i-1] = INDEFR
+	    Memr[DP_APSHARP(apsel)+i-1] = INDEFR
+	    Memi[DP_NIER(nstar)+i-1] = pier
+	}
+end
+
+
+# DP_GNINDICES -- Get the memory allocation fields.
+
+procedure dp_gnindices (indices)
+
+int	indices[ARB]		# index array
+
+begin
+	indices[1] = DP_PAPID
+        indices[2] = DP_PAPXCEN
+	indices[3] = DP_PAPYCEN
+	indices[4] = DP_PAPMAG1
+	indices[5] = DP_PAPSKY
+	indices[6] = DP_PAPGROUP
+	indices[7] = DP_PAPMERR1
+	indices[8] = DP_PAPNITER
+	indices[9] = DP_PAPSHARP
+	indices[10] = DP_PAPCHI
+end
+
+
+define	GRP_CTRINDEF	-1
+define	GRP_CTROFFIMAGE	 0
+define	GRP_CTROK	 1
+
+# DP_NXYCHECK -- Check that the center of the group is defined. -1 is
+# returned if the center of the group is undefined, 0, is returned if the
+# group is entirely off the image, 1 is returned if the group is at least
+# partially on the input image.
+
+int procedure dp_nxycheck (im, x, y, group_size, radius, stat)
+
+pointer	im			# pointer to the input image
+real	x[ARB]			# array of x values
+real	y[ARB]			# array of y values
+int	group_size		# the size of the group
+real	radius			# the fitting radius
+int	stat			# the return status
+
+int	i, nxy
+real	xmin, xmax, ymin, ymax, xx, yy
+
+begin
+	# Initialize.
+	stat = GRP_CTRINDEF
+	xmin = MAX_REAL
+	xmax = -MAX_REAL
+	ymin = MAX_REAL
+	ymax = -MAX_REAL
+
+	# Compute the minimum and maximum x and y values.
+	nxy = 0
+	do i = 1, group_size {
+	    xx = x[i]
+	    yy = y[i]
+	    if (IS_INDEFR(xx) || IS_INDEFR(yy))
+		next
+	    nxy = nxy + 1
+	    if (xx < xmin)
+		xmin = xx
+	    if (xx > xmax)
+		xmax = xx
+	    if (yy < ymin)
+		ymin = yy
+	    if (yy > ymax)
+		ymax = yy
+	}
+	if (nxy <= 0)
+	    return (stat)
+
+	# Test the min and max values. 
+	stat = GRP_CTROFFIMAGE
+	if ((int (xmin - radius) + 1) > IM_LEN(im,1))
+	    return (stat)
+	if (int (xmax + radius) < 1)
+	    return (stat)
+	if ((int (ymin - radius) + 1) > IM_LEN(im,2))
+	    return (stat)
+	if (int (ymax + radius) < 1)
+	    return (stat)
+	
+	# The group is on the image.
+	stat = GRP_CTROK
+	return (stat)
+end
+
+
+# DP_NMSKY -- Compute the mean sky value for the group of stars.
+
+real procedure dp_nmsky (sky, group_size, msky)
+
+real	sky[ARB]		# the array of sky values
+int	group_size		# the size of the group of stars
+real	msky			# the mean sky value
+
+int	i, nsky
+real	sky_sum
+
+begin
+	sky_sum = 0.0
+	nsky = 0
+
+	do i = 1, group_size {
+	    if (IS_INDEFR(sky[i]))
+		next
+	    sky_sum = sky_sum + sky[i]
+	    nsky = nsky + 1
+	}
+
+	if (nsky <= 0)
+	    msky = INDEFR
+	else
+	    msky = sky_sum / nsky
+
+	return (msky)
+end
+
+
+define	MIN_REL_BRIGHT		1.0E-04		# minimum relative brightness
+define	INIT_REL_BRIGHT		0.01		# initial relative brightness
+
+# DP_NMAGINIT --  Initialize the magnitude and magnitude error arrays before
+# fitting the group.
+
+procedure dp_nmaginit (dao, mag, magerr, group_size)
+
+pointer	dao			# pointer to the daophot strucuture
+real	mag[ARB]		# the magnitude array
+real	magerr[ARB]		# the magnitude error array
+int	group_size		# size of the group
+
+int	i
+pointer	psffit
+
+begin
+	psffit = DP_PSFFIT (dao)
+	do i = 1, group_size {
+	    if (IS_INDEFR(mag[i]))
+		mag[i] = INIT_REL_BRIGHT
+	    else {
+	        mag[i] = DAO_RELBRIGHT (psffit, mag[i])
+	        if (mag[i] <= MIN_REL_BRIGHT) 
+		    mag[i] = INIT_REL_BRIGHT
+	    }
+	    magerr[i] = 0.0
+	}
+end
diff --git a/noao/digiphot/daophot/nstar/dpnstarfit.x b/noao/digiphot/daophot/nstar/dpnstarfit.x
new file mode 100644
index 00000000..8b864c86
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpnstarfit.x
@@ -0,0 +1,1383 @@
+include <mach.h>
+include	<imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/nstardef.h"
+
+
+# DP_NSTARFIT -- Fit the stellar group.
+
+int procedure dp_nstarfit (dao, im, nin_group, mean_sky, cdimen, iter,
+	converge)
+
+pointer	dao			# pointer to the daophot structure
+pointer im			# pointer to the input image
+int	nin_group		# original group size
+real	mean_sky		# the mean sky for the group
+int	cdimen			# dimensions of the coefficient matrix
+int	iter			# the current iteration
+bool	converge		# did the fit converge ?
+
+bool 	clip, refit
+int	i, j, k, xpix, ypix, group_size, nterm, ncols, mindex, ifaint, tifaint
+int	ntmin, ixmin, ixmax, iymin, iymax, ifxmin, ifxmax, ifymin, ifymax, flag
+pointer	apsel, psffit, nstar, subim, ypixel, pixel
+real	fitradsq, cutoff, psfradsq, sepcrit, sepmin, perr, peakerr, sky_value
+real	read_noise, mingdata, maxgdata, chigrp, wcrit, xmin, xmax, ymin, ymax
+real	datum, sumres, grpwt, xtemp, ytemp, weight, ds, pred_pixval, sigmasq
+real	relerr, dswt, faint, tfaint, noise
+
+bool	dp_nstmerge(), dp_ntomit(), dp_ntmin(), dp_ncheckc()
+pointer imgs2r()
+real	dp_ntskyval(), dp_ntsubtract()
+
+begin
+	# Define the daophot pointers.
+	psffit = DP_PSFFIT (dao)
+	apsel = DP_APSEL(dao)
+	nstar = DP_NSTAR (dao)
+
+	# Set up some daophot constants. At some point these will be computed
+	# when the NSTAR task is started up instead of at the beginning of
+	# each group fit. For the moment it is convenient and not too
+	# costly to compute them here. Also initialize the fit.
+
+	if (iter == 1) {
+
+	    # Compute the fitting and psf radii.
+	    fitradsq = DP_FITRAD (dao) ** 2
+	    psfradsq = DP_PSFRAD(dao) ** 2
+	    cutoff = CUT_FACTOR * fitradsq
+
+	    # Compute the merging parameters.
+	    if (IS_INDEFR(DP_MERGERAD(dao))) {
+	        sepcrit = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 +
+	            Memr[DP_PSFPARS(psffit)+1] ** 2)
+	        sepmin = min (1.0, FRACTION_MINSEP * sepcrit)
+	    } else {
+		sepcrit = DP_MERGERAD(dao) ** 2
+	        sepmin = min (1.0, FRACTION_MINSEP * sepcrit)
+	    }
+
+	    # Compute the noise parameters.
+	    read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+	    perr = 0.01 * DP_FLATERR(dao)
+	    peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] *
+	        Memr[DP_PSFPARS(psffit)+1])
+
+	    # Compute the minimum and maximum good pixel values.
+	    if (IS_INDEFR(DP_MINGDATA(dao)))
+	        mingdata = -MAX_REAL
+	    else
+	        mingdata = DP_MINGDATA(dao)
+	    if (IS_INDEFR(DP_MAXGDATA(dao)))
+	        maxgdata = MAX_REAL
+	    else
+	        maxgdata = DP_MAXGDATA(dao)
+
+	    # Initialize the fit.
+	    chigrp = 1.0
+	    clip = false
+	    if (DP_RECENTER(dao) == YES) {
+	        nterm = 3 * nin_group
+		ntmin = 3
+	    } else {
+	        nterm = nin_group
+		ntmin = 1
+	    }
+	    if (DP_FITSKY(dao) == YES) {
+		nterm = nterm + 1
+		ntmin = ntmin + 1
+	    }
+	    call aclrr (Memr[DP_NXOLD(nstar)], nterm) 
+	    call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+	}
+
+	# Start a new iteration.
+	group_size = nin_group
+	if (DP_RECENTER(dao) == YES)
+	    nterm = 3 * group_size
+	else
+	    nterm = group_size
+	if (DP_FITSKY(dao) == YES)
+	    nterm = nterm + 1
+
+	# Begin fitting the current group of stars.
+	repeat {
+
+	    # Initialize the convergence criteria.
+	    converge = false
+
+	    # Check that there is at least 1 star in the group.
+	    if (group_size < 1)
+	        return (group_size)
+
+	    # Set up the critical error for star rejection.
+	    if (iter >= NITER_MAX)
+	        wcrit = WCRIT_MAX
+	    else if (iter >= NITER_MED)
+	        wcrit = WCRIT_MED
+	    else if (iter >= NITER_MIN)
+	        wcrit = WCRIT_MIN
+	    else
+	        wcrit = MAX_REAL
+
+	    # Set the sky fitting derivative.
+	    if (DP_FITSKY(dao) == YES)
+	        Memr[DP_NX(nstar)+nterm-1] = -1.0
+
+	    # Initialize arrays.
+	    call aclrr (Memr[DP_APCHI(apsel)], group_size)
+	    call aclrr (Memr[DP_NSUMWT(nstar)], group_size)
+	    call aclrr (Memr[DP_NNUMER(nstar)], group_size)
+	    call aclrr (Memr[DP_NDENOM(nstar)], group_size)
+	    call amovki (NSTERR_OK, Memi[DP_NIER(nstar)], group_size)
+
+	    # Compute the minimum and maximum x and y values.
+	    call alimr (Memr[DP_APXCEN(apsel)], group_size, xmin, xmax) 
+	    call alimr (Memr[DP_APYCEN(apsel)], group_size, ymin, ymax) 
+
+	    # Check to see whether any two stars are within the critical
+	    # difference from each other.
+
+	    if ((group_size > 1) && dp_nstmerge (Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		Memr[DP_APERR(apsel)], group_size, sepcrit, sepmin, wcrit,
+		i, j, k)) {
+
+		# Compute the new centroid and brightness.
+		call dp_nstcen (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], i, j, k)
+
+		# Print out verbose comments.
+		if (DP_VERBOSE (dao) == YES)  {
+		    call printf ("\tStar %-5d has merged with star %-5d\n")
+		        call pargi (Memi[DP_APID(apsel)+k-1])
+		        call pargi (Memi[DP_APID(apsel)+i-1])
+		}
+
+		# Recompute the mean sky.
+		if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES))
+		    mean_sky = (mean_sky * group_size -
+		        Memr[DP_APMSKY(apsel)+k-1]) / (group_size - 1)
+
+		# Now remove the k-th star from the group.
+		call dp_remove (k, group_size, Memi[DP_APID(apsel)],
+		    Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)],
+		    Memi[DP_NIER(nstar)], NSTERR_MERGE)
+
+		# After deleting a star, resize the matrix, release all of
+		# the clamps and back up the iteration counter.
+
+		if (DP_RECENTER(dao) == YES)
+		    nterm = 3 * group_size
+		else
+		    nterm = group_size
+		if (DP_FITSKY(dao) == YES)
+		    nterm = nterm + 1
+		clip = false
+		call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+		call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+		iter = max (1, iter - 1)
+		next
+	    }
+
+	    # Now we can proceed with the iteration. If this is the first
+	    # iteration read in the subraster. Determine the size of the
+	    # subraster we need to extract from the image for this group.
+
+	    if (iter == 1) {
+	        ixmin = max (1, int (xmin - DP_PSFRAD(dao) -
+		    DP_FITRAD(dao)) + 1)
+	        iymin = max (1, int (ymin - DP_PSFRAD(dao) -
+		    DP_FITRAD(dao)) + 1)
+	        ixmax = min (IM_LEN(im,1), int (xmax + DP_PSFRAD(dao) +
+		    DP_FITRAD(dao)))
+	        iymax = min (IM_LEN(im,2), int (ymax + DP_PSFRAD(dao) +
+		    DP_FITRAD(dao)))
+	        subim = imgs2r (im, ixmin, ixmax, iymin, iymax)
+	        ncols = ixmax - ixmin + 1
+	    }
+
+	    # Compute the area on the subraster that is off interest to
+	    # the current iteration.
+
+	    ifxmin = max (ixmin, int (xmin - DP_FITRAD(dao)) + 1)
+	    ifymin = max (iymin, int (ymin - DP_FITRAD(dao)) + 1)
+	    ifxmax = min (ixmax, int (xmax + DP_FITRAD(dao)))
+	    ifymax = min (iymax, int (ymax + DP_FITRAD(dao)))
+
+	    # Zero the normal matrix and the vector of residuals.
+
+	    call aclrr (Memr[DP_NV(nstar)], nterm)
+	    call aclrr (Memr[DP_NC(nstar)], cdimen * cdimen)
+	    call aclri (Memi[DP_NNPIX(nstar)], group_size)
+
+	    sumres = 0.0
+	    grpwt = 0.0
+
+	    # Loop over the pixels in the subraster.
+
+	    ypixel = subim + (ifymin - iymin) * ncols + ifxmin - ixmin - 1
+	    do ypix = ifymin, ifymax {
+
+		ytemp = ypix
+	        pixel = ypixel
+
+	        do xpix = ifxmin, ifxmax {
+
+		    xtemp = xpix
+		    pixel = pixel + 1
+		    datum = Memr[pixel]
+
+		    # Reject pixel if not in good data range.
+		    if ((datum < mingdata) || (datum > maxgdata))
+		        next
+
+		    # If this pixel is within one fitting radius of at
+		    # least one star then include it in the solution.
+		    # While figuring this out, compute the squared distance
+		    # of this pixel from the centroid of each star in the
+		    # group, and sum its contribution to the number
+		    # of pixels within one fitting radius of each star.
+
+		    if (dp_ntomit (Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)],
+			Memi[DP_NSKIP(nstar)], group_size, xtemp, ytemp,
+			cutoff))
+		        next
+
+		    if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == NO)) {
+			sky_value = dp_ntskyval (Memr[DP_APMSKY(apsel)],
+			    Memi[DP_NSKIP(nstar)], group_size)
+			if (IS_INDEFR(sky_value))
+			    sky_value = mean_sky
+		    } else
+			sky_value = mean_sky
+
+		    # Subtract the mean sky from the pixel.
+		    #ds = datum - mean_sky
+		    ds = datum - sky_value
+
+		    # Now loop over the stars and subtract from this pixel
+		    # the light contribution from each star within one psf
+		    # radius.
+		    #
+		    # The condition equation for pixel[xpix,ypix] is
+		    #
+		    # residual = data[xpix,ypix] - sky - sum (scale *
+		    #            psf[xpix-xc,ypix-yc])
+		    #
+		    # The scale, xc, and yc are iterated until
+		    #
+		    #          sum (weight * residual ** 2)
+		    #
+		    # is minimized. Weight will be a function of 1) the
+		    # distance of the pixel from the center of the nearest
+		    # star, 2) the model-predicted brightness of the pixel
+		    # (taking into consideration the readout noise, the
+		    # photons/ADU, and the interpolating error of the PSF,
+		    # and, 3) the size of the residual itself. One is
+		    # necessary to prevent the non-linear least squares
+		    # solution from oscillating: oftimes it will come
+		    # to pass that if you include a pixel in the solution
+		    # then the predicted shift of the centroid will cause
+		    # that pixel to be excluded in the next iteration, and the 
+		    # new predicted shift of the centroid will cause that
+		    # pixel to be included again. This could go on ad
+		    # infinitum. The cure is to have the weight of the pixel
+		    # go asymptotically to zero as its distance from the
+		    # stellar centroid approaches the fitting radius. In
+		    # the case like that just described, the solution can
+		    # then find a real minimum of the sum of the weighted
+		    # squared residuals with that pixel at some low-weight
+		    # position just inside the fitting radius. Two is
+		    # just sensible weighting. Three is a crude attempt
+		    # at making the solution more robust against bad pixels.
+
+		    weight = dp_ntsubtract (dao, im, Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)],
+		        Memr[DP_APMAG(apsel)], Memi[DP_NSKIP(nstar)],
+		        Memr[DP_NX(nstar)], group_size, xtemp, ytemp, ds,
+			psfradsq, fitradsq, DP_RECENTER(dao))
+
+		    # At this point the vector X contains the first
+		    # derivative of the condition equation, for the pixel
+		    # under consideration, with respect to each of the
+		    # fitting parameters for all of these stars. We now need
+		    # to add these values into the normal matrix and the vector
+		    # of residuals.
+
+		    # The expected random error in the pixel is the quadratic
+		    # sum of the Poisson statistics, plus the readout noise,
+		    # plus an estimated error of 0.75% of the total brightness
+		    # of the total brightness for the difficulty of flat-
+		    # fielding and bias subtraction, plus an estimated error
+		    # of the same fraction of the fourth derivative at the
+		    # peak of the profile, to account for the difficulty 
+		    # of accurately interpolating within the point-spread
+		    # function. The fourth derivative of the PSF is
+		    # proportional to H / sigma ** 4 (sigma is the Gaussian
+		    # width parameter for the stellar core); using the geometric
+		    # mean of sigma(x) and sigma(y), this becomes H / [sigma(x)
+		    # * sigma(y)] ** 2. The ratio of the fitting error to this
+		    # quantity is estimated to be approximately 0.027 from a
+		    # good-seeing CTIO frame. (This is probably a function of
+		    # seeing, sampling etc.)
+
+		    # Get the residual from the PSF fit and the pixel
+		    # intensity as predicted by the fit. Pred_pixval = raw data
+		    # minus residual = model predicted value of the intensity at
+		    # this point.
+
+		    pred_pixval = max (0.0, datum - ds)
+		    if ((pred_pixval > maxgdata) && (iter >= 4))
+			next
+
+		    #sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + 
+		        #(perr * pred_pixval) ** 2 + (peakerr *
+		        #(pred_pixval - mean_sky)) ** 2
+		    sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + 
+		        (perr * pred_pixval) ** 2 + (peakerr *
+		        (pred_pixval - sky_value)) ** 2
+		    if (sigmasq <= 0.0)
+			next
+		    relerr = abs (ds) / sqrt (sigmasq)
+
+		    # Add this residual into the weighted sum of the
+		    # absolute relative residuals.
+
+		    sumres = sumres + relerr * weight
+		    grpwt = grpwt + weight
+
+		    # Add into the accumulating sums of the weighted 
+		    # absolute relative residuals and of the image sharpness
+		    # parameter for each of the stars. Include in the
+		    # sharpness index only those pixels within NCORE_SIGMA
+		    # sigma of the centroid of the object. This saves time
+		    # and floating underflows by excluding pixels
+		    # which contribute less than about one part in a
+		    # million to the index.
+
+		    call dp_acsharp (Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memi[DP_NSKIP(nstar)],
+			Memi[DP_NNPIX(nstar)], Memr[DP_NNUMER(nstar)],
+			Memr[DP_NDENOM(nstar)], Memr[DP_NSUMWT(nstar)],
+			Memr[DP_APCHI(apsel)], group_size, xtemp, ytemp,
+			Memr[DP_PSFPARS(psffit)], Memr[DP_PSFPARS(psffit)+1],
+			ds, sigmasq, relerr, weight)
+
+		    # If clipping is in effect, reduce the weight of a bad
+		    # pixel. A pixel having a residual of 2.5 sigma gets
+		    # reduced to half weight and one with a rersidual of 5
+		    # sigma gets weight of 1 / 257.
+
+		    weight = weight / sigmasq
+		    if (clip) 
+		        weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) *
+			    chigrp)) ** DP_CLIPEXP(dao))
+		    dswt = ds * weight
+
+		    # Work this pixel into the normal matrix.
+		    call dp_mataccum (Memr[DP_NX(nstar)], Memi[DP_NSKIP(nstar)],
+		        group_size, Memr[DP_NC(nstar)], Memr[DP_NV(nstar)],
+			cdimen, nterm, weight, dswt, DP_RECENTER(dao),
+			DP_FITSKY(dao))
+
+	        }
+
+	        ypixel = ypixel + ncols
+	    }
+
+	    # Make sure that every star in the group has at least MIN_FIT_PIXEL
+	    # pixels within one fitting radius.
+
+	    refit = dp_ntmin (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APMSKY(apsel)], Memi[DP_NNPIX(nstar)],
+		Memi[DP_NIER(nstar)], group_size, nterm, DP_RECENTER(dao),
+		DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, DP_VERBOSE(dao))
+	    if (group_size < 1)
+	        return (group_size)
+	    if (refit)
+	        next
+
+	    # Reflect the normal matrix across the diagonal.
+	    call dp_mreflect (Memr[DP_NC(nstar)], cdimen, nterm)
+
+	    # Compute the robust estimate of the standard deviation of the
+	    # residuals for the group as a whole, and for each star. This
+	    # estimate is sqrt (PI/2) * weighted mean absolute relative
+	    # residual. Do you like that "absolute relative" stuff? (PBS)
+	    # NO! (LED)
+	    #
+	    #      CHI = CHI_NORM * SUM (weight * resid) / (# of pixels)
+	    # 
+	    # This gets corrected for bias by being multiplied by
+	    #
+	    #	SQRT (# of pixels) / (# of pixels - 3)
+
+	    # Then the value is driven towards unity, depending on
+	    # exactly how many pixels were involved: if CHIOLD is based
+	    # on a total weight of 3, then it is extremely poorly determined
+	    # and we just want to keep CHIOLD = 1. The larger GRPWT is, the
+	    # better determined CHIOLD is, and the less we want to force it
+	    # toward unity. So, just take the weighted average of CHIOLD and 
+	    # unity, with weights GRPWT - 3 and 1, respectively.
+
+	    if (grpwt > ntmin) {
+	        chigrp = CHI_NORM * sumres * sqrt (1.0 / (grpwt * (grpwt -
+	            ntmin)))
+	        chigrp = ((grpwt - ntmin) * chigrp + ntmin) / grpwt
+	    }
+
+	    # CHIOLD has been pulled toward its expected value of unity to
+	    # keep the statistics of a small number of pixels from completely
+	    # dominating the error analysis. Similarly, the photometric 
+	    # errors for the individual stars will be pulled toward unity
+	    # now. Later on, if the number of stars in the group is 
+	    # greated than one, CHI will be nudged toward the group average.
+	    # In order to work optimally, of course, this requires that	
+	    # the # of photons / ADU, the READNOISE and the other noise
+	    # contributors are properly specified.
+
+	    call dp_ntchi (Memr[DP_NSUMWT(nstar)], Memr[DP_APCHI(apsel)],
+	        group_size, ntmin, chigrp, grpwt)
+
+	    # Invert the matrix.
+	    call invers (Memr[DP_NC(nstar)], cdimen, nterm, flag)
+
+	    # Check for a singular matrix.
+	    refit = dp_ncheckc (Memr[DP_NC(nstar)], cdimen, nterm,
+		Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size,
+		DP_RECENTER(dao), DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky,
+		DP_VERBOSE(dao)) 
+	    if (group_size < 1)
+	        return (group_size)
+	    if (refit)
+	        next
+
+	    # Solve for position and scale factor increments.
+	    call mvmul (Memr[DP_NC(nstar)], cdimen, nterm, Memr[DP_NV(nstar)],
+	        Memr[DP_NX(nstar)])
+
+	    if (iter <= 1)
+	        refit = true
+	    else
+	        refit = false
+
+	    # Fit the sky.
+	    if (DP_FITSKY(dao) == YES) {
+	        noise = sqrt (abs (mean_sky / DP_PHOTADU(dao)) + read_noise)
+	        mean_sky = mean_sky - max (-3.0 * noise,
+		    min (Memr[DP_NX(nstar)+nterm-1], 3.0 * noise))
+		if (abs (Memr[DP_NX(nstar)+nterm-1]) > (1.0e-4 * mean_sky))
+		    refit = true
+	    }
+
+	    # In the beginning, the brightness of each star will be permitted
+	    # to change by no more than 2 magnitudes per iteration, and the x,y
+	    # coordinates of each centroid will be permitted to change by
+	    # no more than 0.4 pixels per iteration. Any time that the 
+	    # parameter correction changes sign from one iteration to the
+	    # next, the maximum permissible change will be reduced by a factor
+	    # of two. These clamps are released any time a star disappears.
+
+	    call dp_ntclamp (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+	        Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+	        Memr[DP_NSUMWT(nstar)], group_size, Memr[DP_NC(nstar)], cdimen,
+		Memr[DP_NX(nstar)], Memr[DP_NXOLD(nstar)],
+		Memr[DP_NXCLAMP(nstar)], DP_RECENTER(dao), clip, refit)
+
+	    # Check whether the estimated centroid of the any star has
+	    # moved so far out of the limits of the picture that it has fewer
+	    # than 4 or 5 pixels within one fitting radius.
+
+	    call dp_ntcentroid (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size,
+		ixmin, ixmax, iymin, iymax, fitradsq, DP_FITSKY(dao),
+		DP_GROUPSKY(dao), mean_sky, refit, DP_VERBOSE(dao))
+
+	    if (group_size < 1)
+	        return (group_size)
+
+	    # Update matrix dimensions.
+	    if (DP_RECENTER(dao) == YES)
+	        nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (DP_FITSKY(dao) == YES)
+		nterm = nterm + 1
+
+	    # Now check whether any of the stars is too faint (more than 12.5
+	    # magnitudes fainter than the PSF star). If several stars are too
+	    # faint, delete the faintest one, and set the brightness of the 
+	    # other faint ones to 12.5 magnitudes below the PSF star. That way
+	    # on the next iteration we will see whether these stars want to 
+	    # grow or to disappear.
+
+	    faint = 0.0
+	    ifaint = 0
+	    call dp_ntfmag (Memr[DP_APMAG(apsel)], group_size, tfaint, tifaint)
+
+	    # If at least one star is more than 12.5 magnitudes fainter
+	    # than the PSF then ifaint is the relative index of the faintest
+	    # of them, and faint is the relative brightness of the
+	    # faintest of them.
+
+	    # No very faint star was detected.
+	    if (tifaint <= 0) {
+
+	        # If the solution has not converged and if the number of
+	        # iterations is less than MIN_ITER, perform another iteration
+	        # with no questions asked.
+
+	        if ((refit) && (iter < MIN_ITER))
+		    return (group_size)
+	    
+	        # If the solution doesn't think it has converged, after the
+	        # fourth iteration delete the least certain star if it is less
+	        # less than a one-sigma detection; after the eighth iteration
+	        # delete the least certain star if it is less than a 1.5 sigma
+	        # detection; after the twelfth iteration OR if the solution
+	        # thinks it has converged, delete the least certain star if it 
+	        # is less than a two-sigma detection.
+
+	        call dp_fsnoise (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		    group_size, faint, ifaint)
+
+	        if ((refit) && (iter < DP_MAXITER (dao)) && (faint < wcrit))
+		    return (group_size)
+	    }
+
+	    # Reject the appropriate star.
+	    if ((tifaint > 0) || (faint >= MIN_FAINT)) {
+
+	        # Either (a) the solution thinks it has not converged
+	        # and the faintest star is more uncertain than sqrt(wcrit)
+	        # or (b) the solution thinks it has converged and the
+	        # faintest star is more uncertain than two-sigma.
+
+	        if (DP_VERBOSE (dao) == YES) {
+		    mindex = max (tifaint, ifaint)
+		    call printf (
+		    "\tStar %-5d has been deleted because it is too faint\n")
+		        call pargi (Memi[DP_APID(apsel)+mindex-1])
+	        }
+
+		if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES) &&
+		    (group_size > 1))
+		    mean_sky = (mean_sky * group_size -
+		        Memr[DP_APMSKY(apsel)+max(tifaint,ifaint)-1]) /
+			    (group_size - 1)
+
+	        call dp_remove (max (tifaint, ifaint), group_size,
+	            Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+		    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		    Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)],
+		    NSTERR_FAINT)
+
+	        if (group_size < 1)
+		    return (group_size)
+
+		if (DP_RECENTER(dao) == YES)
+	            nterm = 3 * group_size
+		else
+		    nterm = group_size
+		if (DP_FITSKY(dao) == YES)
+		    nterm = nterm + 1
+	        call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+	        call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+	        clip = false
+	        iter = max (1, iter - 1)
+	        next
+	    }
+
+	    # Solution has either converged or gone to the maximum number
+	    # of iterations.
+
+	    if ((iter < DP_MAXITER(dao)) && (! clip)) {
+
+	        # The first convergence milestone has been reached. Turn on the
+	        # clipper, loosen the clamps and keep on going.
+
+		if (DP_CLIPEXP(dao) > 0)
+	            clip = true
+	        converge = false
+	        call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+	        call amaxkr (Memr[DP_NXCLAMP(nstar)], 0.25,
+	            Memr[DP_NXCLAMP(nstar)], nterm)
+	        return (group_size)
+	    }
+
+	    converge = true		    
+
+	} until (converge)
+
+	return (group_size)
+end
+
+
+# DP_NSTMERGE -- Decide whether two stars in a group should merge.
+
+bool procedure dp_nstmerge (xcen, ycen, mag, magerr, group_size, sepcrit,
+	sepmin, wcrit, i, j, k)
+
+real	xcen[ARB]		# array of x centers	
+real	ycen[ARB]		# array of y centers	
+real	mag[ARB]		# array of magnitudes
+real	magerr[ARB]		# array of magnitude errors
+int	group_size		# group size
+real	sepcrit			# the critical separation squared
+real	sepmin			# the minimum separation
+real	wcrit			# critical error for rejection
+int	i, j, k			# output indices
+
+real	separation
+
+begin
+	do i = 1, group_size {
+  	    do j = 1,  i - 1 {
+
+		# Compute the separation.
+	        separation = (xcen[j] - xcen[i]) ** 2 +
+		    (ycen[j] - ycen[i]) ** 2
+		if (separation > sepcrit)
+		    next
+
+		# Find the fainter of the two stars.
+		k = j
+		if (mag[i] < mag[j]) 
+		    k = i
+		if ((separation <  sepmin) || ((magerr[k] / mag[k])  > wcrit))
+		        return (true)
+	    }
+	}
+
+	return (false)
+end
+
+
+# DP_NSTCEN --  Recompute the centroid and brightness of the i-th star.
+
+procedure dp_nstcen (xcen, ycen, mag, i, j, k)
+
+real	xcen[ARB]		# the x centers
+real	ycen[ARB]		# the y centers
+real	mag[ARB]		# the magnitudes
+int	i, j, k			# array indices
+
+begin
+	# Now eliminate the fainter of the two stars. The k-th
+	# star is now the fainter of the two, the i-th the
+	# brighter.
+
+        if (mag[i] < mag[j])
+            i = j
+
+	# Now replace the centroid of the i-th star with the
+	# weighted mean of the most recent estimates of the 
+	# centroids of the i-th and the k-th stars, and the
+	# brightness of i-th with the sum of the brightnesses.
+
+	xcen[i] = xcen[i] * mag[i] + xcen[k] * mag[k]
+	ycen[i] = ycen[i] * mag[i] + ycen[k] * mag[k]
+	mag[i] = mag[i] + mag[k]
+	xcen[i] = xcen[i] / mag[i]
+	ycen[i] = ycen[i] / mag[i]
+end
+
+
+# DP_NTOMIT -- Check whether a pixel is within one fitting radius of another
+# star. 
+
+bool procedure dp_ntomit (xcen, ycen, rpixsq, skip, group_size, fx, fy, cutoff)
+
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	rpixsq[ARB]		# array of distances squared
+int	skip[ARB]		# array of skip values
+int	group_size		# the group size
+real	fx, fy			# pixel position in image
+real	cutoff			# radius cuttoff for pixel inclusion
+
+bool	omit
+int	i
+
+begin
+	omit = true
+	do i = 1, group_size {
+	    skip[i] = YES
+	    rpixsq[i] = (fx - xcen[i]) ** 2 + (fy - ycen[i]) ** 2
+	    if (rpixsq[i] > cutoff)
+		next
+	    skip[i] = NO
+	    omit = false
+	}	 
+
+	return (omit)
+end
+
+
+# DP_NTSKYVAL -- Compute the average sky value to use for a particular
+# pixel by averaging the sky values of all stars for which the
+# pixel is within one fitting radius.
+
+real procedure dp_ntskyval (sky, skip, nstar)
+
+real	sky[ARB]		# array of sky values
+int	skip[ARB]		# array of skip values
+int	nstar			# the number of stars
+
+int	i, npts
+real	sum
+
+begin
+	sum = 0.0
+	npts = 0
+	do i = 1, nstar {
+	    if (skip[i] == YES)
+		next
+	    sum = sum + sky[i]
+	    npts = npts + 1
+	}
+	if (npts <= 0)
+	    return (INDEFR)
+	else
+	    return (sum / npts)
+end
+
+
+# DP_NTSUBTRACT -- Procedure to subtract the contribution of a particular
+# pixel from a particular star.
+
+real procedure dp_ntsubtract (dao, im, xcen, ycen, rpixsq, mag, skip, x,
+	group_size, fx, fy, ds, psfradsq, fitradsq, recenter)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	rpixsq[ARB]		# array of distances squared
+real	mag[ARB]		# array of magnitudes
+int	skip[ARB]		# array of skip values
+real	x[ARB]			# x accumulation vector array
+int	group_size		# size of the group
+real	fx, fy			# center of pixel in image
+real	ds			# pixel value
+real	psfradsq		# psf radius squared
+real	fitradsq		# fit radius squared
+int	recenter		# recenter the coordinates
+
+int	i, i3, k
+pointer	psffit
+real	weight, dx, dy, deltax, deltay, val, dvdx, dvdy, rsq
+real	dp_usepsf()
+
+begin
+	psffit = DP_PSFFIT(dao)
+
+	weight = 0.0
+	do i = 1, group_size {
+	    if (rpixsq[i] >= psfradsq)
+		next
+	    dx = fx - xcen[i]
+	    dy = fy - ycen[i]
+	    call dp_wpsf (dao, im, xcen[i], ycen[i], deltax, deltay, 1)
+	    deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+	    deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+	    val = dp_usepsf (DP_PSFUNCTION(psffit), dx, dy,
+	        DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+		Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit),
+		DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), deltax, deltay,
+		dvdx, dvdy)
+	    ds = ds - mag[i] * val
+	    if (skip[i] == YES)
+		next
+	    rsq = rpixsq[i] / fitradsq
+	    weight = max (weight, 5.0 / (5.0 + rsq / (1.0 - rsq)))
+	    if (recenter == YES) {
+	        i3 = 3 * i
+	        k = i3 - 2
+	        x[k] = -val
+	        k = i3 - 1
+	        x[k] = -mag[i] * dvdx
+	        x[i3] = -mag[i] * dvdy
+	    } else
+		x[i] = -val
+	}
+
+	return (weight)
+end
+
+
+# DP_ACSHARP -- Procedure to accumulate sums of the weighted absolute
+# relative residuals and the image sharpness parameter for each of the
+# stars.
+
+procedure dp_acsharp (xcen, ycen, skip, npix, numer, denom, sumwt, chi,
+	group_size, fx, fy, fwhmx, fwhmy, ds, sigmasq, relerr, weight)
+
+real	xcen[ARB]		# array of object x centers
+real	ycen[ARB]		# array of object y centers
+int	skip[ARB]		# array of skip values
+int	npix[ARB]		# array of numbers of pixels
+real	numer[ARB]		# numerator array
+real	denom[ARB]		# denominator array
+real	sumwt[ARB]		# array of summed weights
+real	chi[ARB]		# array of chis
+int	group_size		# group size paramter.
+real	fx, fy			# position of data in image
+real	fwhmx, fwhmy		# gaussian core widths in x and y
+real	ds			# the data residual
+real	sigmasq			# the sigma squared
+real	relerr			# the relative error
+real	weight			# the weight
+
+int	i
+real	rhosq, dfdsig
+
+begin
+	do i = 1, group_size {
+
+	    if (skip[i] == YES)
+		next
+
+	    # Accumulate the number of pixels, chi and sum of the weights.
+	    npix[i] = npix[i] + 1
+	    chi[i] = chi[i] + relerr * weight
+	    sumwt[i] = sumwt[i] + weight
+
+	    # Include in the sharpness index only those pixels
+	    # within NCORE_SIGMASQ of the centroid of the
+	    # object. (This saves time and floating underflows
+	    # by excluding pixels which contribute very little
+	    # to the index.
+
+	    rhosq = ((xcen[i] - fx) / fwhmx) ** 2 + ((ycen[i] - fy) /
+	        fwhmy) ** 2
+	    if (rhosq > NCORE_SIGMASQ)
+		next
+	    rhosq = 0.6931472 * rhosq
+	    dfdsig = exp (-rhosq) * (rhosq - 1.0)
+	    numer[i] = numer[i] + dfdsig * ds / sigmasq
+	    denom[i] = denom[i] + (dfdsig ** 2) / sigmasq
+	}			
+end
+
+
+# DP_MATACCUM -- Procedure to accumulate the data into the matrices.
+
+procedure dp_mataccum (x, skip, group_size, c, v, cdimen, nterm, weight, dswt,
+	recenter, fitsky)
+
+real	x[ARB]			# x array
+int	skip[ARB]		# skip vector
+int	group_size		# size of the group
+real	c[cdimen,ARB]		# coefficient matrix
+real	v[ARB]			# vector array
+int	cdimen			# dimensions of the coefficient matrix
+int	nterm			# the number of terms
+real	weight			# weight
+real	dswt			# data weight
+int	recenter		# recenter the coordinates
+int	fitsky			# fit the sky value
+
+int	i, i3, i3m2, k, j, l
+
+begin
+	if (fitsky == YES) {
+	    c[nterm,nterm] = c[nterm,nterm] + weight
+	    v[nterm] = v[nterm] - dswt
+	}
+
+	do i = 1, group_size {
+	    if (skip[i] == YES)
+		next
+	    if (recenter == YES) {
+	        i3 = i * 3
+	        i3m2 = i3 - 2
+	        do k = i3m2, i3 {
+		    if (fitsky == YES)
+			c[nterm,k] = c[nterm,k] - x[k] * weight
+		    v[k] = v[k] + x[k] * dswt
+		}
+	        do j = 1, i {
+		    if (skip[j] == YES)
+		        next
+		    do k = i3m2, i3 {
+		        do l = 3 * j - 2, min (k, 3 * j)
+			    c[k,l] = c[k,l] + x[k] * x[l] * weight
+		    }
+	        }
+	    } else {
+		v[i] = v[i] + x[i] * dswt
+		if (fitsky == YES)
+			c[nterm,i] = c[nterm,i] - x[i] * weight
+		do j = 1, i {
+		    if (skip[j] == YES)
+		        next
+		    c[i,j] = c[i,j] + x[i] * x[j] * weight
+		}
+	    }
+	}
+end
+
+
+# DP_NTMIN -- Make sure that every star in the group has at least 
+# MIN_NPIX pixels within one fitting radius.
+
+bool procedure dp_ntmin (ids, xcen, ycen, mag, sky, npix, nier, group_size,
+	nterm, recenter, fitsky, groupsky, mean_sky, verbose)
+
+int	ids[ARB]		# array of ids
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of magnitudes
+real	sky[ARB]		# array of sky values
+int	npix[ARB]		# array of pixel numbers
+int	nier[ARB]		# array of error codes
+int	group_size		# size of the group
+int	nterm			# number of terms
+int	recenter		# recenter the objects
+int	fitsky			# fit the sky value
+int	groupsky		# use group sky value
+real	mean_sky		# the current mean sky value
+int	verbose			# verbose flag
+
+bool	redo
+int	i
+
+begin
+	redo = false
+	do i = 1, group_size {
+	    if (npix[i] >= MIN_NPIX)
+	        next
+	    redo = true
+	    if (verbose == YES) {
+		call printf (
+		    "\tStar %-5d has been deleted: too few good pixels\n")
+		    call pargi (ids[i])
+	    }
+	    if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1)
+	    call dp_remove (i, group_size, ids, xcen, ycen, mag, sky,
+		nier, NSTERR_NOPIX)
+	    if (group_size <= 0)
+		return (redo)
+	    if (recenter == YES)
+	        nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (fitsky == YES)
+		nterm = nterm + 1
+	}
+
+	return (redo)
+end
+
+
+# DP_MREFLECT -- Reflect the normal matrix around the diagonal.
+
+procedure dp_mreflect (c, cdimen, nterm)
+
+real	c[cdimen,ARB]		# coefficient matrix
+int	cdimen			# dimension of the c matrix
+int	nterm			# number of terms
+
+int	l, k
+
+begin
+	# Reflect the normal matrix across the diagonal.
+	do l = 2, nterm {
+	    do k = 1, l - 1
+		c[k,l] = c[l,k]
+	}
+end
+
+
+# DP_NTCHI -- Compute the chi value for each star.
+
+procedure dp_ntchi (sumwt, chi, group_size, ntmin, chigrp, grpwt)
+
+real	sumwt[ARB]		# sum of the weights
+real	chi[ARB]		# the chis:wq
+int	group_size		# size of the group
+int	ntmin			# minimum number of points
+real	chigrp			# the group chi
+real	grpwt			# the group weight
+
+int	i
+
+begin
+	do i = 1, group_size {
+	    if (sumwt[i] > ntmin) {
+		chi[i] = CHI_NORM * chi[i] / sqrt ((sumwt[i] - ntmin) *
+		    sumwt[i])
+		sumwt[i] = ((sumwt[i] - ntmin) * chi[i] + MIN_SUMWT) /
+		    sumwt[i]
+	    } else { 
+		chi[i] = chigrp
+		sumwt[i] = grpwt
+	    }
+	}
+end
+
+
+# DP_NCHECKC -- Check the inverted matrix for singularity.
+
+bool procedure dp_ncheckc (c, cdimen, nterm, ids, xcen, ycen, mag, sky,
+	nier, group_size, recenter, fitsky, groupsky, mean_sky, verbose)
+
+real	c[cdimen,ARB]		# coefficient matrix
+int	cdimen			# dimension of the c matrix
+int	nterm			# number of terms
+int	ids[ARB]		# array of ids
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of magnitudes
+real	sky[ARB]		# array of sky values
+int	nier[ARB]		# array of error codes
+int	group_size		# size of the group
+int	recenter		# recenter the objects
+int	fitsky			# fit the sky value
+int	groupsky		# use group sky value
+real	mean_sky		# the current mean sky value
+int	verbose			# verbose flag
+
+bool	redo
+int	j, starno, i
+
+begin
+	redo = false
+	do j = 1, nterm {
+	    if (c[j,j] > 0.0)
+		next
+	    redo = true
+	    if ((j == nterm) && (fitsky == YES))
+		starno = 0
+	    else if (recenter == YES)
+		starno = (j + 2) / 3
+	    else
+		starno = j
+	    if (starno == 0) {
+		if (verbose == YES) {
+		    do i = 1, group_size {
+		        call printf (
+			    "\tStar %-5d has been deleted: singular matrix\n")
+			    call pargi (ids[i])
+		    }
+		}
+		call amovki (NSTERR_SINGULAR, nier, group_size)
+		group_size = 0
+	    } else {
+		if (verbose == YES) {
+		    call printf (
+			"\tStar %-5d has been deleted: singular matrix\n")
+			call pargi (ids[starno])
+		}
+		if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		    mean_sky = (mean_sky * group_size - sky[starno]) /
+			(group_size - 1)
+		call dp_remove (starno, group_size, ids, xcen, ycen, mag,
+		    sky, nier, NSTERR_SINGULAR)
+	    }
+	    if (group_size <= 0)
+		return (redo)
+	    if (recenter == YES)
+		nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (fitsky == YES)
+		nterm = nterm + 1
+	}
+
+	return (redo)
+end
+
+
+# DP_NTCLAMP -- Restrict the amount the solution can vary on each iteration.
+
+procedure dp_ntclamp (xcen, ycen, mag, magerr, sumwt, group_size, c, cdimen,
+	x, xold, clamp, recenter, clip, redo)
+
+real	xcen[ARB]		# x centers array
+real	ycen[ARB]		# y centers array
+real	mag[ARB]		# magnitude array
+real	magerr[ARB]		# magnitude errors array
+real	sumwt[ARB]		# array of weight sums
+int	group_size		# size of the group
+real	c[cdimen, ARB]		# coefficient matrix
+int	cdimen			# dimensions of c
+real	x[ARB]			# x vector
+real	xold[ARB]		# old x vector
+real	clamp[ARB]		# clamp on solution matrix
+int	recenter		# recenter the objects
+bool	clip			# clip the matrix
+bool	redo			# redo the solution
+
+int	i, l, j, k
+real	df
+
+begin
+	do i = 1, group_size {
+
+
+	    # If any correction has changed sign since the last
+	    # iteration, reduce the maximum permissible change by
+	    # a factor of two.
+
+	    # Note that the sign of the correction is such that it
+	    # must be SUBTRACTED from the current value of the 
+	    # parameter to get the improved parameter value. This being
+	    # the case, if the correction to the brightness is
+	    # negative (the least-squares thinks that the star should
+	    # be brighter) a change of 1 magnitude is a change of a
+	    # factor of 2.5; if the brightness correction is positive
+	    # (the star should be fainter) a change of 1 magnitude
+	    # is a change of 60%.
+
+	    if (recenter == YES) {
+
+	        l = 3 * i
+	        k = l - 1
+	        j = l - 2
+
+	        if ((xold[j] * x[j]) < 0.0) 
+		    clamp[j] = 0.5 * clamp[j]
+	        mag[i] = mag[i] - x[j] / (1.0 + max (x[j] /
+		    (MAX_DELTA_FAINTER * mag[i]), -x[j] / (MAX_DELTA_BRIGHTER *
+		    mag[i])) / clamp[j])
+	        xold[j] = x[j]
+
+	        if ((xold[k] * x[k]) < 0.0)
+		    clamp[k] = 0.5 * clamp[k]
+	        if ((xold[l] * x[l]) < 0.0)
+		    clamp[l] = 0.5 * clamp[l]
+	        xcen[i] = xcen[i] - x[k] / (1.0 + abs(x[k]) / (clamp[k] *
+	            MAX_DELTA_PIX))
+	        ycen[i] = ycen[i] - x[l] / (1.0 + abs(x[l]) / (clamp[l] *
+	            MAX_DELTA_PIX))
+	        xold[k] = x[k]
+	        xold[l] = x[l]
+	        magerr[i] =sumwt[i] *  sqrt (c[j,j])
+
+	    } else {
+
+	        if ((xold[i] * x[i]) < 0.0) 
+		    clamp[i] = 0.5 * clamp[i]
+	        mag[i] = mag[i] - x[i] / (1.0 + max (x[i] /
+		    (MAX_DELTA_FAINTER * mag[i]), -x[i] / (MAX_DELTA_BRIGHTER *
+		    mag[i])) / clamp[i])
+	        xold[i] = x[i]
+	        magerr[i] =sumwt[i] *  sqrt (c[i,i])
+	    }
+
+
+	    # There are two milestones in the convergence process: the fits
+	    # proceed normally until each star's magnitude changes by less
+	    # than its standard error or MAX_NEW_ERRMAG magnitudes, whichever
+	    # is greater, and its x and y centroids change by less than 0.02 
+	    # pixel. At this point the least squares begins to apply
+	    # down-weighting of pixels with large residuals as described
+	    # above.  The fits then continue until each star's
+	    # magnitude changes by less than MAX (MAX_NEW_ERRMAG * std. error,
+	    # MAX_NEW_RELBRIGHT2 magnitude), ad its centroids change by
+	    # less than 0.002 pixel.
+
+	    if (redo)
+		next
+
+	    if (clip) {
+		if (abs (x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+		        MAX_NEW_RELBRIGHT2 * mag[i])) {
+		    redo = true
+		} else if (recenter == YES) {
+		    df = (MAX_NEW_ERRMAG * sumwt[i]) ** 2
+		    if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR2))
+			redo = true
+		    else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR2))
+			redo = true
+		}
+	    } else {
+		if (abs (x[j]) > max (magerr[i], MAX_NEW_RELBRIGHT1 *
+		        mag[i])) {
+		    redo = true
+		} else if (recenter == YES) {
+		    df = sumwt[i] ** 2
+		    if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR1))
+			redo = true
+		    else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR1))
+			redo = true
+		}
+	    }
+	}
+end
+
+
+# DP_NTCENTROID -- Check the new centroids to see if they have moved too
+# far off the edge of the image.
+
+procedure dp_ntcentroid (ids, xcen, ycen, mag, sky, nier, group_size, ixmin,
+	ixmax, iymin, iymax, fitradsq, fitsky, groupsky, mean_sky, redo,
+	verbose)
+
+int	ids[ARB]			# array of ids
+real	xcen[ARB]			# array of x centers
+real	ycen[ARB]			# array of y centers
+real	mag[ARB]			# array of magnitudes
+real	sky[ARB]			# array of sky values
+int	nier[ARB]			# array of error codes
+int	group_size			# size of the group
+int	ixmin,ixmax			# subraster x limits
+int	iymin,iymax			# subraster y limits
+real	fitradsq			# fit radius squared
+int	fitsky				# fit the sky value
+int	groupsky			# use the group sky value
+real	mean_sky			# the mean sky value
+bool	redo				# redo fit
+int	verbose				# verbose mode
+
+int	i
+real	dx, dy
+
+begin
+	# Check whether the centroid of any star has moved so far outside
+	# the picture that it has fewer than four or five pixels within
+	# one fitting radius.
+
+	do i = 1, group_size {
+
+	    # If the centroid of the star is outside the picture in x or
+	    # y, then DX or DY is its distance from the center of the edge 
+	    # pixel; otherwise DX and DY are zero.
+
+	    dx = max (ixmin - xcen[i], xcen[i] - ixmax, 0.0)
+	    dy = max (iymin - ycen[i], ycen[i] - iymax, 0.0)
+	    if ((dx <= MAX_PIX_INCREMENT) && (dy <= MAX_PIX_INCREMENT))
+		    next
+	    if (((dx + 1.0) ** 2 + (dy + 1.0) ** 2) < fitradsq)
+		next
+
+	    # Print a warning message about the star.
+	    if (verbose == YES) {
+		call printf (
+		"\tStar %-5d has been deleted: new center too far off image\n")
+		    call pargi (i)
+	    }
+
+	    # Adjust the sky.
+	    if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1)
+
+	    # Delete it.
+	    call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, nier,
+	        NSTERR_OFFIMAGE)
+	    if (group_size < 1)
+		break
+	    redo = true
+	}
+end
+
+
+# DP_NTFMAG -- Check for faint stars.
+
+procedure dp_ntfmag (mag, group_size, faint, ifaint)
+
+real	mag[ARB]		# array of magnitudes
+int	group_size		# size of the group
+real	faint			# faintest magnitude
+int	ifaint			# index of faintest magnitude
+
+int	i
+
+begin
+	faint = 1.0
+	ifaint = 0
+	do i = 1, group_size {
+	    if (mag[i] > MIN_REL_BRIGHT)
+		next
+	    if (mag[i] <= faint) {
+		faint = mag[i]
+		ifaint = i
+	    }
+	    mag[i] = MIN_REL_BRIGHT
+	}
+end
+
+
+# DP_FSNOISE -- Compute the smallest signal to noise ratio.
+
+procedure dp_fsnoise (mag, magerr, group_size, faint, ifaint)
+
+real	mag[ARB]		# array of magnitudes
+real	magerr[ARB]		# array of magnitude errors
+int	group_size		# size of group
+real	faint			# faint value
+int	ifaint			# faint index
+
+int	i
+real	weight
+
+begin
+	faint = 0.0
+	ifaint = 0
+	do i = 1, group_size {
+	    weight = magerr[i] / mag[i]
+	    if (weight < faint)
+		next
+	    faint = weight
+	    ifaint = i
+	}
+end
+
+
+# DP_REMOVE -- Remove the i-th star from the list of stars in the current
+# group by moving it to the end of the group.
+
+procedure dp_remove (i, nstar, ids, xcen, ycen, mag, sky, nier, pier)
+
+int	i			# the star to be removed
+int	nstar			# the number of stars in the group
+int	ids[ARB]		# the array of star ids
+real	xcen[ARB]		# the array of star x positions
+real	ycen[ARB]		# the array of star y positions
+real	mag[ARB]		# the array of magnitudes
+real	sky[ARB]		# the array of sky values.
+int	nier[ARB]		# array of error codes
+int	pier			# error code for deleted star
+
+int	ihold, phold
+real	xhold, yhold, shold, mhold
+
+begin
+	nier[i] = pier
+	if (i != nstar) {
+
+	    ihold = ids[i]
+	    xhold = xcen[i]
+	    yhold = ycen[i]
+	    shold = sky[i]
+	    mhold = mag[i]
+	    phold = nier[i]
+
+	    ids[i] = ids[nstar]
+	    xcen[i] = xcen[nstar]
+	    ycen[i] = ycen[nstar]
+ 	    sky[i] = sky[nstar]
+	    mag[i] = mag[nstar]
+	    nier[i] = nier[nstar]
+
+	    ids[nstar] = ihold
+	    xcen[nstar] = xhold
+	    ycen[nstar] = yhold
+	    sky[nstar] = shold
+	    mag[nstar] = mhold
+	    nier[nstar] = phold
+	}
+	nstar = nstar - 1
+end
diff --git a/noao/digiphot/daophot/nstar/dpntwrite.x b/noao/digiphot/daophot/nstar/dpntwrite.x
new file mode 100644
index 00000000..248a9f28
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpntwrite.x
@@ -0,0 +1,600 @@
+include <tbset.h>
+include <time.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/nstardef.h"
+
+
+# DP_TPNEWNSTAR -- Create a new NSTAR output ST table.
+
+procedure dp_tpnewnstar (dao, nst, colpoint)
+
+pointer	dao			# pointer to the daophot structure
+pointer	nst			# pointer to output photometry file
+pointer	colpoint[ARB]		# array  of column pointers
+
+
+int	i
+pointer	sp, colnames, colunits, colformat, col_dtype, col_len
+
+begin
+	# Allocate space for table definition.
+	call smark (sp)
+	call salloc (colnames, NST_NOUTCOL * (SZ_COLNAME + 1), TY_CHAR)
+	call salloc (colunits, NST_NOUTCOL * (SZ_COLUNITS + 1), TY_CHAR)
+	call salloc (colformat, NST_NOUTCOL * (SZ_COLFMT + 1), TY_CHAR)
+	call salloc (col_dtype, NST_NOUTCOL, TY_INT)
+	call salloc (col_len, NST_NOUTCOL, TY_INT)
+
+	# Set up the column definitions.
+	call strcpy (ID, Memc[colnames], SZ_COLNAME)
+	call strcpy (GROUP, Memc[colnames+SZ_COLNAME+1], SZ_COLNAME)
+	call strcpy (XCENTER, Memc[colnames+2*SZ_COLNAME+2], SZ_COLNAME)
+	call strcpy (YCENTER, Memc[colnames+3*SZ_COLNAME+3], SZ_COLNAME)
+	call strcpy (MAG, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME)
+	call strcpy (MAGERR, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME)
+	call strcpy (SKY, Memc[colnames+6*SZ_COLNAME+6], SZ_COLNAME)
+	call strcpy (NITER, Memc[colnames+7*SZ_COLNAME+7], SZ_COLNAME)
+	call strcpy (SHARP, Memc[colnames+8*SZ_COLNAME+8], SZ_COLNAME)
+	call strcpy (CHI, Memc[colnames+9*SZ_COLNAME+9], SZ_COLNAME)
+	call strcpy (PIER, Memc[colnames+10*SZ_COLNAME+10], SZ_COLNAME)
+	call strcpy (PERROR, Memc[colnames+11*SZ_COLNAME+11], SZ_COLNAME)
+
+	# Se up the column format definitions.
+	call strcpy ("%6d", Memc[colformat], SZ_COLFMT)
+	call strcpy ("%6d", Memc[colformat+SZ_COLFMT+1], SZ_COLFMT)
+	call strcpy ("%10.3f", Memc[colformat+2*SZ_COLFMT+2], SZ_COLFMT)
+	call strcpy ("%10.3f", Memc[colformat+3*SZ_COLFMT+3], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT)
+	call strcpy ("%14.3f", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT)
+	call strcpy ("%15.7g", Memc[colformat+6*SZ_COLFMT + 6], SZ_COLFMT)
+	call strcpy ("%6d", Memc[colformat+7*SZ_COLFMT+7], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+8*SZ_COLFMT+8], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+9*SZ_COLFMT+9], SZ_COLFMT)
+	call strcpy ("%6d", Memc[colformat+10*SZ_COLFMT+10], SZ_COLFMT)
+	call strcpy ("%13s", Memc[colformat+11*SZ_COLFMT+11], SZ_COLFMT)
+
+	# Set up the column unit definitions.
+	call strcpy ("NUMBER", Memc[colunits], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+SZ_COLUNITS+1], SZ_COLUNITS)
+	call strcpy ("PIXELS", Memc[colunits+2*SZ_COLUNITS+2], SZ_COLUNITS)
+	call strcpy ("PIXELS", Memc[colunits+3*SZ_COLUNITS+3], SZ_COLUNITS)
+	call strcpy ("MAGNITIDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS)
+	call strcpy ("MAGNITIDES", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS)
+	call strcpy ("COUNTS", Memc[colunits+6*SZ_COLUNITS+6], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+7*SZ_COLUNITS+7], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+8*SZ_COLUNITS+8], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+9*SZ_COLUNITS+9], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+10*SZ_COLUNITS+10], SZ_COLUNITS)
+	call strcpy ("PERRORS", Memc[colunits+11*SZ_COLUNITS+11], SZ_COLUNITS)
+
+	# Define the column datatypes.
+	Memi[col_dtype] = TY_INT
+	Memi[col_dtype+1] = TY_INT
+	Memi[col_dtype+2] = TY_REAL
+	Memi[col_dtype+3] = TY_REAL
+	Memi[col_dtype+4] = TY_REAL
+	Memi[col_dtype+5] = TY_REAL
+	Memi[col_dtype+6] = TY_REAL
+	Memi[col_dtype+7] = TY_INT
+	Memi[col_dtype+8] = TY_REAL
+	Memi[col_dtype+9] = TY_REAL
+	Memi[col_dtype+10] = TY_INT
+	Memi[col_dtype+11] = -13
+
+	# Define columnlengths.
+	do i = 1, NST_NOUTCOL 
+	    Memi[col_len+i-1] = 1
+	
+	# Define and create the table.
+	call tbcdef (nst, colpoint, Memc[colnames], Memc[colunits],
+	    Memc[colformat], Memi[col_dtype], Memi[col_len], NST_NOUTCOL)
+	call tbtcre (nst)
+
+	# Write out some header parameters.
+	call dp_tnstarpars (dao, nst)
+
+	call sfree (sp)
+	
+end
+
+
+define NST_NAME1STR "#N%4tID%10tGROUP%16tXCENTER%26tYCENTER%36tMAG%48t\
+MERR%62tMSKY%80t\\\n"
+define NST_UNIT1STR "#U%4t##%10t##%16tpixels%26tpixels%36tmagnitudes%48t\
+magnitudes%62tcounts%80t\\\n"
+define NST_FORMAT1STR "#F%4t%%-9d%10t%%-6d%16t%%-10.3f%26t%%-10.3f%36t\
+%%-12.3f%48t%%-14.3f%62t%%-15.7g%80t \n"
+
+define NST_NAME2STR "#N%12tNITER%18tSHARPNESS%30tCHI%42tPIER%48tPERROR%80t\\\n"
+define NST_UNIT2STR "#U%12t##%18t##%30t##%42t##%48tperrors%80t\\\n"
+define NST_FORMAT2STR "#F%12t%%-17d%18t%%-12.3f%30t%%-12.3f%42t%%-6d\
+%48t%%-13s%80t \n"
+
+# DP_XPNEWNSTAR -- Create a new NSTAR output text file.
+
+procedure dp_xpnewnstar (dao, nst)
+
+pointer	dao			# pointer to the daophot structure
+int	nst			# the output file descriptor
+
+
+begin
+	# Write out some header parameters.
+	call dp_xnstarpars (dao, nst)
+
+	# Write out the header banner.
+	call fprintf (nst, "#\n")
+	call fprintf (nst, NST_NAME1STR)
+	call fprintf (nst, NST_UNIT1STR)
+	call fprintf (nst, NST_FORMAT1STR)
+	call fprintf (nst, "#\n")
+	call fprintf (nst, NST_NAME2STR)
+	call fprintf (nst, NST_UNIT2STR)
+	call fprintf (nst, NST_FORMAT2STR)
+	call fprintf (nst, "#\n")
+end
+
+
+# DP_XNSTARPARS -- Add parameters to the header of the output NSTAR text file.
+
+procedure dp_xnstarpars (dao, nst)
+
+pointer	dao			# pointer to the daophot structure
+int	nst			# the output file descriptor
+
+pointer	psffit, sp, outstr, date, time, comment
+bool	itob()
+int	envfind()
+
+begin
+	# Get pointers
+	psffit = DP_PSFFIT(dao)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (outstr, SZ_LINE, TY_CHAR)
+	call salloc (date, SZ_DATE, TY_CHAR)
+	call salloc (time, SZ_DATE, TY_CHAR)
+	call salloc (comment, SZ_LINE, TY_CHAR)
+	Memc[comment] = EOS
+
+	# Write the id.
+	if (envfind ("version", Memc[outstr], SZ_LINE) <=0)
+	    call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE)
+	call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE)
+	call dp_sparam (nst, "IRAF", Memc[outstr], "version", Memc[comment])
+	if (envfind ("userid", Memc[outstr], SZ_LINE) > 0)
+	    call dp_sparam (nst, "USER", Memc[outstr], "name", Memc[comment])
+	call gethost (Memc[outstr], SZ_LINE)
+	call dp_sparam (nst, "HOST", Memc[outstr], "computer", Memc[comment])
+	call dp_date (Memc[date], Memc[time], SZ_DATE)
+	call dp_sparam (nst, "DATE", Memc[date], "yyyy-mm-dd", Memc[comment])
+	call dp_sparam (nst, "TIME", Memc[time], "hh:mm:ss", Memc[comment])
+	call dp_sparam (nst, "PACKAGE", "daophot", "name", Memc[comment])
+	call dp_sparam (nst, "TASK", "nstar", "name", Memc[comment])
+
+	# Write out the file names.
+	call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_FNAME)
+	call dp_sparam (nst, "IMAGE", Memc[outstr], "imagename",
+	    Memc[comment])
+	call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_FNAME)
+	call dp_sparam (nst, "GRPFILE", Memc[outstr], "filename", Memc[comment])
+	call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_FNAME)
+	call dp_sparam (nst, "PSFIMAGE", Memc[outstr], "imagename",
+	    Memc[comment])
+	call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME)
+	call dp_sparam (nst, "NSTARFILE", Memc[outstr], "filename",
+	    Memc[comment])
+	if (DP_OUTREJFILE(dao) == EOS)
+	    call dp_sparam (nst, "REJFILE", "\"\"", "filename", Memc[comment])
+	else {
+	    call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_FNAME)
+	    call dp_sparam (nst, "REJFILE", Memc[outstr], "filename",
+	        Memc[comment])
+	}
+
+	# Write out the data dependent parameters.
+	call dp_rparam (nst, "SCALE", DP_SCALE(dao), "units/pix",
+	    Memc[comment])
+	call dp_rparam (nst, "DATAMIN", DP_MINGDATA(dao), "counts",
+	    Memc[comment])
+	call dp_rparam (nst, "DATAMAX", DP_MAXGDATA(dao), "counts",
+	    Memc[comment])
+	call dp_rparam (nst, "GAIN", DP_PHOTADU(dao), "number", Memc[comment])
+	call dp_rparam (nst, "READNOISE", DP_READNOISE(dao), "electrons",
+	    Memc[comment])
+
+	# Write out the observing parameters.
+	call dp_sparam (nst, "OTIME", DP_OTIME(dao), "timeunit", Memc[comment])
+	call dp_rparam (nst, "XAIRMASS", DP_XAIRMASS(dao), "number",
+	    Memc[comment])
+	call dp_sparam (nst, "IFILTER", DP_IFILTER(dao), "filter",
+	    Memc[comment])
+
+	# Write out the fitting parameters.
+	call dp_bparam (nst, "RECENTER", itob (DP_RECENTER(dao)), "switch",
+	    Memc[comment])
+	call dp_bparam (nst, "FITSKY", itob (DP_FITSKY(dao)), "switch",
+	    Memc[comment])
+	call dp_bparam (nst, "GRPSKY", itob (DP_GROUPSKY(dao)), "switch",
+	    Memc[comment])
+	call dp_rparam (nst, "PSFMAG", DP_PSFMAG(psffit), "magnitude",
+	    Memc[comment])
+	call dp_rparam (nst, "PSFRAD", DP_SPSFRAD(dao), "scaleunit",
+	    Memc[comment])
+	call dp_rparam (nst, "FITRAD", DP_SFITRAD(dao), "scaleunit",
+	    Memc[comment])
+	call dp_iparam (nst, "MAXITER", DP_MAXITER(dao), "number",
+	    Memc[comment])
+	call dp_iparam (nst, "MAXGROUP", DP_MAXGROUP(dao), "number",
+	    Memc[comment])
+	call dp_rparam (nst, "FLATERROR", DP_FLATERR(dao), "percentage",
+	    Memc[comment])
+	call dp_rparam (nst, "PROFERROR", DP_PROFERR(dao), "percentage",
+	    Memc[comment])
+	call dp_iparam (nst, "CLIPEXP", DP_CLIPEXP(dao), "number",
+	    Memc[comment])
+	call dp_rparam (nst, "CLIPRANGE", DP_CLIPRANGE(dao), "sigma",
+	    Memc[comment])
+	call dp_rparam (nst, "MERGERAD", DP_SMERGERAD(dao), "scaleunit",
+	    Memc[comment])
+
+	call sfree(sp)
+end
+
+
+# DP_TNSTARPARS -- Add parameters to the header of the output NSTAR ST table.
+
+procedure dp_tnstarpars (dao, nst)
+
+pointer	dao			# pointer to the daophot structure
+pointer	nst			# pointer to the output photometry table
+
+pointer	psffit, sp, outstr, date, time
+bool	itob()
+int	envfind()
+
+begin
+	# Get pointers
+	psffit = DP_PSFFIT(dao)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (outstr, SZ_LINE, TY_CHAR)
+	call salloc (date, SZ_DATE, TY_CHAR)
+	call salloc (time, SZ_DATE, TY_CHAR)
+
+	# Write the id.
+	if (envfind ("version", Memc[outstr], SZ_LINE) <=0)
+	    call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE)
+	call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE)
+	call tbhadt (nst, "IRAF", Memc[outstr])
+	if (envfind ("userid", Memc[outstr], SZ_LINE) > 0)
+	    call tbhadt (nst, "USER", Memc[outstr])
+	call gethost (Memc[outstr], SZ_LINE)
+	call tbhadt (nst, "HOST", Memc[outstr])
+	call dp_date (Memc[date], Memc[time], SZ_DATE)
+	call tbhadt (nst, "DATE", Memc[date])
+	call tbhadt (nst, "TIME", Memc[time])
+	call tbhadt (nst, "PACKAGE", "daophot")
+	call tbhadt (nst, "TASK", "nstar")
+
+	# Write out the file names.
+	call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_FNAME)
+	call tbhadt (nst, "IMAGE", Memc[outstr])
+	call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_FNAME)
+	call tbhadt (nst, "GRPFILE", Memc[outstr])
+	call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_FNAME)
+	call tbhadt (nst, "PSFIMAGE", Memc[outstr])
+	call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME)
+	call tbhadt (nst, "NSTARFILE", Memc[outstr])
+	if (DP_OUTREJFILE(dao) == EOS)
+	    call tbhadt (nst, "REJFILE", "\"\"")
+	else {
+	    call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME)
+	    call tbhadt (nst, "REJFILE", Memc[outstr])
+	}
+
+	# Write out the data dependent parameters.
+	call tbhadr (nst, "SCALE", DP_SCALE(dao))
+	call tbhadr (nst, "DATAMIN", DP_MINGDATA(dao))
+	call tbhadr (nst, "DATAMAX", DP_MAXGDATA(dao))
+	call tbhadr (nst, "GAIN", DP_PHOTADU(dao))
+	call tbhadr (nst, "READNOISE", DP_READNOISE(dao))
+
+	# Write out the observing parameters.
+	call tbhadt (nst, "OTIME", DP_OTIME(dao))
+	call tbhadr (nst, "XAIRMASS", DP_XAIRMASS(dao))
+	call tbhadt (nst, "IFILTER", DP_IFILTER(dao))
+
+	# Write out the fitting parameters.
+	call tbhadb (nst, "RECENTER", itob (DP_RECENTER(dao)))
+	call tbhadb (nst, "FITSKY", itob (DP_FITSKY(dao)))
+	call tbhadb (nst, "GRPSKY", itob (DP_GROUPSKY(dao)))
+	call tbhadr (nst, "PSFMAG", DP_PSFMAG(psffit))
+	call tbhadr (nst, "PSFRAD", DP_SPSFRAD(dao))
+	call tbhadr (nst, "FITRAD", DP_SFITRAD(dao))
+	call tbhadi (nst, "MAXITER", DP_MAXITER(dao))
+	call tbhadi (nst, "MAXGROUP", DP_MAXGROUP(dao))
+	call tbhadr (nst, "FLATERROR", DP_FLATERR(dao))
+	call tbhadr (nst, "PROFERROR", DP_PROFERR(dao))
+	call tbhadi (nst, "CLIPEXP", DP_CLIPEXP(dao))
+	call tbhadr (nst, "CLIPRANGE", DP_CLIPRANGE(dao))
+	call tbhadr (nst, "MERGERAD", DP_SMERGERAD(dao))
+
+	call sfree(sp)
+end
+
+
+# DP_TNTWRITE -- Write out the NSTAR results to an ST table.
+
+procedure dp_tntwrite (dao, im, nst, rej, niter, old_size, output_row,
+	routput_row, colpoint) 
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+pointer	nst			# output photometry file descriptor
+int	rej			# output rejections file descriptor
+int	niter			# number of iterations
+int	old_size		# original size of group
+int	output_row		# output photometry file row number
+int	routput_row		# output rejections file row number
+pointer	colpoint[ARB]		# column pointer array
+
+int	i, id, nkeep, nreject, pier, plen, iter
+pointer	psffit, nstar, apsel, sp, perror
+real	xcen, ycen, mag, errmag, sharp
+int	dp_gnsterr()
+
+begin
+	# Get some daophot pointers.
+	nstar = DP_NSTAR(dao)
+	apsel = DP_APSEL(dao)
+	psffit = DP_PSFFIT(dao)
+
+	# Fill in the INDEFS.
+	nkeep = DP_NNUM(nstar)
+	nreject = old_size - nkeep
+	if (nreject  > 0) {
+	    call amovkr (INDEFR, Memr[DP_APMAG(apsel)+nkeep], nreject)
+	    call amovkr (INDEFR, Memr[DP_APERR(apsel)+nkeep], nreject)
+	    call amovkr (INDEFR, Memr[DP_APCHI(apsel)+nkeep], nreject)
+	}
+
+	call smark (sp)
+	call salloc (perror, SZ_FNAME, TY_CHAR)
+
+	# Now write out the results.
+	do i = 1, old_size {
+
+	    # Get the results.
+	    id = Memi[DP_APID (apsel)+i-1]
+	    xcen = Memr[DP_APXCEN (apsel)+i-1]
+	    ycen = Memr[DP_APYCEN (apsel)+i-1]
+	    if (IS_INDEFR(xcen) || IS_INDEFR(ycen))
+		next
+	    call dp_wout (dao, im, xcen, ycen, xcen, ycen, 1)
+	    mag = Memr[DP_APMAG(apsel)+i-1]
+	    errmag = Memr[DP_APERR(apsel)+i-1]
+	    if (! IS_INDEFR(mag)) {
+		if (! IS_INDEFR(errmag))
+	    	    errmag = 1.085736 * errmag / mag
+	        sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] *
+		    Memr[DP_PSFPARS(psffit)+1] * Memr[DP_NNUMER(nstar)+i-1] /
+		    (mag * DP_PSFHEIGHT(psffit) * Memr[DP_NDENOM(nstar)+i-1])
+	        if ((sharp < MIN_SHARP) || (sharp > MAX_SHARP))
+		    sharp = INDEFR
+	        mag = DP_PSFMAG (psffit) - 2.5 * log10 (mag)
+	    } else
+		sharp = INDEFR
+	    pier = Memi[DP_NIER(nstar)+i-1]
+	    if (pier == NSTERR_OK)
+		iter = niter
+	    else
+		iter = 0
+	    plen = dp_gnsterr (pier, Memc[perror], SZ_FNAME)
+
+	    # Write the results to the standard output.
+	    if (DP_VERBOSE (dao) == YES) {
+		call printf (
+		"\tID: %5d  XCEN: %8.2f  YCEN: %8.2f  MAG: %8.2f\n")
+		    call pargi (id)
+		    call pargr (xcen)
+		    call pargr (ycen)
+		    call pargr (mag)
+	    }
+
+	    # Write the output row to the proper table.
+	    if ((rej != NULL) && (pier != NSTERR_OK)) {
+	        routput_row = routput_row + 1
+	        call tbrpti (rej, colpoint[1], id, 1, routput_row)
+	        call tbrpti (rej, colpoint[2], DP_NGNUM(nstar), 1, routput_row)
+	        call tbrptr (rej, colpoint[3], xcen, 1, routput_row)
+	        call tbrptr (rej, colpoint[4], ycen, 1, routput_row)
+	        call tbrptr (rej, colpoint[5], mag, 1, routput_row)
+	        call tbrptr (rej, colpoint[6], errmag, 1, routput_row)
+	        call tbrptr (rej, colpoint[7], Memr[DP_APMSKY(apsel)+i-1],
+		    1, routput_row)
+	        call tbrpti (rej, colpoint[8], iter, 1, routput_row)
+	        call tbrptr (rej, colpoint[9], sharp, 1, routput_row)
+	        call tbrptr (rej, colpoint[10], Memr[DP_APCHI(apsel)+i-1],
+		    1, routput_row)
+	        call tbrpti (rej, colpoint[11], pier, 1, routput_row)
+		call tbrptt (rej, colpoint[12], Memc[perror], plen, 1,
+		    routput_row)
+	    } else {
+	        output_row = output_row + 1
+	        call tbrpti (nst, colpoint[1], id, 1, output_row)
+	        call tbrpti (nst, colpoint[2], DP_NGNUM(nstar), 1, output_row)
+	        call tbrptr (nst, colpoint[3], xcen, 1, output_row)
+	        call tbrptr (nst, colpoint[4], ycen, 1, output_row)
+	        call tbrptr (nst, colpoint[5], mag, 1, output_row)
+	        call tbrptr (nst, colpoint[6], errmag, 1, output_row)
+	        call tbrptr (nst, colpoint[7], Memr[DP_APMSKY(apsel)+i-1],
+		    1, output_row)
+	        call tbrpti (nst, colpoint[8], iter, 1, output_row)
+	        call tbrptr (nst, colpoint[9], sharp, 1, output_row)
+	        call tbrptr (nst, colpoint[10], Memr[DP_APCHI(apsel)+i-1],
+		    1, output_row)
+	        call tbrpti (nst, colpoint[11], pier, 1, output_row)
+		call tbrptt (nst, colpoint[12], Memc[perror], plen, 1,
+		    output_row)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+define NST_DATA1STR "%-9d%10t%-6d%16t%-10.3f%26t%-10.3f%36t%-12.3f%48t\
+%-14.3f%62t%-15.7g%80t\\\n"
+define NST_DATA2STR "%12t%-6d%18t%-12.3f%30t%-12.3f%42t%-6d%48t%-13.13s%80t \n"
+
+# DP_XNTWRITE -- Write out the NSTAR results to a text file.
+
+procedure dp_xntwrite (dao, im, nst, rej, niter, old_size) 
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+int	nst			# the output photometry file descriptor
+int	rej			# the output rejections file descriptor
+int	niter			# the number of the iteration
+int	old_size		# old size of group
+
+int	i, id, nkeep, nreject, pier, plen, iter
+pointer	nstar, psffit, apsel, sp, perror
+real	xcen, ycen, mag, errmag, sharp
+int	dp_gnsterr()
+
+begin
+	# Get some daophot pointers.
+	nstar = DP_NSTAR(dao)
+	psffit = DP_PSFFIT(dao)
+	apsel = DP_APSEL(dao)
+
+	# Fill in the results for the rejected stars with INDEFS.
+	nkeep = DP_NNUM(nstar)
+	nreject = old_size - nkeep
+	if (nreject  > 0) {
+	    call amovkr (INDEFR, Memr[DP_APMAG(apsel)+nkeep], nreject)
+	    call amovkr (INDEFR, Memr[DP_APERR(apsel)+nkeep], nreject)
+	    call amovkr (INDEFR, Memr[DP_APCHI(apsel)+nkeep], nreject)
+	}
+
+	call smark (sp)
+	call salloc (perror, SZ_FNAME, TY_CHAR)
+
+	# Now write out the results.
+	do i = 1, old_size {
+
+	    # Get the results.
+	    id = Memi[DP_APID (apsel)+i-1]
+	    xcen = Memr[DP_APXCEN (apsel)+i-1]
+	    ycen = Memr[DP_APYCEN (apsel)+i-1]
+	    if (IS_INDEFR(xcen) || IS_INDEFR(ycen))
+		next
+	    call dp_wout (dao, im, xcen, ycen, xcen, ycen, 1)
+	    mag = Memr[DP_APMAG(apsel)+i-1]
+	    errmag = Memr[DP_APERR(apsel)+i-1]
+	    if (! IS_INDEFR(mag)) {
+		if (! IS_INDEFR(errmag))
+	            errmag = 1.085736 * errmag / mag
+	        sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] *
+		    Memr[DP_PSFPARS(psffit)+1] * Memr[DP_NNUMER(nstar)+i-1] /
+		    (mag * DP_PSFHEIGHT(psffit) * Memr[DP_NDENOM(nstar)+i-1])
+	        if ((sharp < MIN_SHARP) || (sharp > MAX_SHARP))
+		    sharp = INDEFR
+	        mag = DP_PSFMAG (psffit) - 2.5 * log10 (mag)
+	    } else
+		sharp = INDEFR
+	    pier = Memi[DP_NIER(nstar)+i-1]
+	    if (pier == NSTERR_OK)
+		iter = niter
+	    else
+		iter = 0
+	    plen = dp_gnsterr (pier, Memc[perror], SZ_FNAME)
+
+	    # Write the results to the standard output.
+	    if (DP_VERBOSE (dao) == YES) {
+		call printf (
+		"\tID: %5d  XCEN: %8.2f  YCEN: %8.2f  MAG: %8.2f\n")
+	            call pargi (id)
+		    call pargr (xcen)
+		    call pargr (ycen)
+		    call pargr (mag)
+	    }
+
+	    # Write the results to the output file.
+	    if ((rej != NULL) && (pier != NSTERR_OK)) {
+	        call fprintf (rej, NST_DATA1STR)
+		    call pargi (id)
+		    call pargi (DP_NGNUM(nstar))
+		    call pargr (xcen)
+		    call pargr (ycen)
+	            call pargr (mag)
+		    call pargr (errmag)
+	            call pargr (Memr[DP_APMSKY(apsel)+i-1])
+	        call fprintf (rej, NST_DATA2STR)
+		    call pargi (iter)
+	            call pargr (sharp)
+	            call pargr (Memr[DP_APCHI(apsel)+i-1])
+		    call pargi (pier)
+		    call pargstr (Memc[perror])
+	    } else {
+	        call fprintf (nst, NST_DATA1STR)
+		    call pargi (id)
+		    call pargi (DP_NGNUM(nstar))
+		    call pargr (xcen)
+		    call pargr (ycen)
+	            call pargr (mag)
+		    call pargr (errmag)
+	            call pargr (Memr[DP_APMSKY(apsel)+i-1])
+	        call fprintf (nst, NST_DATA2STR)
+		    call pargi (iter)
+	            call pargr (sharp)
+	            call pargr (Memr[DP_APCHI(apsel)+i-1])
+		    call pargi (pier)
+		    call pargstr (Memc[perror])
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# DP_GNSTERR -- Set the NSTAR task error code string.
+ 
+int procedure dp_gnsterr (ier, perror, maxch)
+ 
+int     ier             # the integer error code
+char    perror          # the output error code string
+int     maxch           # the maximum size of the error code string
+ 
+int     plen
+int     gstrcpy()
+ 
+begin
+        switch (ier) {
+        case NSTERR_OK:
+            plen = gstrcpy ("No_error", perror, maxch)
+	case NSTERR_BIGGROUP:
+            plen = gstrcpy ("Big_group", perror, maxch)
+        case NSTERR_INDEFSKY:
+            plen = gstrcpy ("Bad_sky", perror, maxch)
+        case NSTERR_NOPIX:
+            plen = gstrcpy ("Npix_too_few", perror, maxch)
+        case NSTERR_SINGULAR:
+            plen = gstrcpy ("Singular", perror, maxch)
+        case NSTERR_FAINT:
+            plen = gstrcpy ("Too_faint", perror, maxch)
+        case NSTERR_MERGE:
+            plen = gstrcpy ("Merged", perror, maxch)
+        case NSTERR_OFFIMAGE:
+            plen = gstrcpy ("Off_image", perror, maxch)
+        default:
+            plen = gstrcpy ("No_error", perror, maxch)
+        }
+ 
+        return (plen)
+end
diff --git a/noao/digiphot/daophot/nstar/mkpkg b/noao/digiphot/daophot/nstar/mkpkg
new file mode 100644
index 00000000..327f3487
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/mkpkg
@@ -0,0 +1,24 @@
+# NSTAR task
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	dpggroup.x	"../../lib/ptkeysdef.h"	../lib/daophotdef.h	  \
+			../lib/apseldef.h
+	dpmemnstar.x	../lib/daophotdef.h	../lib/nstardef.h
+	dpnconfirm.x	../lib/daophotdef.h
+	dpnstar.x	<imhdr.h>		<tbset.h>                 \
+			<mach.h>                ../lib/daophotdef.h       \
+			../lib/apseldef.h       ../lib/nstardef.h
+	dpnstarfit.x	<imhdr.h>		../lib/daophotdef.h       \
+			../lib/apseldef.h	../lib/nstardef.h         \
+			<mach.h>
+	dpntwrite.x	<tbset.h>		<time.h>		  \
+			../lib/daophotdef.h	../lib/apseldef.h         \
+			../lib/nstardef.h
+	t_nstar.x	<fset.h>		<imhdr.h>                 \
+	                ../lib/daophotdef.h
+	;
diff --git a/noao/digiphot/daophot/nstar/t_nstar.x b/noao/digiphot/daophot/nstar/t_nstar.x
new file mode 100644
index 00000000..8f6c542e
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/t_nstar.x
@@ -0,0 +1,308 @@
+include	<fset.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+
+# T_NSTAR  -- Procedure to fit the PSF to multiple stars.
+
+procedure t_nstar ()
+
+pointer	image				# input image descriptor
+pointer	groupfile			# input group file descriptor
+pointer	psfimage			# input PSF image descriptor
+pointer	nstarfile			# output photometry file descriptor
+pointer	rejfile				# output rejections file descriptor
+int	verbose				# print messages 
+int	verify				# verify the critical parameters
+int	update				# update the parameter set
+int	cache				# cache the input image pixels
+
+pointer	sp, outfname, im, psfim, dao, str
+int	imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist
+int	rlist, lrlist, root, grp, nst, rejfd, wcs, req_size, old_size
+int	buf_size, memstat
+bool	ap_text
+
+pointer	immap(), tbtopn()
+int	strlen(), strncmp(), fnldir(), fstati(), open(), btoi()
+int	access(), imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb()
+int	fntgfnb(), clgwrd(), sizeof(), dp_memstat()
+bool	clgetb(), itob()
+
+begin
+	# Set the standard output to flush on newline.
+	if (fstati (STDOUT, F_REDIR) == NO)
+	    call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (groupfile, SZ_FNAME, TY_CHAR)
+	call salloc (psfimage, SZ_FNAME, TY_CHAR)
+	call salloc (nstarfile, SZ_FNAME, TY_CHAR)
+	call salloc (rejfile, SZ_FNAME, TY_CHAR)
+	call salloc (outfname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the input and output file names.
+	call clgstr ("image", Memc[image], SZ_FNAME)
+	call clgstr ("groupfile", Memc[groupfile], SZ_FNAME)
+	call clgstr ("psfimage", Memc[psfimage], SZ_FNAME)
+	call clgstr ("nstarfile", Memc[nstarfile], SZ_FNAME)
+	call clgstr ("rejfile", Memc[rejfile], SZ_FNAME)
+
+	# Get the task mode parameters.
+	verbose = btoi (clgetb ("verbose"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+	cache = btoi (clgetb ("cache"))
+
+	# Get the lists.
+	imlist = imtopen (Memc[image])
+	limlist = imtlen (imlist)
+	alist = fntopnb (Memc[groupfile], NO)
+	lalist = fntlenb (alist)
+	pimlist = imtopen (Memc[psfimage])
+	lpimlist = imtlen (pimlist)
+	olist = fntopnb (Memc[nstarfile], NO)
+	lolist = fntlenb (olist)
+	rlist = fntopnb (Memc[rejfile], NO)
+	lrlist = fntlenb (rlist)
+
+	# Test that the lengths of the photometry file, psf image, and
+	# output file lists are the same as the length of the input image
+	# list.
+
+	if ((limlist != lalist) && (strncmp (Memc[groupfile], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and group file list lengths")
+	}
+
+	if ((limlist != lpimlist) && (strncmp (Memc[psfimage], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and psf file list lengths")
+	}
+
+	if ((limlist != lolist) && (strncmp (Memc[nstarfile], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and nstar file list lengths")
+	}
+
+	if ((lrlist > 0) && (limlist != lrlist) && (strncmp (Memc[rejfile],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and rejections file list lengths")
+	}
+
+	# Open the daophot structure and get some parameters.
+	call dp_gppars (dao)	
+
+	# Set some parameters.
+	call dp_seti (dao, VERBOSE, verbose)
+
+	# Verify and update the parameters if appropriate.
+	if (verify == YES) {
+	    call dp_nconfirm (dao)
+	    if (update == YES)
+		call dp_pppars (dao)
+	}
+
+        # Get the wcs information.
+        wcs = clgwrd ("wcsin", Memc[str], SZ_FNAME, WCSINSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the input coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSIN, wcs)
+        wcs = clgwrd ("wcsout", Memc[str], SZ_FNAME, WCSOUTSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the output coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSOUT, wcs)
+        wcs = clgwrd ("wcspsf", Memc[str], SZ_FNAME, WCSPSFSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the psf coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSPSF, wcs)
+
+	# Initialize the photometry structure.
+	call dp_apsetup (dao)
+
+	# Initialize the PSF structure.
+	call dp_fitsetup (dao)
+
+	# Initialize the nstar structure.
+	call dp_nstarsetup (dao)
+
+	# Loop over the images.
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open the input image.
+	    im = immap (Memc[image], READ_ONLY, 0)		
+	    call dp_imkeys (dao, im)
+	    call dp_sets (dao, INIMAGE, Memc[image])
+
+            # Cache the input image pixels.
+            req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) *
+                sizeof (IM_PIXTYPE(im))
+            memstat = dp_memstat (cache, req_size, old_size)
+            if (memstat == YES)
+                call dp_pcache (im, INDEFI, buf_size)
+
+	    # Open the input group table.
+	    if (fntgfnb (alist, Memc[groupfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[groupfile], SZ_FNAME)
+	    root = fnldir (Memc[groupfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[groupfile+root],
+	        DEF_LENDEFNAME) == 0 || (root == strlen (Memc[groupfile])))
+	        call dp_inname (Memc[image], Memc[outfname], "grp",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	    	call strcpy (Memc[groupfile], Memc[outfname], SZ_FNAME)
+	    ap_text = itob (access (Memc[outfname], 0, TEXT_FILE))
+	    if (ap_text)
+	        grp = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+	    else
+	        grp = tbtopn (Memc[outfname], READ_ONLY, 0)
+	    call dp_sets (dao, INPHOTFILE, Memc[outfname])
+
+	    # Open and read the PSF image.
+	    if (imtgetim (pimlist, Memc[psfimage], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[psfimage], SZ_FNAME)
+	    root = fnldir (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[psfimage+root],
+	        DEF_LENDEFNAME) == 0 || (root == strlen (Memc[psfimage])))
+	        call dp_iimname (Memc[image], Memc[outfname], "psf",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    psfim = immap (Memc[outfname], READ_ONLY, 0)
+	    call dp_readpsf (dao, psfim)
+	    call dp_sets (dao, PSFIMAGE, Memc[outfname])
+	
+	    # Open the output NSTAR file. If the output is DEF_DEFNAME,
+	    # dir$default or a directory specification then the extension
+	    # "nst" is added to the image name and a suitable version number
+	    # is appended to the output name.
+
+	    if (fntgfnb (olist, Memc[nstarfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[nstarfile], SZ_FNAME)
+	    root = fnldir (Memc[nstarfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[nstarfile+root],
+	        DEF_LENDEFNAME) == 0 || (root == strlen (Memc[nstarfile]))) {
+	        call dp_outname (Memc[image], Memc[outfname], "nst",
+		    Memc[outfname], SZ_FNAME)
+	    } else
+	        call strcpy (Memc[nstarfile], Memc[outfname], SZ_FNAME)
+	    if (DP_TEXT(dao) == YES)
+	        nst = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        nst = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTPHOTFILE, Memc[outfname])
+
+	    if (lrlist <= 0) {
+		rejfd = NULL
+		Memc[outfname] = EOS
+	    } else {
+	        if (fntgfnb (rlist, Memc[rejfile], SZ_FNAME) == EOF)
+		    call strcpy (DEF_DEFNAME, Memc[rejfile], SZ_FNAME)
+	        root = fnldir (Memc[rejfile], Memc[outfname], SZ_FNAME)
+	        if (strncmp (DEF_DEFNAME, Memc[rejfile+root],
+		    DEF_LENDEFNAME) == 0 || (root == strlen (Memc[rejfile])))
+	            call dp_outname (Memc[image], Memc[outfname], "nrj",
+		        Memc[outfname], SZ_FNAME)
+	        else
+	            call strcpy (Memc[rejfile], Memc[outfname], SZ_FNAME)
+	        if (DP_TEXT(dao) == YES)
+	            rejfd = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	        else
+	            rejfd = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    }
+	    call dp_sets (dao, OUTREJFILE, Memc[outfname])
+
+	    # Do the PSF fitting.
+	    call dp_nphot (dao, im, grp, nst, rejfd, ap_text)
+
+	    # Close the input image.
+	    call imunmap (im)
+
+	    # Close the group file. 
+	    if (ap_text)
+		call close (grp)
+	    else
+	        call tbtclo (grp)
+
+	    # Close the PSF image.
+	    call imunmap (psfim)
+
+	    # Close the output photometry file.
+	    if (DP_TEXT(dao) == YES)
+		call close (nst)
+	    else
+	        call tbtclo (nst)
+
+	    # Close the output rejections file.
+	    if (rejfd != NULL) {
+	        if (DP_TEXT(dao) == YES)
+		    call close (rejfd)
+	        else
+	            call tbtclo (rejfd)
+	    }
+
+            # Uncache memory.
+            call fixmem (old_size)
+
+	}
+
+	# Close the image/file lists.
+	call imtclose (imlist)
+	call fntclsb (alist)
+	call imtclose (pimlist)
+	call fntclsb (olist)
+	call fntclsb (rlist)
+
+	# Close the nstar structure.
+	call dp_nsclose (dao)
+
+	# Close the PSF structure.
+	call dp_fitclose (dao)
+
+	# Close the photometry structure.
+	call dp_apclose (dao)
+
+	# Free the daophot structure.
+	call dp_free (dao)
+
+	call sfree(sp)
+end