Initial commit

1 files changed, 530 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/daolib/dpgetapert.x b/noao/digiphot/daophot/daolib/dpgetapert.x
new file mode 100644
index 00000000..8268373f
--- /dev/null
+++ b/noao/digiphot/daophot/daolib/dpgetapert.x
@@ -0,0 +1,530 @@
+include	<tbset.h>
+include "../../lib/ptkeysdef.h"
+include	"../lib/daophotdef.h"
+include	"../lib/apseldef.h"
+
+# DP_WGETAPERT -- Read the aperture photometry results and transform the
+# input coordinates to the appropriate coordinate system. Works with
+# either the "old" APPHOT files or the new ST Tables.
+
+procedure dp_wgetapert (dao, im, apd, max_nstars, old_ap)
+
+pointer	dao			# pointer to the DAOPHOT structure
+pointer	im			# the input image descriptor
+int	apd			# input photometry file descriptor
+int	max_nstars		# maximum number of stars
+bool	old_ap			# YES indicates old APPHOT file
+
+pointer	apsel
+int	dp_stati()
+
+begin
+	# Get the stars.
+	call dp_getapert (dao, apd, max_nstars, old_ap)
+
+	# Transform the coordinates if necessary.
+	apsel = DP_APSEL (dao)
+	if (dp_stati (dao, WCSIN) != WCS_LOGICAL)
+	    call dp_win (dao, im, Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], DP_APNUM(apsel))
+end
+
+
+# DP_GETAPERT -- Read the aperture photometry results. Works with
+# either the "old" APPHOT files or the new ST Tables.
+
+procedure dp_getapert (dao, apd, max_nstars, old_ap)
+
+pointer	dao			# pointer to the DAOPHOT structure
+int	apd			# input Photometry file descriptor
+int	max_nstars		# maximum number of stars
+bool	old_ap			# YES indicates old APPHOT file
+
+int	nstars
+pointer	apsel
+int	tbpsta(), dp_goldap(), dp_gtabphot()
+
+begin
+	# Get APSEL pointer.
+	apsel = DP_APSEL (dao)
+
+	# Get the required memory.
+	Memi[DP_APRESULT(apsel)] = DP_PAPID
+	Memi[DP_APRESULT(apsel)+1] = DP_PAPXCEN
+	Memi[DP_APRESULT(apsel)+2] = DP_PAPYCEN
+	Memi[DP_APRESULT(apsel)+3] = DP_PAPMAG1
+	Memi[DP_APRESULT(apsel)+4] = DP_PAPSKY
+
+	# Get the results.
+	if (old_ap) {
+	    call dp_memapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT,
+	        max_nstars)
+	    nstars = dp_goldap (apd, dao, max_nstars)
+	} else {
+	    call dp_memapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT,
+	        tbpsta (apd, TBL_NROWS))
+	    nstars = dp_gtabphot (apd, dao, tbpsta (apd, TBL_NROWS))
+	}
+
+	# Reallocate to save space if appropopriate.
+	if (nstars < max_nstars)
+	    call dp_rmemapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT,
+	        nstars)
+end
+
+
+# DP_MEMAPSEL -- Procedure to allocate memory for the apselect strucuture.
+
+procedure dp_memapsel (dao, fields, nfields, max_nstars)
+
+pointer	dao		# pointer to the daophot strucuture
+int	fields[ARB]	# array of fields
+int	nfields		# number of fields to allocate space for
+int	max_nstars	# maximum number of stars
+
+int	i
+pointer	apsel
+
+begin
+	apsel = DP_APSEL(dao)
+
+	# Allocate space for results.
+	do i = 1, nfields {
+	    switch (fields[i]) {
+	    case DP_PAPID:
+	        if (DP_APID(apsel) != NULL)
+	            call mfree (DP_APID(apsel), TY_INT)
+	        call malloc (DP_APID(apsel), max_nstars, TY_INT)
+
+	    case DP_PAPXCEN:
+	        if (DP_APXCEN(apsel) != NULL)
+	            call mfree (DP_APXCEN(apsel), TY_REAL)
+	        call malloc (DP_APXCEN(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPYCEN:
+	        if (DP_APYCEN(apsel) != NULL)
+	            call mfree (DP_APYCEN(apsel), TY_REAL)
+	        call malloc (DP_APYCEN(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPSKY:
+	        if (DP_APMSKY(apsel) != NULL)
+	            call mfree (DP_APMSKY(apsel), TY_REAL)
+	        call malloc (DP_APMSKY(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPMAG1:
+		if (DP_APMAG(apsel) != NULL)
+	    	    call mfree (DP_APMAG(apsel), TY_REAL)
+		call malloc (DP_APMAG(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPGROUP:
+		if (DP_APGROUP(apsel) != NULL)
+	    	    call mfree (DP_APGROUP(apsel), TY_INT)
+		#call malloc (DP_APGROUP(apsel), max_nstars, TY_INT)
+
+	    case DP_PAPMERR1:
+		if (DP_APERR(apsel) != NULL)
+	    	    call mfree (DP_APERR(apsel), TY_REAL)
+		call malloc (DP_APERR(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPNITER:
+		if (DP_APNITER(apsel) != NULL)
+	    	    call mfree (DP_APNITER(apsel), TY_INT)
+		#call malloc (DP_APNITER(apsel), max_nstars, TY_INT)
+
+	    case DP_PAPCHI:
+		if (DP_APCHI(apsel) != NULL)
+	    	    call mfree (DP_APCHI(apsel), TY_REAL)
+		call malloc (DP_APCHI(apsel), max_nstars, TY_REAL)
+
+	    case DP_PAPSHARP:
+		if (DP_APSHARP(apsel) != NULL)
+	    	    call mfree (DP_APSHARP(apsel), TY_REAL)
+		call malloc (DP_APSHARP(apsel), max_nstars, TY_REAL)
+	    }
+	}
+end
+
+
+# DP_GOLDAP -- Read in the photometry from an old style APPHOT file
+
+int procedure  dp_goldap (apd, dao, max_nstars)
+
+int	apd		# the input file descriptor
+pointer	dao		# pointer to the daophot structure
+int	max_nstars	# maximum number of stars
+
+int	nstars, bufsize, stat
+pointer	apsel, apkey, sp, fields
+int	dp_apsel()
+
+begin
+	# Define the point to the apselect structure.
+	apsel = DP_APSEL (dao)
+
+	# Allocate some temporary space.
+	call smark (sp)
+	call salloc (fields, SZ_LINE, TY_CHAR)
+
+	# Initialize the keyword structure.
+	call pt_kyinit (apkey)
+
+	# Set up the fields to be retrieved.
+	call dp_gappsf (Memi[DP_APRESULT(apsel)], Memc[fields], NAPRESULT)
+
+	# Now read in the results.
+	nstars = 0
+	bufsize = max_nstars
+	repeat {
+
+	    # Read in a group of stars.
+	    while (nstars < bufsize) {
+	        stat = dp_apsel (apkey, apd, Memc[fields],
+		    Memi[DP_APRESULT(apsel)], Memi[DP_APID(apsel)+nstars],
+	            Memr[DP_APXCEN(apsel)+nstars],
+		    Memr[DP_APYCEN(apsel)+nstars],
+	            Memr[DP_APMSKY(apsel)+nstars],
+		    Memr[DP_APMAG(apsel)+nstars])
+		if (stat == EOF)
+		    break
+	        nstars = nstars + 1
+	    }
+
+	    # Check the buffer size.
+	    if (stat == EOF)
+		break
+	    bufsize = bufsize + max_nstars
+	    call dp_rmemapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT,
+	        bufsize)
+
+	}
+	DP_APNUM(apsel) = nstars
+
+	# Free the keyword structure.
+	call pt_kyfree (apkey)
+	call sfree (sp)
+
+	return (nstars)
+end
+
+
+# DP_GTABPHOT -- Read in the complete photometry from an ST table.
+
+int procedure dp_gtabphot (tp, dao, max_nstars)
+
+pointer	tp			# table descriptor
+pointer	dao			# pointer to daophot structure
+int	max_nstars		# maximum number of stars
+
+bool	nullflag
+int	record, index, nrow
+pointer	apsel, idpt, xcenpt, ycenpt, magpt, skypt
+int	tbpsta()
+
+begin
+	# Define the point to the apselect structure.
+	apsel = DP_APSEL (dao)
+
+	# Find the column pointers
+	call tbcfnd (tp, ID, idpt, 1)
+	if (idpt == NULL)
+	    call tbcfnd (tp, "ID", idpt, 1)
+	if (idpt == NULL)
+	    call printf ("Error reading ID.\n")
+
+	call tbcfnd (tp, XCENTER, xcenpt, 1)
+	if (xcenpt == NULL)
+	    call tbcfnd (tp, "XCENTER", xcenpt, 1)
+	if (xcenpt == NULL)
+	    call printf ("Error reading XCENTER.\n")
+
+	call tbcfnd (tp, YCENTER, ycenpt, 1)
+	if (ycenpt == NULL)
+	    call tbcfnd (tp, "YCENTER", ycenpt, 1)
+	if (ycenpt == NULL)
+	    call printf ("Error reading YCENTER.\n")
+
+	call tbcfnd (tp, MAG, magpt, 1)
+	if (magpt == NULL)
+	    call tbcfnd (tp, APMAG, magpt, 1)
+   	if (magpt == NULL)
+	    call printf ("Error reading MAG.\n")
+
+	call tbcfnd (tp, SKY, skypt, 1)
+	if (skypt == NULL)
+	    call tbcfnd (tp, SKY, skypt, 1)
+	if (skypt == NULL)
+	    call printf ("Error reading SKY.\n")
+
+
+	# Get the results ignoring any record with ID = NULL.
+	nrow = min (tbpsta (tp, TBL_NROWS), max_nstars)
+	index = 0
+	do record = 1, nrow {
+
+	    # Check the ID record.
+	    call tbrgti (tp, idpt, Memi[DP_APID(apsel)+index], nullflag, 1,
+	        record)
+	    if (nullflag)
+		next
+
+	    # Read the remaining records.
+	    call tbrgtr (tp, xcenpt, Memr[DP_APXCEN(apsel)+index], nullflag, 1,
+	        record)
+	    call tbrgtr (tp, ycenpt, Memr[DP_APYCEN(apsel)+index], nullflag, 1,
+	        record)
+	    call tbrgtr (tp, magpt, Memr[DP_APMAG(apsel)+index], nullflag, 1,
+	        record)
+	    call tbrgtr (tp, skypt, Memr[DP_APMSKY(apsel)+index], nullflag, 1,
+	        record)
+	    index = index + 1
+	}
+
+	DP_APNUM(apsel) = index
+
+	return (index)
+end
+
+
+# DP_RMEMAPSEL -- Procedure to reallocate memory for the apselect strucuture.
+
+procedure dp_rmemapsel (dao, fields, nfields, max_nstars)
+
+pointer	dao		# pointer to the daophot strucuture
+int	fields[ARB]	# integer fields
+int	nfields		# number of fields
+int	max_nstars	# maximum number of stars
+
+int	i
+pointer	apsel
+
+begin
+	# Reallocate space for results.
+	apsel = DP_APSEL(dao)
+
+	do i = 1,  nfields {
+	    switch (fields[i]) {
+	    case DP_PAPID:
+	        call realloc (DP_APID(apsel), max_nstars, TY_INT)
+	    case DP_PAPXCEN:
+	        call realloc (DP_APXCEN(apsel), max_nstars, TY_REAL)
+	    case DP_PAPYCEN:
+	        call realloc (DP_APYCEN(apsel), max_nstars, TY_REAL)
+	    case DP_PAPSKY:
+	        call realloc (DP_APMSKY(apsel), max_nstars, TY_REAL)
+	    case DP_PAPGROUP:
+		#call realloc (DP_APGROUP(apsel), max_nstars, TY_INT)
+	    case DP_PAPMAG1:
+	        call realloc (DP_APMAG(apsel), max_nstars, TY_REAL)
+	    case DP_PAPMERR1:
+	        call realloc (DP_APERR(apsel), max_nstars, TY_REAL)
+	    case DP_PAPNITER:
+	        #call realloc (DP_APNITER(apsel), max_nstars, TY_INT)
+	    case DP_PAPSHARP:
+	        call realloc (DP_APSHARP(apsel), max_nstars, TY_REAL)
+	    case DP_PAPCHI:
+	        call realloc (DP_APCHI(apsel), max_nstars, TY_REAL)
+	    }
+	}
+end
+
+
+# DP_GAPPSF -- Set up the structures necessary for retrieving the
+# aperture phtometery results needed for the PSF task.
+
+procedure dp_gappsf (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
+	sel_fields[1] = EOS
+
+	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 (APMAG)
+	    }
+	}
+
+	if (sel_fields[1] != EOS)
+	    sel_fields[strlen(sel_fields)] = EOS
+end
+
+
+# DP_APSEL -- Select records from an apphot/daophot text file.
+
+int procedure dp_apsel (key, fd, fields, indices, id, 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 of fields
+int	id		# star id 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)
+	}
+
+	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 table 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 <= 0)
+			    break
+		    }
+
+		    # Construct the output record by moving selected fields
+		    # into data structures.
+		    call dp_getap (key, indices, id, x, y, sky, mag)
+		    record = record + 1
+		    first_rec = NO
+
+		    # Record is complete so exist the loop.
+		    break
+	        }
+	    }
+
+	}
+
+	if (nchars == EOF || nselect <= 0) {
+	    first_rec = YES
+	    nunique = 0
+	    uunique = 0
+	    funique = 0
+	    nselect = 0
+	    call mfree (line, TY_CHAR)
+	    return (EOF)
+	} else
+	    return (record)
+end
+
+
+# DP_GETAP  -- Decode the selected daophot/apphot values.
+
+procedure dp_getap (key, indices, id, x, y, sky, mag)
+
+pointer	key		# pointer to keys strucuture
+int	indices[ARB]	# index array
+int	id		# star 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
+
+	    # Extract the appropriate field.
+	    call amovc (Memc[kip], buffer, maxch)
+	    buffer[maxch+1] = EOS
+
+	    # Decode the output value.
+	    ip = 1
+	    switch (indices[i]) {
+	    case DP_PAPID: 
+		if (ctoi (buffer, ip, id) <= 0)
+		    id = 0
+	    case DP_PAPXCEN:
+		if (ctor (buffer, ip, x) <= 0)
+		    x = INDEFR
+	    case DP_PAPYCEN:
+		if (ctor (buffer, ip, y) <= 0)
+		    y = INDEFR
+	    case DP_PAPSKY:
+	        if (ctor (buffer, ip, sky) <= 0)
+		    sky = INDEFR
+	    case DP_PAPMAG1:
+		if (ctor (buffer, ip, mag) <= 0)
+		    mag = INDEFR
+	    default:
+		call printf ("Error reading the APPHOT results.\n")
+	    }
+
+	}
+end