aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/daophot/substar
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/daophot/substar')
-rw-r--r--noao/digiphot/daophot/substar/dpgimbufr.x137
-rw-r--r--noao/digiphot/daophot/substar/dprestars.x116
-rw-r--r--noao/digiphot/daophot/substar/dpsconfirm.x18
-rw-r--r--noao/digiphot/daophot/substar/dpsubstar.x200
-rw-r--r--noao/digiphot/daophot/substar/mkpkg17
-rw-r--r--noao/digiphot/daophot/substar/t_substar.x286
6 files changed, 774 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/substar/dpgimbufr.x b/noao/digiphot/daophot/substar/dpgimbufr.x
new file mode 100644
index 00000000..d4fac359
--- /dev/null
+++ b/noao/digiphot/daophot/substar/dpgimbufr.x
@@ -0,0 +1,137 @@
+include <imhdr.h>
+
+# DP_GIMBUFR -- Maintain buffer of image lines. A new buffer is created when
+# the buffer pointer is null. No changing of buffer size is allowed, although
+# this should be added. The minimum number of image reads is used.
+
+define flush_ 91
+
+procedure dp_gimbufr (inim, outim, line1, line2, buf, flush)
+
+pointer inim # input image pointer
+pointer outim # output image pointer
+int line1 # first image line of buffer
+int line2 # last image line of buffer
+pointer buf # buffer
+bool flush # flush the current contents of the buffer
+
+int i, ncols, nlines, llast1, llast2, nllast, lp, lout
+pointer buf1, buf2
+pointer imgl2r(), impl2r()
+
+begin
+ nlines = line2 - line1 + 1
+ ncols = IM_LEN (inim, 1)
+ lp = 0
+
+ if (flush)
+ goto flush_
+
+ # If the buffer pointer is undefined then allocate memory for the
+ # buffer. If the number of columns or lines requested changes
+ # reallocate the buffer. Initialize the last line values to force
+ # a full buffer image read.
+
+ if (buf == NULL) {
+ call malloc (buf, ncols * nlines, TY_REAL)
+ #llast1 = line1 - nlines
+ #llast2 = line2 - nlines
+ llast1 = 0
+ llast2 = 0
+ } else if ((nlines > nllast)) {
+ call eprintf ("Buffer requested is larger than previous one\n")
+ return
+ }
+
+ #call printf ("line1=%d line2=%d llast1=%d llast2=%d\n")
+ #call pargi (line1)
+ #call pargi (line2)
+ #call pargi (llast1)
+ #call pargi (llast2)
+
+ # Write out the lines that are not needed any more.
+ if (line1 > llast1 && llast1 > 0) {
+ buf2 = buf
+ lout = min (llast2, line1 - 1)
+ do i = llast1, lout{
+ buf1 = impl2r (outim, i)
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ #call printf ("Writing line: %d\n")
+ #call pargi (i)
+ buf2 = buf2 + ncols
+ }
+ }
+
+ # Write out any skipped image lines.
+ if (line1 > llast2) {
+ do i = llast2 + 1, line1 - 1 {
+ buf2 = imgl2r (inim, i)
+ buf1 = impl2r (outim, i)
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ #call printf ("Copying line: %d\n")
+ #call pargi (i)
+ }
+ }
+
+ # Now move the remaining lines to the begining of the buffer.
+ if (llast2 >= line1 ) {
+ buf2 = buf + ncols * (line1 - llast1)
+ buf1 = buf
+ do i = line1, llast2 {
+ lp = lp + 1
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ #call printf ("Moving line: %d\n")
+ #call pargi (i)
+ buf2 = buf2 + ncols
+ buf1 = buf1 + ncols
+ }
+ }
+
+ # Read only the image lines with are different from the last buffer.
+ buf1 = buf + ncols * lp
+ lout = max (line1, llast2 + 1)
+ do i = lout, line2 {
+ #call printf ("Reading line: %d\n")
+ #call pargi (i)
+ buf2 = imgl2r (inim, i)
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ buf1 = buf1 + ncols
+ }
+
+ # Save the buffer parameters.
+ llast1 = line1
+ llast2 = line2
+ nllast = nlines
+
+ # Quit
+ return
+
+flush_
+ # If requested to flush the current contents of the buffer we
+ # write out lines llast1 to llast2 and then set buf == NULL.
+
+ # Flush the data buffer.
+ if (buf != NULL) {
+ buf2 = buf
+ do i = llast1, llast2 {
+ buf1 = impl2r (outim, i)
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ #call printf ("Writing line: %d\n")
+ #call pargi (i)
+ buf2 = buf2 + ncols
+ }
+ }
+
+ # Copy any remaining image lines.
+ do i = llast2 + 1, IM_LEN(inim,2) {
+ buf2 = imgl2r (inim, i)
+ buf1 = impl2r (outim, i)
+ call amovr (Memr[buf2], Memr[buf1], ncols)
+ #call printf ("Copying line: %d\n")
+ #call pargi (i)
+ }
+
+ call mfree (buf, TY_REAL)
+ buf = NULL
+ nllast = 0
+end
diff --git a/noao/digiphot/daophot/substar/dprestars.x b/noao/digiphot/daophot/substar/dprestars.x
new file mode 100644
index 00000000..3b317960
--- /dev/null
+++ b/noao/digiphot/daophot/substar/dprestars.x
@@ -0,0 +1,116 @@
+include <tbset.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+
+
+# DP_RESTARS -- Read in the IDS of the stars to be excluded, find these stars
+# in the photometry list, and set their magnitudes to INDEF.
+
+int procedure dp_restars (dao, im, ext, text_file)
+
+pointer dao # pointer to the daophot structure
+pointer im # the input image descriptor
+int ext # the exclude list file descriptor
+bool text_file # text or table file ?
+
+real tx, ty, rjunk
+pointer apsel, sp, fields, indices, key
+int i, nrow, idno, nexcl, starno
+int tbpsta(), dp_apsel(), dp_exfind()
+
+begin
+ # Get some pointers.
+ apsel = DP_APSEL(dao)
+
+ # Get some working space.
+ call smark (sp)
+ call salloc (fields, SZ_LINE, TY_CHAR)
+ call salloc (indices, 1, TY_INT)
+
+ # Initialize the read.
+ if (text_file) {
+ call pt_kyinit (key)
+ Memi[indices] = DP_PAPID
+ call dp_gappsf (Memi[indices], Memc[fields], 1)
+ } else {
+ call dp_tptinit (ext, Memi[indices])
+ nrow = tbpsta (ext, TBL_NROWS)
+ }
+
+ i = 1
+ nexcl = 0
+ repeat {
+
+ # Read the next star.
+ if (text_file) {
+ if (dp_apsel (key, ext, Memc[fields], Memi[indices], idno,
+ rjunk, rjunk, rjunk, rjunk) == EOF)
+ break
+ } else {
+ if (i > nrow)
+ break
+ call dp_tptread (ext, Memi[indices], idno, rjunk, rjunk, rjunk,
+ i)
+ }
+
+ # Subtract star from the photometry list.
+ if (idno > 0) {
+ starno = dp_exfind (Memi[DP_APID(apsel)],
+ Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+ Memr[DP_APMAG(apsel)], DP_APNUM(apsel), idno)
+ if (starno > 0) {
+ if (DP_VERBOSE(dao) == YES) {
+ call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+starno-1],
+ Memr[DP_APYCEN(apsel)+starno-1], tx, ty, 1)
+ call printf (
+ "EXCLUDING - Star:%5d X =%8.2f Y =%8.2f Mag =%8.2f\n")
+ call pargi (Memi[DP_APID(apsel)+starno-1])
+ call pargr (tx)
+ call pargr (ty)
+ call pargr (Memr[DP_APMAG(apsel)+starno-1])
+ }
+ nexcl = nexcl + 1
+ } else if (DP_VERBOSE(dao) == YES) {
+ call printf ("EXCLUDING - Star:%5d not found\n")
+ call pargi (idno)
+ }
+ }
+
+ i = i + 1
+ }
+
+ if (text_file)
+ call pt_kyfree (key)
+ call sfree (sp)
+
+ return (nexcl)
+end
+
+
+# DP_EXFIND -- Find the star to be exclude in the photometry list.
+
+int procedure dp_exfind (ids, xcen, ycen, mags, nstars, idex)
+
+int ids[ARB] # array of stellar ids
+real xcen[ARB] # array of x coordinates
+real ycen[ARB] # array of y coordinates
+real mags[ARB] # array of magnitudes
+int nstars # number of stars in photometry list
+int idex # id of star to be excluded
+
+int i, found
+
+begin
+ found = 0
+ do i = 1, nstars {
+ if (ids[i] != idex)
+ next
+ found = i
+ break
+ }
+
+ if (found > 0)
+ mags[i] = INDEFR
+
+ return (found)
+end
diff --git a/noao/digiphot/daophot/substar/dpsconfirm.x b/noao/digiphot/daophot/substar/dpsconfirm.x
new file mode 100644
index 00000000..3ed5bc27
--- /dev/null
+++ b/noao/digiphot/daophot/substar/dpsconfirm.x
@@ -0,0 +1,18 @@
+# DP_SCONFIRM -- Procedure to confirm the critical substar parameters.
+
+procedure dp_sconfirm (dao)
+
+pointer dao # pointer to the daophot structure
+
+begin
+ call printf ("\n")
+
+ # Confirm the psf radius.
+ call dp_vpsfrad (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/substar/dpsubstar.x b/noao/digiphot/daophot/substar/dpsubstar.x
new file mode 100644
index 00000000..b33358dc
--- /dev/null
+++ b/noao/digiphot/daophot/substar/dpsubstar.x
@@ -0,0 +1,200 @@
+include <mach.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+
+define EXPAND 8
+
+# DP_SUBSTAR -- Subtract the scaled and shifted PSF from the data
+
+procedure dp_substar (dao, inim, exfd, ex_text, outim)
+
+pointer dao # pointer to the DAOPHOT structure
+pointer inim # pointer to the input image
+int exfd # exclude file descriptor
+bool ex_text # text or table exclude file
+pointer outim # pointer to the output image
+
+real pradius, psfradsq, x, y, dxfrom_psf, dyfrom_psf, mag, tx, ty
+real rel_bright, maxgdata
+pointer apsel, psffit, buf, sp, index
+int i, id, line1, line2, nline_buf, x1, x2, y1, y2
+int lowy, highy, offset, nstars, ier
+int dp_restars()
+
+begin
+ # Get the daophot pointers.
+ apsel = DP_APSEL (dao)
+ psffit = DP_PSFFIT (dao)
+
+ # Exit gracefully if there are no stars.
+ #if (DP_APNUM(apsel) <= 0) {
+ #call printf ("The number of stars in the photometry list is %d\n")
+ #call pargi (DP_APNUM(apsel))
+ #return
+ #}
+
+ # Check for stars to be excluded.
+ if (exfd != NULL) {
+ if (dp_restars (dao, inim, exfd, ex_text) <= 0)
+ ;
+ }
+
+ # Compute the size of subraster to read from the PSF image.
+ if (DP_PSFSIZE(dao) == 0)
+ pradius = DP_PSFRAD(dao)
+ else
+ pradius = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0
+ psfradsq = pradius * pradius
+
+ # Set the maximum good bad limit.
+ if (IS_INDEFR (DP_MAXGDATA(dao)))
+ maxgdata = MAX_REAL
+ else
+ maxgdata = DP_MAXGDATA(dao)
+
+ # Get some working memory.
+ call smark (sp)
+ call salloc (index, DP_APNUM (apsel), TY_INT)
+
+ # Sort the photometry on increasing Y.
+ if (DP_APNUM(apsel) > 0)
+ call quick (Memr[DP_APYCEN(apsel)], DP_APNUM(apsel), Memi[index],
+ ier)
+
+ # Initialize the boundary of the buffer.
+ buf = NULL
+ line1 = 0
+ line2 = 0
+ nline_buf = EXPAND * pradius
+
+ nstars = 0
+ do i = 1, DP_APNUM (apsel) {
+
+ # Get the data for the next star.
+ offset = Memi[index+i-1] - 1
+ x = Memr[DP_APXCEN(apsel)+offset]
+ y = Memr[DP_APYCEN(apsel)+i-1]
+ id = Memi[DP_APID(apsel)+offset]
+ mag = Memr[DP_APMAG (apsel)+offset]
+ call dp_wpsf (dao, inim, x, y, dxfrom_psf, dyfrom_psf, 1)
+ dxfrom_psf = (dxfrom_psf - 1.0) / DP_PSFX(psffit) - 1.0
+ dyfrom_psf = (dyfrom_psf - 1.0) / DP_PSFY(psffit) - 1.0
+
+ # Reject star is the magnitude is INDEF.
+ if (IS_INDEFR(x) || IS_INDEFR(y) || IS_INDEFR(mag)) {
+ if (DP_VERBOSE(dao) == YES) {
+ if (IS_INDEFR(x) || IS_INDEFR(y)) {
+ tx = x
+ ty = y
+ } else
+ call dp_wout (dao, inim, x, y, tx, ty, 1)
+ call printf (
+ "REJECTING - Star:%5d X =%8.2f Y =%8.2f Mag =%8.2f\n")
+ call pargi (id)
+ call pargr (tx)
+ call pargr (ty)
+ call pargr (mag)
+ }
+ next
+ }
+
+ # Print out the verbose message.
+ if (DP_VERBOSE(dao) == YES) {
+ call dp_wout (dao, inim, x, y, tx, ty, 1)
+ call printf (
+ "SUBTRACTING - Star:%5d X =%8.2f Y =%8.2f Mag =%8.2f\n")
+ call pargi (id)
+ call pargr (tx)
+ call pargr (ty)
+ call pargr (mag)
+ }
+
+ # Determine the range of lines required.
+ lowy = max (1, int (y - pradius) + 1)
+ highy = min (IM_LEN (inim, 2), int (y + pradius))
+ if (highy > line2) {
+ line1 = max (1, lowy)
+ line2 = min (line1 + nline_buf, IM_LEN (inim, 2))
+ call dp_gimbufr (inim, outim, line1, line2, buf, false)
+ }
+
+ # Change coordinates to reference frame of buffer.
+ y = y - line1 + 1.0
+ y1 = max (1, int (y - pradius) + 1)
+ y2 = min (line2 - line1 + 1, int (y + pradius))
+ x1 = max (1, int (x - pradius) + 1)
+ x2 = min (IM_LEN (inim, 1), int (x + pradius))
+
+ # Computee the relative brightness.
+ rel_bright = DAO_RELBRIGHT (psffit, mag)
+
+ # Subtract this star.
+ call dp_sstar (dao, Memr[buf], int (IM_LEN(inim,1)), nline_buf,
+ x1, x2, y1, y2, x, y, psfradsq, rel_bright, dxfrom_psf,
+ dyfrom_psf, maxgdata)
+
+ nstars = nstars + 1
+ }
+
+ # Flush the remaining lines in the image buffer.
+ call dp_gimbufr (inim, outim, y1, y2, buf, true)
+
+ # Summarize data on the number of stars subtracted.
+ if (DP_VERBOSE(dao) == YES) {
+ call printf (
+ "\nA total of %d stars were subtracted out of a possible %d\n")
+ call pargi (nstars)
+ call pargi (DP_APNUM(apsel))
+ }
+
+ # Free memory.
+ call sfree (sp)
+end
+
+
+# DP_SSTAR -- Subtract the star from the image.
+
+procedure dp_sstar (dao, data, nx, ny, x1, x2, y1, y2, xstar, ystar, psfradsq,
+ rel_bright, dxfrom_psf, dyfrom_psf, maxgdata)
+
+pointer dao # pointer to the daophot structure
+real data[nx,ny] # sata buffer
+int nx, ny # size of buffer
+int x1, x2, y1, y2 # area of interest
+real xstar, ystar # position of star to subtract
+real psfradsq # PSF radius ** 2
+real rel_bright # relative brightness of star
+real dxfrom_psf, dyfrom_psf # not currently used
+real maxgdata # maximum good data
+
+int ix, iy
+pointer psffit
+real dx, dy, dxsq, dysq, radsq, dvdx, dvdy
+real dp_usepsf()
+
+begin
+ psffit = DP_PSFFIT(dao)
+ do iy = y1, y2 {
+ dy = real (iy) - ystar
+ dysq = dy * dy
+ do ix = x1, x2 {
+ if (data[ix,iy] > maxgdata)
+ next
+ dx = real (ix) - xstar
+ dxsq = dx * dx
+ radsq = dxsq + dysq
+ if (radsq >= psfradsq) {
+ if (dx > 0.0)
+ break
+ next
+ }
+ data[ix,iy] = data[ix,iy] - rel_bright *
+ 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),
+ dxfrom_psf, dyfrom_psf, dvdx, dvdy)
+ }
+ }
+end
diff --git a/noao/digiphot/daophot/substar/mkpkg b/noao/digiphot/daophot/substar/mkpkg
new file mode 100644
index 00000000..d51cf73e
--- /dev/null
+++ b/noao/digiphot/daophot/substar/mkpkg
@@ -0,0 +1,17 @@
+# SUBSTAR task
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+ dpsconfirm.x
+ dpsubstar.x <mach.h> <imhdr.h> \
+ ../lib/daophotdef.h ../lib/apseldef.h
+ dpgimbufr.x <imhdr.h>
+ dprestars.x <tbset.h> ../lib/daophotdef.h \
+ ../lib/apseldef.h
+ t_substar.x <fset.h> <imhdr.h> \
+ ../lib/daophotdef.h
+ ;
diff --git a/noao/digiphot/daophot/substar/t_substar.x b/noao/digiphot/daophot/substar/t_substar.x
new file mode 100644
index 00000000..11ffc72b
--- /dev/null
+++ b/noao/digiphot/daophot/substar/t_substar.x
@@ -0,0 +1,286 @@
+include <fset.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+
+# T_SUBSTAR -- Procedure to subtract DAOPHOT photometry from an image.
+
+procedure t_substar ()
+
+pointer image # name of the image
+pointer photfile # input photometry file
+pointer exfile # input exclude file
+pointer psfimage # name of the output PSF
+pointer subimage # subtracted image
+
+pointer sp, input, output, dao, outfname, str
+int psffd, photfd, root, verify, update, wcs
+int imlist, limlist, alist, lalist, pimlist, lpimlist, simlist, lsimlist
+int exfd, elist, lelist, cache, req_size, old_size, buf_size, memstat
+bool ap_text, ex_text
+
+pointer immap(), tbtopn()
+int open(), fnldir(), strlen(), strncmp(), access(), fstati(), btoi()
+int imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb(), fntgfnb()
+int 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)
+
+ # Get some working memory.
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (photfile, SZ_FNAME, TY_CHAR)
+ call salloc (exfile, SZ_FNAME, TY_CHAR)
+ call salloc (psfimage, SZ_FNAME, TY_CHAR)
+ call salloc (subimage, SZ_FNAME, TY_CHAR)
+ call salloc (outfname, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Get the various task parameters.
+ call clgstr ("image", Memc[image], SZ_FNAME)
+ call clgstr ("photfile", Memc[photfile], SZ_FNAME)
+ call clgstr ("exfile", Memc[exfile], SZ_FNAME)
+ call clgstr ("psfimage", Memc[psfimage], SZ_FNAME)
+ call clgstr ("subimage", Memc[subimage], SZ_FNAME)
+ 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[photfile], NO)
+ lalist = fntlenb (alist)
+ elist = fntopnb (Memc[exfile], NO)
+ lelist = fntlenb (elist)
+ pimlist = imtopen (Memc[psfimage])
+ lpimlist = imtlen (pimlist)
+ simlist = imtopen (Memc[subimage])
+ lsimlist = imtlen (simlist)
+
+ # Test that the lengths of the photometry file, psf image and
+ # subtracted image lists are the same as the length of the input
+ # image list.
+
+ if ((limlist != lalist) && (strncmp (Memc[photfile], DEF_DEFNAME,
+ DEF_LENDEFNAME) != 0)) {
+ call imtclose (imlist)
+ call fntclsb (alist)
+ call fntclsb (elist)
+ call imtclose (pimlist)
+ call imtclose (simlist)
+ call sfree (sp)
+ call error (0,
+ "Incompatible image and photometry file list lengths")
+ }
+
+ if ((lelist != 0) && (limlist != lelist) && (strncmp (Memc[exfile],
+ DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+ call imtclose (imlist)
+ call fntclsb (alist)
+ call fntclsb (elist)
+ call imtclose (pimlist)
+ call imtclose (simlist)
+ call sfree (sp)
+ call error (0,
+ "Incompatible image and exclude file list lengths")
+ }
+
+ if ((limlist != lpimlist) && (strncmp (Memc[psfimage], DEF_DEFNAME,
+ DEF_LENDEFNAME) != 0)) {
+ call imtclose (imlist)
+ call fntclsb (alist)
+ call fntclsb (elist)
+ call imtclose (pimlist)
+ call imtclose (simlist)
+ call sfree (sp)
+ call error (0,
+ "Incompatible image and psf file list lengths")
+ }
+
+ if ((limlist != lsimlist) && (strncmp (Memc[subimage], DEF_DEFNAME,
+ DEF_LENDEFNAME) != 0)) {
+ call imtclose (imlist)
+ call fntclsb (alist)
+ call fntclsb (elist)
+ call imtclose (pimlist)
+ call imtclose (simlist)
+ call sfree (sp)
+ call error (0,
+ "Incompatible image and subtracted image list lengths")
+ }
+
+ # Initialize the DAOPHOT structure and get the pset parameters.
+ call dp_gppars (dao)
+ call dp_seti (dao, VERBOSE, btoi (clgetb ("verbose")))
+
+ # Verify the critical parameters.
+ if (verify == YES) {
+ call dp_sconfirm (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 PSF structure.
+ call dp_fitsetup (dao)
+
+ # Initialize the star list.
+ call dp_apselsetup (dao)
+
+ # Loop over the images
+ while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+ # Open input and output images
+ input = immap (Memc[image], READ_ONLY, 0)
+ call dp_sets (dao, INIMAGE, Memc[image])
+
+ # Cache the input image pixels.
+ req_size = MEMFUDGE * (2 * IM_LEN(input,1) * IM_LEN(input,2) *
+ sizeof (IM_PIXTYPE(input)))
+ memstat = dp_memstat (cache, req_size, old_size)
+ if (memstat == YES)
+ call dp_pcache (input, INDEFI, buf_size)
+
+ # If the output image name is DEF_DEFNAME, dir$default or a
+ # directory specification then the extension "sub" is added to
+ # the image name and a suitable version number is appended to the
+ # output name.
+
+ if (imtgetim (simlist, Memc[subimage], SZ_FNAME) == EOF)
+ call strcpy (DEF_DEFNAME, Memc[subimage], SZ_FNAME)
+ root = fnldir (Memc[subimage], Memc[outfname], SZ_FNAME)
+ if (strncmp (DEF_DEFNAME, Memc[subimage + root],
+ DEF_LENDEFNAME) == 0 || root == strlen (Memc[subimage])) {
+ call dp_oimname (Memc[image], Memc[outfname], "sub",
+ Memc[outfname], SZ_FNAME)
+ output = immap (Memc[outfname], NEW_COPY, input)
+ } else {
+ call strcpy (Memc[subimage], Memc[outfname], SZ_FNAME)
+ output = immap (Memc[outfname], NEW_COPY, input)
+ }
+ call dp_sets (dao, OUTIMAGE, Memc[outfname])
+ if (memstat == YES)
+ call dp_pcache (output, INDEFI, buf_size)
+
+ # Open input photometry table and read in the photometry.
+ if (fntgfnb (alist, Memc[photfile], SZ_FNAME) == EOF)
+ call strcpy (DEF_DEFNAME, Memc[photfile], SZ_FNAME)
+ root = fnldir (Memc[photfile], Memc[outfname], SZ_FNAME)
+ if (strncmp (DEF_DEFNAME, Memc[photfile+root],
+ DEF_LENDEFNAME) == 0 || root == strlen (Memc[photfile]))
+ call dp_inname (Memc[image], Memc[outfname], "nst",
+ Memc[outfname], SZ_FNAME)
+ else
+ call strcpy (Memc[photfile], Memc[outfname], SZ_FNAME)
+ ap_text = itob (access (Memc[outfname], 0, TEXT_FILE))
+ if (ap_text)
+ photfd = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+ else
+ photfd = tbtopn (Memc[outfname], READ_ONLY, 0)
+ call dp_wgetapert (dao, input, photfd, DP_MAXNSTAR(dao), ap_text)
+ call dp_sets (dao, INPHOTFILE, Memc[outfname])
+
+ # Open the input exclude file.
+ if (lelist == 0) {
+ exfd = NULL
+ Memc[outfname] = EOS
+ } else {
+ if (fntgfnb (elist, Memc[exfile], SZ_FNAME) == EOF)
+ call strcpy (DEF_DEFNAME, Memc[exfile], SZ_FNAME)
+ root = fnldir (Memc[exfile], Memc[outfname], SZ_FNAME)
+ if (strncmp (DEF_DEFNAME, Memc[exfile+root],
+ DEF_LENDEFNAME) == 0 || root == strlen (Memc[exfile]))
+ call dp_inname (Memc[image], Memc[outfname], "pst",
+ Memc[outfname], SZ_FNAME)
+ else
+ call strcpy (Memc[exfile], Memc[outfname], SZ_FNAME)
+ ex_text = itob (access (Memc[outfname], 0, TEXT_FILE))
+ if (ex_text)
+ exfd = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+ else
+ exfd = tbtopn (Memc[outfname], READ_ONLY, 0)
+ }
+ call dp_sets (dao, COORDS, Memc[outfname])
+
+ # Read in the PSF
+ 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)
+ psffd = immap (Memc[outfname], READ_ONLY, 0)
+ call dp_readpsf (dao, psffd)
+ call dp_sets (dao, PSFIMAGE, Memc[outfname])
+
+ # Now go and subtract those stars!
+ call dp_substar (dao, input, exfd, ex_text, output)
+
+ # Close the input and output images.
+ call imunmap (input)
+ call imunmap (output)
+
+ # Close the photometry file.
+ if (ap_text)
+ call close (photfd)
+ else
+ call tbtclo (photfd)
+
+ # Close the exclude file.
+ if (ex_text)
+ call close (exfd)
+ else
+ call tbtclo (exfd)
+
+ # Close the PSF image.
+ call imunmap (psffd)
+
+ # Uncache memory
+ call fixmem (old_size)
+
+ }
+
+ # Close the lists.
+ call imtclose (imlist)
+ call fntclsb (alist)
+ call fntclsb (elist)
+ call imtclose (pimlist)
+ call imtclose (simlist)
+
+ # Free the daophot structures.
+ call dp_apclose (dao)
+ call dp_fitclose (dao)
+ call dp_free (dao)
+
+ call sfree (sp)
+end