aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/apphot/aputil
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/apphot/aputil')
-rw-r--r--noao/digiphot/apphot/aputil/apbsmooth.x33
-rw-r--r--noao/digiphot/apphot/aputil/apclip.x23
-rw-r--r--noao/digiphot/apphot/aputil/apdate.x29
-rw-r--r--noao/digiphot/apphot/aputil/apdiverr.x7
-rw-r--r--noao/digiphot/apphot/aputil/apfnames.x280
-rw-r--r--noao/digiphot/apphot/aputil/apgscur.x88
-rw-r--r--noao/digiphot/apphot/aputil/apgtools.x141
-rw-r--r--noao/digiphot/apphot/aputil/aplucy.x47
-rw-r--r--noao/digiphot/apphot/aputil/apmapr.x21
-rw-r--r--noao/digiphot/apphot/aputil/apmeds.x154
-rw-r--r--noao/digiphot/apphot/aputil/apmoments.x138
-rw-r--r--noao/digiphot/apphot/aputil/apnlfuncs.x375
-rw-r--r--noao/digiphot/apphot/aputil/appcache.x83
-rw-r--r--noao/digiphot/apphot/aputil/aprmwhite.x22
-rw-r--r--noao/digiphot/apphot/aputil/aptopt.x232
-rw-r--r--noao/digiphot/apphot/aputil/apvectors.x515
-rw-r--r--noao/digiphot/apphot/aputil/mkpkg25
17 files changed, 2213 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/aputil/apbsmooth.x b/noao/digiphot/apphot/aputil/apbsmooth.x
new file mode 100644
index 00000000..28c53c00
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apbsmooth.x
@@ -0,0 +1,33 @@
+# AP_BSMOOTH - Box car smooth a histogram 1, 2 or 3 times.
+
+procedure ap_bsmooth (hgm, shgm, nbins, nker, iter)
+
+real hgm[ARB] # the original histogram
+real shgm[ARB] # the smoothed histogram
+int nbins # length of the histogram
+int nker # half width of box kernel
+int iter # number of iterations
+
+pointer sp, work1, work2
+
+begin
+ # Smooth the histogram
+ switch (iter) {
+ case 1:
+ call ap_sboxr (hgm, shgm, nbins, nker)
+ case 2:
+ call smark (sp)
+ call salloc (work1, nbins, TY_REAL)
+ call ap_sboxr (hgm, Memr[work1], nbins, nker)
+ call ap_sboxr (Memr[work1], shgm, nbins, nker)
+ call sfree (sp)
+ default:
+ call smark (sp)
+ call salloc (work1, nbins, TY_REAL)
+ call salloc (work2, nbins, TY_REAL)
+ call ap_sboxr (hgm, Memr[work1], nbins, nker)
+ call ap_sboxr (Memr[work1], Memr[work2], nbins, nker)
+ call ap_sboxr (Memr[work2], shgm, nbins, nker)
+ call sfree (sp)
+ }
+end
diff --git a/noao/digiphot/apphot/aputil/apclip.x b/noao/digiphot/apphot/aputil/apclip.x
new file mode 100644
index 00000000..62a52c30
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apclip.x
@@ -0,0 +1,23 @@
+# AP_CLIP -- Clip the ends of a sorted pixel distribution by a certain
+# percent.
+
+int procedure ap_clip (skypix, index, npix, loclip, hiclip, loindex, hiindex)
+
+real skypix[ARB] # the unsorted array of sky pixels
+int index[ARB] # the array of sorted indices
+int npix # the number of sky pixels
+real loclip, hiclip # the clipping factors in percent
+int loindex, hiindex # the clipping indices
+
+begin
+ # Sort the pixels.
+ call apqsort (skypix, index, index, npix)
+
+ # Determine the clipping factors.
+ loindex = nint (0.01 * loclip * npix) + 1
+ hiindex = npix - nint (0.01 * hiclip * npix)
+ if ((hiindex - loindex + 1) <= 0)
+ return (npix)
+ else
+ return (loindex - 1 + npix - hiindex)
+end
diff --git a/noao/digiphot/apphot/aputil/apdate.x b/noao/digiphot/apphot/aputil/apdate.x
new file mode 100644
index 00000000..21faea88
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apdate.x
@@ -0,0 +1,29 @@
+include <time.h>
+
+# APDATE -- Fetch the date and time strings required for the apphot output
+# files.
+
+procedure apdate (date, time, maxch)
+
+char date[ARB] # the date string
+char time[ARB] # the time string
+int maxch # the maximum number of character in the string
+
+int tm[LEN_TMSTRUCT]
+long ctime
+long clktime()
+#long lsttogmt()
+
+begin
+ ctime = clktime (long(0))
+ #ctime = lsttogmt (ctime)
+ call brktime (ctime, tm)
+ call sprintf (date, maxch, "%04d-%02d-%02d")
+ call pargi (TM_YEAR(tm))
+ call pargi (TM_MONTH(tm))
+ call pargi (TM_MDAY(tm))
+ call sprintf (time, maxch, "%02d:%02d:%02d")
+ call pargi (TM_HOUR(tm))
+ call pargi (TM_MIN(tm))
+ call pargi (TM_SEC(tm))
+end
diff --git a/noao/digiphot/apphot/aputil/apdiverr.x b/noao/digiphot/apphot/aputil/apdiverr.x
new file mode 100644
index 00000000..d3c36444
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apdiverr.x
@@ -0,0 +1,7 @@
+# AP_DIVERR -- Error condition procedure which returns a value of unity.
+
+real procedure ap_diverr ()
+
+begin
+ return (1.0)
+end
diff --git a/noao/digiphot/apphot/aputil/apfnames.x b/noao/digiphot/apphot/aputil/apfnames.x
new file mode 100644
index 00000000..95e56167
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apfnames.x
@@ -0,0 +1,280 @@
+
+# APINNAME -- Construct an apphot input file name.
+# If input is null or a directory, a name is constructed from the root
+# of the image name and the extension. The disk is searched to avoid
+# name collisions.
+
+procedure apinname (image, input, ext, name, maxch)
+
+char image[ARB] # image name
+char input[ARB] # input directory or name
+char ext[ARB] # extension
+char name[ARB] # input name
+int maxch # maximum size of name
+
+int ndir, nimdir, clindex, clsize
+pointer sp, root, str
+int fnldir(), strlen()
+
+begin
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ ndir = fnldir (input, name, maxch)
+ if (strlen (input) == ndir) {
+ call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+ Memc[str], SZ_FNAME, clindex, clsize)
+ nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+ if (clindex >= 0) {
+ call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ call pargstr (ext)
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargstr (ext)
+ }
+ call apiversion (name, name, maxch)
+ } else
+ call strcpy (input, name, maxch)
+
+ call sfree (sp)
+end
+
+
+# APOUTNAME -- Construct an apphot output file name.
+# If output is null or a directory, a name is constructed from the root
+# of the image name and the extension. The disk is searched to avoid
+# name collisions.
+
+procedure apoutname (image, output, ext, name, maxch)
+
+char image[ARB] # image name
+char output[ARB] # output directory or name
+char ext[ARB] # extension
+char name[ARB] # output name
+int maxch # maximum size of name
+
+int ndir, nimdir, clindex, clsize
+pointer sp, root, str
+int fnldir(), strlen()
+
+begin
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ ndir = fnldir (output, name, maxch)
+ if (strlen (output) == ndir) {
+ call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+ Memc[str], SZ_FNAME, clindex, clsize)
+ nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+ if (clindex >= 0) {
+ call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ call pargstr (ext)
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargstr (ext)
+ }
+ call apoversion (name, name, maxch)
+ } else
+ call strcpy (output, name, maxch)
+
+ call sfree (sp)
+end
+
+
+# APTMPIMAGE -- Generate a temporary image name either by calling a system
+# routine or by appending the image name to a user specified prefix.
+
+int procedure aptmpimage (image, prefix, tmp, name, maxch)
+
+char image[ARB] # image name
+char prefix[ARB] # user supplied prefix
+char tmp[ARB] # user supplied temporary root
+char name[ARB] # output name
+int maxch # max number of chars
+
+int npref, ndir, nimdir, clindex, clsize
+pointer sp, root, str
+int fnldir(), strlen()
+
+begin
+ npref = strlen (prefix)
+ ndir = fnldir (prefix, name, maxch)
+ if (npref == ndir) {
+ call mktemp (tmp, name[ndir+1], maxch)
+ return (NO)
+ } else {
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call strcpy (prefix, name, npref)
+ call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+ Memc[str], SZ_FNAME, clindex, clsize)
+ nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+ if (clindex >= 0) {
+ call sprintf (name[npref+1], maxch, "%s%d")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ } else {
+ call sprintf (name[npref+1], maxch, "%s")
+ call pargstr (Memc[root+nimdir])
+ }
+ call sfree (sp)
+ return (YES)
+ }
+end
+
+
+# APIMROOT -- Fetch the root image name minus the directory specification
+# and the section notation.
+
+procedure apimroot (image, root, maxch)
+
+char image[ARB] # image specification
+char root[ARB] # output root name
+int maxch # maximum number of characters
+
+pointer sp, imroot, kernel, section, str
+int clindex, clsize, nchars
+int fnldir()
+
+begin
+ call smark (sp)
+ call salloc (imroot, SZ_PATHNAME, TY_CHAR)
+ call salloc (kernel, SZ_FNAME, TY_CHAR)
+ call salloc (section, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_PATHNAME, TY_CHAR)
+
+ call imparse (image, Memc[imroot], SZ_PATHNAME, Memc[kernel], SZ_FNAME,
+ Memc[section], SZ_FNAME, clindex, clsize)
+ nchars = fnldir (Memc[imroot], Memc[str], SZ_PATHNAME)
+ if (clindex >= 0) {
+ call sprintf (root, maxch, "%s[%d]%s%s")
+ call pargstr (Memc[imroot+nchars])
+ call pargi (clindex)
+ call pargstr (Memc[kernel])
+ call pargstr (Memc[section])
+ } else {
+ call sprintf (root, maxch, "%s%s%s")
+ call pargstr (Memc[imroot+nchars])
+ call pargstr (Memc[kernel])
+ call pargstr (Memc[section])
+ }
+
+ call sfree (sp)
+end
+
+
+
+# APFROOT -- Fetch the file name minus the directory specification,
+
+procedure apfroot (filename, root, maxch)
+
+char filename[ARB] # input file name
+char root[ARB] # output root file name
+int maxch # maximum number of characters
+
+pointer sp, str
+int nchars
+int fnldir()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_PATHNAME, TY_CHAR)
+
+ nchars = fnldir (filename, Memc[str], SZ_PATHNAME)
+ call strcpy (filename[nchars+1], root, maxch)
+
+ call sfree (sp)
+end
+
+
+# APOVERSION -- Compute the next available version number of a given file
+# name template and output the new file name.
+
+procedure apoversion (template, filename, maxch)
+
+char template[ARB] # name template
+char filename[ARB] # output name
+int maxch # maximum number of characters
+
+char period
+int newversion, version, len
+pointer sp, list, name
+int fntgfnb() strldx(), ctoi(), fntopnb()
+
+begin
+ # Allocate temporary space
+ call smark (sp)
+ call salloc (name, maxch, TY_CHAR)
+ period = '.'
+ list = fntopnb (template, NO)
+
+ # Loop over the names in the list searchng for the highest version.
+ newversion = 0
+ while (fntgfnb (list, Memc[name], maxch) != EOF) {
+ len = strldx (period, Memc[name])
+ len = len + 1
+ if (ctoi (Memc[name], len, version) <= 0)
+ next
+ newversion = max (newversion, version)
+ }
+
+ # Make new output file name.
+ len = strldx (period, template)
+ call strcpy (template, filename, len)
+ call sprintf (filename[len+1], maxch, "%d")
+ call pargi (newversion + 1)
+
+ call fntclsb (list)
+ call sfree (sp)
+end
+
+
+# APIVERSION -- Compute the highest available version number of a given file
+# name template and output the file name.
+
+procedure apiversion (template, filename, maxch)
+
+char template[ARB] # name template
+char filename[ARB] # output name
+int maxch # maximum number of characters
+
+char period
+int newversion, version, len
+pointer sp, list, name
+int fntgfnb() strldx(), ctoi(), fntopnb()
+
+begin
+ # Allocate temporary space
+ call smark (sp)
+ call salloc (name, maxch, TY_CHAR)
+ period = '.'
+ list = fntopnb (template, NO)
+
+ # Loop over the names in the list searchng for the highest version.
+ newversion = 0
+ while (fntgfnb (list, Memc[name], maxch) != EOF) {
+ len = strldx (period, Memc[name])
+ len = len + 1
+ if (ctoi (Memc[name], len, version) <= 0)
+ next
+ newversion = max (newversion, version)
+ }
+
+ # Make new output file name.
+ len = strldx (period, template)
+ call strcpy (template, filename, len)
+ call sprintf (filename[len+1], maxch, "%d")
+ call pargi (newversion)
+
+ call fntclsb (list)
+ call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/apgscur.x b/noao/digiphot/apphot/aputil/apgscur.x
new file mode 100644
index 00000000..5cc8a6df
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apgscur.x
@@ -0,0 +1,88 @@
+include <gset.h>
+include <fset.h>
+
+# APGSCUR -- Read the x and y coordinates of an object from a file and move
+# the cursor to those coordinates.
+
+int procedure apgscur (sl, gd, xcur, ycur, prev_num, req_num, num)
+
+int sl # coordinate file descriptor
+pointer gd # pointer to graphics stream
+real xcur, ycur # x cur and y cur
+int prev_num # previous number
+int req_num # requested number
+int num # list number
+
+int stdin, nskip, ncount
+pointer sp, fname
+int fscan(), nscan(), strncmp()
+errchk greactivate, gdeactivate, gscur
+
+begin
+ if (sl == NULL)
+ return (EOF)
+
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+
+ # Find the number of objects to be skipped.
+ call fstats (sl, F_FILENAME, Memc[fname], SZ_FNAME)
+ if (strncmp ("STDIN", Memc[fname], 5) == 0) {
+ stdin = YES
+ nskip = 1
+ } else {
+ stdin = NO
+ if (req_num <= prev_num) {
+ call seek (sl, BOF)
+ nskip = req_num
+ } else
+ nskip = req_num - prev_num
+ }
+
+ # Initialize.
+ ncount = 0
+ num = prev_num
+
+ # Read the coordinates and write the cursor.
+ repeat {
+
+ # Print the prompt if file is STDIN.
+ if (stdin == YES) {
+ call printf ("Type object x and y coordinates: ")
+ call flush (STDOUT)
+ }
+
+ # Fetch the coordinates.
+ if (fscan (sl) != EOF) {
+ call gargr (xcur)
+ call gargr (ycur)
+ if (nscan () == 2) {
+ ncount = ncount + 1
+ num = num + 1
+ }
+ } else
+ ncount = EOF
+
+ # Move the cursor.
+ if (gd != NULL && (ncount == nskip || ncount == EOF)) {
+ iferr {
+ call greactivate (gd, 0)
+ call gscur (gd, xcur, ycur)
+ call gdeactivate (gd, 0)
+ } then
+ ;
+ }
+
+ } until (ncount == EOF || ncount == nskip)
+
+ call sfree (sp)
+
+ if (ncount == EOF) {
+ return (EOF)
+ } else if (nskip == req_num) {
+ num = ncount
+ return (ncount)
+ } else {
+ return (num)
+ }
+end
diff --git a/noao/digiphot/apphot/aputil/apgtools.x b/noao/digiphot/apphot/aputil/apgtools.x
new file mode 100644
index 00000000..e35cb543
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apgtools.x
@@ -0,0 +1,141 @@
+include <gset.h>
+include <pkg/gtools.h>
+
+# AP_GTINIT -- Initialize the gtools package for the apphot routines.
+
+pointer procedure ap_gtinit (image, wx, wy)
+
+char image[ARB] # the image name
+real wx, wy # center of sky subraster
+
+pointer sp, gt, str
+pointer gt_init()
+
+begin
+ # Allocate working space.
+ gt = gt_init ()
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Set the plot title.
+ call sprintf (Memc[str], SZ_LINE, "Image: %s: %.2f %.2f\n")
+ call pargstr (image)
+ call pargr (wx)
+ call pargr (wy)
+ call gt_sets (gt, GTTITLE, Memc[str])
+
+ call sfree (sp)
+ return (gt)
+end
+
+
+# AP_GTFREE -- Free the gtools package.
+
+procedure ap_gtfree (gt)
+
+pointer gt # pointer to gtools structure
+
+begin
+ call gt_free (gt)
+end
+
+
+# AP_PLOTRAD -- Plot the radial profile of a list of pixels.
+
+procedure ap_plotrad (gd, gt, r, i, npts, polymark)
+
+pointer gd # pointer to graphics stream
+pointer gt # the GTOOLS pointer
+real r[ARB] # the radii array
+real i[ARB] # the intensity array
+int npts # number of points
+char polymark[ARB] # polyline type
+
+begin
+ call gt_sets (gt, GTTYPE, "mark")
+ call gt_sets (gt, GTMARK, polymark)
+ call gt_plot (gd, gt, r, i, npts)
+end
+
+
+# AP_PLOTPTS -- Plot the radial profile of a list of pixels excluding points
+# that are outside the plotting window altogether.
+
+procedure ap_plotpts (gd, gt, r, i, npts, xmin, xmax, ymin, ymax, polymark)
+
+pointer gd # pointer to graphics stream
+pointer gt # the GTOOLS pointer
+real r[ARB] # the radii array
+real i[ARB] # the intensity array
+int npts # number of points
+real xmin, xmax # the x plot limits
+real ymin, ymax # the y plot limits
+char polymark[ARB] # polyline type
+
+int j
+
+begin
+ call gt_sets (gt, GTTYPE, "mark")
+ call gt_sets (gt, GTMARK, polymark)
+ do j = 1, npts {
+ if (r[j] < xmin || r[j] > xmax)
+ next
+ if (i[j] < ymin || i[j] > ymax)
+ next
+ call gt_plot (gd, gt, r[j], i[j], 1)
+
+ }
+end
+
+
+# AP_RSET -- Set up the parameters for the radial profile plot.
+
+procedure ap_rset (gd, gt, xmin, xmax, ymin, ymax, xscale)
+
+pointer gd # pointer to GRAPHICS stream
+pointer gt # pointer to GTOOLS structure
+real xmin, xmax # min and max of x vector
+real ymin, ymax # min and max of y vector
+real xscale # image scale
+
+pointer sp, str, title
+real vx1, vx2, vy1, vy2, aspect
+real gstatr()
+
+begin
+ call smark (sp)
+ call salloc (title, 3 * SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Reset the aspect ratio.
+ aspect = gstatr (gd, G_ASPECT)
+ call gsetr (gd, G_ASPECT, 0.75)
+
+ # Construct the title.
+ call sysid (Memc[title], 3 * SZ_LINE)
+ call strcat ("\n", Memc[title], 3 * SZ_LINE)
+ call gt_gets (gt, GTTITLE, Memc[str], SZ_LINE)
+ call strcat (Memc[str], Memc[title], 3 * SZ_LINE)
+ call strcat ("\n\n", Memc[title], 3 * SZ_LINE)
+
+ # Draw three axes.
+ call gseti (gd, G_XDRAWAXES, 2)
+ call gswind (gd, xmin / xscale, xmax / xscale, ymin, ymax)
+ call glabax (gd, Memc[title], "", "Intensity")
+
+ # Draw the bottom x axis
+ call gseti (gd, G_YDRAWAXES, 0)
+ call gseti (gd, G_XDRAWAXES, 1)
+ call ggview (gd, vx1, vx2, vy1, vy2)
+ call gsview (gd, vx1, vx2, vy1, vy2)
+ call gswind (gd, xmin, xmax, ymin, ymax)
+ call glabax (gd, "",
+ "Radial Distance (lower-pixels, upper-scale units)", "")
+
+ # Reset to standard gio parameters.
+ call gseti (gd, G_YDRAWAXES, 3)
+ call gseti (gd, G_XDRAWAXES, 3)
+ call gsetr (gd, G_ASPECT, aspect)
+
+ call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/aplucy.x b/noao/digiphot/apphot/aputil/aplucy.x
new file mode 100644
index 00000000..f70a39ed
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aplucy.x
@@ -0,0 +1,47 @@
+# AP_LUCY_SMOOTH - Lucy smooth a histogram with a box shaped function.
+
+procedure ap_lucy_smooth (hgm, shgm, nbins, nker, iter)
+
+real hgm[ARB] # the original histogram
+real shgm[ARB] # smoothed histogram
+int nbins # length of the histogram
+int nker # length of box kernel
+int iter # number of iterations
+
+int i, j, nsmooth
+pointer sp, work1, work2
+real ap_diverr()
+
+begin
+ # Allocate space and clear the work arrays.
+ call smark (sp)
+ call salloc (work1, nbins, TY_REAL)
+ call salloc (work2, nbins, TY_REAL)
+ call aclrr (Memr[work1], nbins)
+ call aclrr (Memr[work2], nbins)
+ call aclrr (shgm, nbins)
+
+ # Compute the smoothing length and initialize output array.
+ nsmooth = nbins - nker + 1
+ #call ap_aboxr (hgm, shgm[1+nker/2], nsmooth, nker)
+ call ap_sboxr (hgm, shgm, nbins, nker)
+
+ # Iterate on the solution.
+ do i = 1, iter {
+ #call ap_aboxr (shgm, Memr[work1+nker/2], nsmooth, nker)
+ call ap_sboxr (shgm, Memr[work1], nbins, nker)
+ #do j = 1 + nker / 2, nsmooth + nker / 2 {
+ do j = 1, nbins {
+ if (Memr[work1+j-1] == 0.0)
+ Memr[work1+j-1] = ap_diverr ()
+ else
+ #Memr[work1+j-1] = hgm[j] / Memr[work1+j-1]
+ Memr[work1+j-1] = shgm[j] / Memr[work1+j-1]
+ }
+ #call ap_aboxr (Memr[work1], Memr[work2+nker/2], nsmooth, nker)
+ call ap_sboxr (Memr[work1], Memr[work2], nbins, nker)
+ call amulr (shgm, Memr[work2], shgm, nbins)
+ }
+
+ call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/apmapr.x b/noao/digiphot/apphot/aputil/apmapr.x
new file mode 100644
index 00000000..0d4159ab
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmapr.x
@@ -0,0 +1,21 @@
+# APMAPR -- Vector linear transformation. Map the range of pixel values
+# a1, a2 from a into the range b1, b2 into b. It is assumed that a1 < a2
+# and b1 < b2.
+
+real procedure apmapr (a, a1, a2, b1, b2)
+
+real a # the value to be mapped
+real a1, a2 # the numbers specifying the input data range
+real b1, b2 # the numbers specifying the output data range
+
+real minout, maxout, aoff, boff
+real scalar
+
+begin
+ scalar = (real (b2) - real (b1)) / (real (a2) - real (a1))
+ minout = min (b1, b2)
+ maxout = max (b1, b2)
+ aoff = a1
+ boff = b1
+ return (max(minout, min(maxout, real((a - aoff) * scalar) + boff)))
+end
diff --git a/noao/digiphot/apphot/aputil/apmeds.x b/noao/digiphot/apphot/aputil/apmeds.x
new file mode 100644
index 00000000..0a206c8a
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmeds.x
@@ -0,0 +1,154 @@
+# APSMED -- Compute the averaged median given a sorted array and an averaging
+# half width.
+
+real procedure apsmed (pix, index, npts, medcut)
+
+real pix[ARB] # array of sky pixels
+int index[ARB] # sorted index array
+int npts # number of pixels
+int medcut # averaging half width
+
+int med, j, nmed, medlo, medhi
+real sumed
+
+begin
+ med = (npts + 1) / 2
+ if (mod (med, 2) == 1) {
+ medlo = max (1, med - medcut)
+ medhi = min (npts, med + medcut)
+ } else {
+ medlo = max (1, med - medcut)
+ medhi = min (npts, med + medcut + 1)
+ }
+
+ sumed = 0.0
+ nmed = 0
+ do j = medlo, medhi {
+ sumed = sumed + pix[index[j]]
+ nmed = nmed + 1
+ }
+
+ return (sumed / nmed)
+end
+
+
+# APIMED -- Compute the index of new median value. Weight is an arbitrary
+# weight array which is assumed to be zero if the pixels has been rejected
+# and is positive otherwise.
+
+int procedure apimed (weight, index, lo, hi, nmed)
+
+real weight[ARB] # array of weights
+int index[ARB] # array of sorted indices
+int lo, hi # ranges of weights
+int nmed # number of good sky pixels
+
+int npts, med
+
+begin
+ npts = 0
+ for (med = lo; med <= hi && npts < nmed; med = med + 1) {
+ if (weight[index[med]] > 0.0)
+ npts = npts + 1
+ }
+
+ if (npts == 0)
+ return (0)
+ else
+ return (med)
+end
+
+
+# APWSMED -- Compute the new averaged median given a sorted input array,
+# an averaging half-width, and assuming that there has been pixel rejection.
+
+real procedure apwsmed (pix, index, weight, npix, med, medcut)
+
+real pix[ARB] # pixel values array
+int index[ARB] # sorted indices array
+real weight[ARB] # the weights array
+int npix # number of pixels
+int med # index of median value
+int medcut # of median cut
+
+int j, nmed, maxmed
+real sumed
+
+begin
+ sumed = pix[index[med]]
+ if (mod (med, 2) == 1)
+ maxmed = 2 * medcut + 1
+ else
+ maxmed = 2 * medcut + 2
+
+ nmed = 1
+ for (j = med - 1; j >= 1; j = j - 1) {
+ if (nmed >= medcut + 1)
+ break
+ if (weight[index[j]] > 0.0) {
+ sumed = sumed + pix[index[j]]
+ nmed = nmed + 1
+ }
+ }
+ for (j = med + 1; j <= npix; j = j + 1) {
+ if (nmed >= maxmed)
+ break
+ if (weight[index[j]] > 0.0) {
+ sumed = sumed + pix[index[j]]
+ nmed = nmed + 1
+ }
+ }
+
+ return (sumed / nmed)
+end
+
+
+# APMEDR -- Vector median selection. The selection is carried out in a temporary
+# array, leaving the input vector unmodified. Especially demanding applications
+# may wish to call the asok routine directory to avoid the call to the memory
+# allocator.
+
+real procedure apmedr (a, index, npix)
+
+real a[ARB] # input array of values
+int index[ARB] # sorted index array
+int npix # number of pixels
+
+int i
+pointer sp, aa
+real median
+real asokr() # select the Kth smallest element from A
+
+begin
+ switch (npix) {
+ case 1, 2:
+ return (a[1])
+
+ case 3:
+ if (a[1] < a[2]) {
+ if (a[2] < a[3])
+ return (a[2])
+ else if (a[1] < a[3])
+ return (a[3])
+ else
+ return (a[1])
+ } else {
+ if (a[2] > a[3])
+ return (a[2])
+ else if (a[1] < a[3])
+ return (a[1])
+ else
+ return (a[3])
+ }
+
+ default:
+ call smark (sp)
+ call salloc (aa, npix, TY_REAL)
+ do i = 1, npix
+ Memr[aa+i-1] = a[index[i]]
+ median = asokr (Memr[aa], npix, (npix + 1) / 2)
+ call sfree (sp)
+
+ return (median)
+ }
+end
diff --git a/noao/digiphot/apphot/aputil/apmoments.x b/noao/digiphot/apphot/aputil/apmoments.x
new file mode 100644
index 00000000..5a1bfbd4
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmoments.x
@@ -0,0 +1,138 @@
+# APMOMENTS -- Procedure to compute the first three moments of a
+# distribution given the appropriate sums.
+
+procedure apmoments (sumpx, sumsqpx, sumcbpx, npix, sky_zero, mean, sigma, skew)
+
+double sumpx # sum of the pixels
+double sumsqpx # sum of pixels squared
+double sumcbpx # sum of cubes of pixels
+int npix # number of pixels
+real sky_zero # sky zero point for moment analysis
+real mean # mean of pixels
+real sigma # sigma of pixels
+real skew # skew of pixels
+
+double dmean, dsigma, dskew
+
+begin
+ # Recompute the moments.
+ dmean = sumpx / npix
+ dsigma = sumsqpx / npix - dmean ** 2
+ if (dsigma <= 0.0d0) {
+ sigma = 0.0
+ skew = 0.0
+ } else {
+ dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+ sigma = sqrt (dsigma)
+ if (dskew < 0.0d0)
+ skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+ else if (dskew > 0.0d0)
+ skew = dskew ** (1.0d0 / 3.0d0)
+ else
+ skew = 0.0
+ }
+ mean = sky_zero + dmean
+end
+
+
+# APFMOMENTS -- Procedure to compute the first three moments of a distribution
+# given the data. The sums are returned as well.
+
+procedure apfmoments (pix, npix, sky_zero, sumpx, sumsqpx, sumcbpx, mean,
+ sigma, skew)
+
+real pix[npix] # array of pixels
+int npix # number of pixels
+real sky_zero # sky zero point for moment analysis
+double sumpx # sum of the pixels
+double sumsqpx # sum of pixels squared
+double sumcbpx # sum of cubes of pixels
+real mean # mean of pixels
+real sigma # sigma of pixels
+real skew # skew of pixels
+
+double dpix, dmean, dsigma, dskew
+int i
+
+begin
+ # Zero and accumulate the sums.
+ sumpx = 0.0d0
+ sumsqpx = 0.0d0
+ sumcbpx = 0.0d0
+ do i = 1, npix {
+ dpix = pix[i] - sky_zero
+ sumpx = sumpx + dpix
+ sumsqpx = sumsqpx + dpix * dpix
+ sumcbpx = sumcbpx + dpix * dpix * dpix
+ }
+
+ # Compute the moments.
+ dmean = sumpx / npix
+ dsigma = sumsqpx / npix - dmean ** 2
+ if (dsigma <= 0.0) {
+ sigma = 0.0
+ skew = 0.0
+ } else {
+ dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+ sigma = sqrt (dsigma)
+ if (dskew < 0.0d0)
+ skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+ else if (dskew > 0.0d0)
+ skew = dskew ** (1.0d0 / 3.0d0)
+ else
+ skew = 0.0
+ }
+ mean = sky_zero + dmean
+end
+
+
+# APFIMOMENTS -- Procedure to compute the first three moments of a distribution
+# given the data. The sums are returned as well.
+
+procedure apfimoments (pix, index, npix, sky_zero, sumpx, sumsqpx, sumcbpx,
+ mean, sigma, skew)
+
+real pix[ARB] # array of pixels
+int index[ARB] # the index array
+int npix # number of pixels
+real sky_zero # sky zero point for moment analysis
+double sumpx # sum of the pixels
+double sumsqpx # sum of pixels squared
+double sumcbpx # sum of cubes of pixels
+real mean # mean of pixels
+real sigma # sigma of pixels
+real skew # skew of pixels
+
+double dpix, dmean, dsigma, dskew
+int i
+
+begin
+ # Zero and accumulate the sums.
+ sumpx = 0.0d0
+ sumsqpx = 0.0d0
+ sumcbpx = 0.0d0
+ do i = 1, npix {
+ dpix = pix[index[i]] - sky_zero
+ sumpx = sumpx + dpix
+ sumsqpx = sumsqpx + dpix * dpix
+ sumcbpx = sumcbpx + dpix * dpix * dpix
+ }
+
+ # Compute the moments.
+ dmean = sumpx / npix
+ dsigma = sumsqpx / npix - dmean ** 2
+ if (dsigma <= 0.0d0) {
+ sigma = 0.0
+ skew = 0.0
+ } else {
+ dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+ sigma = sqrt (dsigma)
+ if (dskew < 0.0d0)
+ skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+ else if (dskew > 0.0d0)
+ skew = dskew ** (1.0d0 / 3.0d0)
+ else
+ skew = 0.0
+ }
+ mean = dmean + sky_zero
+end
diff --git a/noao/digiphot/apphot/aputil/apnlfuncs.x b/noao/digiphot/apphot/aputil/apnlfuncs.x
new file mode 100644
index 00000000..290b9e90
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apnlfuncs.x
@@ -0,0 +1,375 @@
+
+include <math.h>
+include <mach.h>
+
+# CGAUSS1D - Compute the value of a 1-D Gaussian function on a constant
+# background.
+
+procedure cgauss1d (x, nvars, p, np, z)
+
+real x[ARB] # variables, x[1] = position coordinate
+int nvars # the number of variables, not used
+real p[ARB] # p[1]=amplitude p[2]=center p[3]=variance p[4]=sky
+int np # number of parameters np = 4
+real z # function return
+
+real r2
+
+begin
+ if (p[3] == 0.)
+ r2 = 36.0
+ else
+ r2 = (x[1] - p[2]) ** 2 / (2. * p[3])
+ if (abs (r2) > 25.0)
+ z = p[4]
+ else
+ z = p[1] * exp (-r2) + p[4]
+end
+
+
+# CDGAUSS1D -- Compute the value a 1-D Gaussian profile on a constant
+# background and its derivatives.
+
+procedure cdgauss1d (x, nvars p, dp, np, z, der)
+
+real x[ARB] # variables, x[1] = position coordinate
+int nvars # the number of variables, not used
+real p[ARB] # p[1]=amplitude, p[2]=center, p[3]=sky p[4]=variance
+real dp[ARB] # parameter derivatives
+int np # number of parameters np=4
+real z # function value
+real der[ARB] # derivatives
+
+real dx, r2
+
+begin
+ dx = x[1] - p[2]
+ if (p[3] == 0.)
+ r2 = 36.0
+ else
+ r2 = dx * dx / (2.0 * p[3])
+ if (abs (r2) > 25.0) {
+ z = p[4]
+ der[1] = 0.0
+ der[2] = 0.0
+ der[3] = 0.0
+ der[4] = 1.0
+ } else {
+ der[1] = exp (-r2)
+ z = p[1] * der[1]
+ der[2] = z * dx / p[3]
+ der[3] = z * r2 / p[3]
+ der[4] = 1.0
+ z = z + p[4]
+ }
+end
+
+
+# GAUSSR -- Compute the value of a 2-D radially symmetric Gaussian profile
+# which is assumed to be sitting on a constant background.
+
+# Parameter Allocation:
+# 1 Amplitude
+# 2 X-center
+# 3 Y-center
+# 4 Variance
+# 5 Sky
+
+procedure gaussr (x, nvars, p, np, z)
+
+real x[ARB] # the input variables
+int nvars # the number of variables
+real p[np] # parameter vector
+int np # number of parameters
+real z # function return
+
+real dx, dy, r2
+
+begin
+ dx = x[1] - p[2]
+ dy = x[2] - p[3]
+ if (p[4] == 0.)
+ r2 = 36.0
+ else
+ r2 = (dx * dx + dy * dy) / (2.0 * p[4])
+ if (abs (r2) > 25.0)
+ z = p[5]
+ else
+ z = p[1] * exp (- r2) + p[5]
+end
+
+
+# DGAUSSR -- Compute the value of a 2-D Gaussian profile and its derivatives
+# which assumed to be sitting on top of a constant background.
+
+procedure dgaussr (x, nvars, p, dp, np, z, der)
+
+real x[ARB] # the input variables
+int nvars # the number of variables
+real p[np] # parameter vector
+real dp[np] # dummy array of parameter increments
+int np # number of parameters
+real z # function return
+real der[np] # derivatives
+
+real dx, dy, r2
+
+begin
+ dx = x[1] - p[2]
+ dy = x[2] - p[3]
+ if (p[4] == 0.)
+ r2 = 36.0
+ else
+ r2 = (dx * dx + dy * dy) / (2.0 * p[4])
+
+ if (abs (r2) > 25.0) {
+ z = p[5]
+ der[1] = 0.0
+ der[2] = 0.0
+ der[3] = 0.0
+ der[4] = 0.0
+ der[5] = 1.0
+ } else {
+ der[1] = exp (-r2)
+ z = p[1] * der[1]
+ der[2] = z * dx / p[4]
+ der[3] = z * dy / p[4]
+ der[4] = z * r2 / p[4]
+ z = z + p[5]
+ der[5] = 1.0
+ }
+end
+
+
+# ELGAUSS -- Compute the value of a 2-D elliptical Gaussian function which
+# is assumed to be sitting on top of a constant background.
+
+# Parameter Allocation:
+# 1 Amplitude
+# 2 X-center
+# 3 Y-center
+# 4 Variance-x
+# 5 Variance-y
+# 6 Theta-rotation
+# 7 Sky
+
+procedure elgauss (x, nvars, p, np, z)
+
+real x[ARB] # input variables, x[1] = x, x[2] = y
+int nvars # number of variables, not used
+real p[np] # parameter vector
+int np # number of parameters
+real z # function return
+
+real dx, dy, crot, srot, xt, yt, r2
+
+begin
+ dx = x[1] - p[2]
+ dy = x[2] - p[3]
+ crot = cos (p[6])
+ srot = sin (p[6])
+ xt = (dx * crot + dy * srot)
+ yt = (-dx * srot + dy * crot)
+ if (p[4] == 0. || p[5] == 0.)
+ r2 = 36.0
+ else
+ r2 = (xt ** 2 / p[4] + yt ** 2 / p[5]) / 2.0
+ if (abs (r2) > 25.0)
+ z = p[7]
+ else
+ z = p[1] * exp (-r2) + p[7]
+end
+
+
+# DELGAUSS -- Compute the value of a 2-D elliptical Gaussian assumed to
+# sitting on top of a constant background and its derivatives.
+
+procedure delgauss (x, nvars, p, dp, np, z, der)
+
+real x[ARB] # input variables, x[1] = x, x[2] = y
+int nvars # number of variables, not used
+real p[np] # parameter vector
+real dp[np] # delta of parameters
+int np # number of parameters
+real z # function value
+real der[np] # function return
+
+real crot, srot, crot2, srot2, sigx2, sigy2, a, b, c
+real dx, dy, dx2, dy2, r2
+
+begin
+ crot = cos (p[6])
+ srot = sin (p[6])
+ crot2 = crot ** 2
+ srot2 = srot ** 2
+ sigx2 = p[4]
+ sigy2 = p[5]
+ if (sigx2 == 0. || sigy2 == 0.)
+ r2 = 36.0
+ else {
+ a = (crot2 / sigx2 + srot2 / sigy2)
+ b = 2.0 * crot * srot * (1.0 / sigx2 - 1.0 /sigy2)
+ c = (srot2 / sigx2 + crot2 / sigy2)
+
+ dx = x[1] - p[2]
+ dy = x[2] - p[3]
+ dx2 = dx ** 2
+ dy2 = dy ** 2
+ r2 = 0.5 * (a * dx2 + b * dx * dy + c * dy2)
+ }
+
+ if (abs (r2) > 25.0) {
+ z = p[7]
+ der[1] = 0.0
+ der[2] = 0.0
+ der[3] = 0.0
+ der[4] = 0.0
+ der[5] = 0.0
+ der[6] = 0.0
+ der[7] = 1.0
+ } else {
+ der[1] = exp (-r2)
+ z = p[1] * der[1]
+ der[2] = z * (2.0 * a * dx + b * dy)
+ der[3] = z * (b * dx + 2.0 * c * dy)
+ der[4] = z * (crot2 * dx2 + 2.0 * crot * srot * dx * dy +
+ srot2 * dy2) / (2.0 * sigx2 * sigx2)
+ der[5] = z * (srot2 * dx2 - 2.0 * crot * srot * dx * dy +
+ crot2 * dy2) / (2.0 * sigy2 * sigy2)
+ der[6] = z * (b * dx2 + 2.0 * (c - a) * dx * dy - b * dy2)
+ z = z + p[7]
+ der[7] = 1.0
+ }
+end
+
+
+
+# GAUSS1D -- Compute the profile of a 1d Gaussian with a background value
+# of zero.
+
+procedure gauss1d (x, nvars, p, np, z)
+
+real x[ARB] # list of variables, x[1] = position coordinate
+int nvars # number of variables
+real p[ARB] # p[1]=amplitude p[2]=center p[3]=sigma
+int np # number of parameters == 3
+real z # function return
+
+real r
+
+begin
+ if (p[3] == 0.)
+ r = 6.0
+ else
+ r = (x[1] - p[2]) / (p[3] * SQRTOF2)
+ if (abs (r) > 5.0)
+ z = 0.0
+ else
+ z = p[1] * exp (- r ** 2)
+end
+
+
+# DGAUSS1D -- Compute the function value and derivatives of a 1-D Gaussian
+# function with a background value of zero.
+
+procedure dgauss1d (x, nvars, p, dp, np, z, der)
+
+real x[ARB] # list of variables, x[1] = position coordinate
+int nvars # number of variables
+real p[ARB] # p[1]=amplitude, p[2]=center, p[3]=sigma
+real dp[ARB] # parameter derivatives
+int np # number of parameters
+real z # function value
+real der[ARB] # derivatives
+
+real r
+
+begin
+ if (p[3] == 0.)
+ r = 6.0
+ else
+ r = (x[1] - p[2]) / (SQRTOF2 * p[3])
+ if (abs (r) > 5.0) {
+ z = 0.0
+ der[1] = 0.0
+ der[2] = 0.0
+ der[3] = 0.0
+ } else {
+ der[1] = exp (- r ** 2)
+ z = der[1] * p[1]
+ der[2] = z * r * SQRTOF2 / p[3]
+ der[3] = der[2] * SQRTOF2 * r
+ }
+end
+
+
+# GAUSSKEW - Compute the value of a 1-D skewed Gaussian profile.
+# The background value is assumed to be zero.
+
+procedure gausskew (x, nvars, p, np, z)
+
+real x[ARB] # list of variables, x[1] = position coordinate
+int nvars # number of variables, not used
+real p[ARB] # p[1]=amplitude p[2]=center p[3]=variance p[4]=skew
+int np # number of parameters == 3
+real z # function return
+
+real dx, r2, r3
+
+begin
+ dx = (x[1] - p[2])
+ if (p[3] == 0.)
+ r2 = 36.0
+ else {
+ r2 = dx ** 2 / (2.0 * p[3])
+ r3 = r2 * dx / sqrt (2.0 * abs (p[3]))
+ }
+ if (abs (r2) > 25.0)
+ z = 0.0
+ else
+ z = (1.0 + p[4] * r3) * p[1] * exp (-r2)
+end
+
+
+# DGAUSSKEW -- Compute the value of a 1-D skewed Gaussian and its derivatives.
+# The background value is assumed to be zero.
+
+procedure dgausskew (x, nvars, p, dp, np, z, der)
+
+real x[ARB] # list of variables, x[1] = position coordinate
+int nvars # number of variables, not used
+real p[ARB] # p[1]=amplitude, p[2]=center, p[3]=variance, p[4]=skew
+real dp[ARB] # parameter derivatives
+int np # number of parameters
+real z # function value
+real der[ARB] # derivatives
+
+real dx, d1, d2, d3, r, r2, r3, rint
+
+begin
+ dx = x[1] - p[2]
+ if (p[3] == 0.)
+ r2 = 36.0
+ else
+ r2 = dx ** 2 / (2.0 * p[3])
+ if (abs (r2) > 25.0) {
+ z = 0.0
+ der[1] = 0.0
+ der[2] = 0.0
+ der[3] = 0.0
+ der[4] = 0.0
+ } else {
+ r = dx / sqrt (2.0 * abs (p[3]))
+ r3 = r2 * r
+ d1 = exp (-r2)
+ z = d1 * p[1]
+ d2 = z * dx / p[3]
+ d3 = z * r2 / p[3]
+ rint = 1.0 + p[4] * r3
+ der[1] = d1 * rint
+ der[2] = d2 * (rint - 1.5 * p[4] * r)
+ der[3] = d3 * (rint - 1.5 * p[4] * r)
+ der[4] = z * r3
+ z = z * rint
+ }
+end
diff --git a/noao/digiphot/apphot/aputil/appcache.x b/noao/digiphot/apphot/aputil/appcache.x
new file mode 100644
index 00000000..d1059583
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/appcache.x
@@ -0,0 +1,83 @@
+include <imhdr.h>
+include <imset.h>
+
+# AP_MEMSTAT -- Figure out if there is enough memory to cache the image
+# pixels. If it is necessary to request more memory and the memory is
+# avalilable return YES otherwise return NO.
+
+int procedure ap_memstat (cache, req_size, old_size)
+
+int cache #I cache memory ?
+int req_size #I the requested working set size in chars
+int old_size #O the original working set size in chars
+
+int cur_size, max_size
+int begmem()
+
+begin
+ # Find the default working set size.
+ cur_size = begmem (0, old_size, max_size)
+
+ # If cacheing is disabled return NO regardless of the working set size.
+ if (cache == NO)
+ return (NO)
+
+ # If the requested working set size is less than the current working
+ # set size return YES.
+ if (req_size <= cur_size)
+ return (YES)
+
+ # Reset the current working set size.
+ cur_size = begmem (req_size, old_size, max_size)
+ if (req_size <= cur_size) {
+ return (YES)
+ } else {
+ return (NO)
+ }
+end
+
+
+# AP_PCACHE -- Cache the image pixels im memory by resetting the default image
+# buffer size. If req_size is INDEF the size of the image is used to determine
+# the size of the image i/o buffers.
+
+procedure ap_pcache (im, req_size, buf_size)
+
+pointer im #I the input image point
+int req_size #I the requested working set size in chars
+int buf_size #O the new image buffer size
+
+int def_size, new_imbufsize
+int sizeof(), imstati()
+
+begin
+ # Find the default buffer size.
+ def_size = imstati (im, IM_BUFSIZE)
+
+ # Return if the image is not 2-dimensional.
+ if (IM_NDIM(im) != 2) {
+ buf_size = def_size
+ return
+ }
+
+ # Compute the new required image i/o buffer size in chars.
+ if (IS_INDEFI(req_size)) {
+ new_imbufsize = IM_LEN(im,1) * IM_LEN(im,2) *
+ sizeof (IM_PIXTYPE(im))
+ } else {
+ new_imbufsize = req_size
+ }
+
+ # If the default image i/o buffer size is already bigger than
+ # the requested size do nothing.
+ if (def_size >= new_imbufsize) {
+ buf_size = def_size
+ return
+ }
+
+ # Reset the image i/o buffer.
+ call imseti (im, IM_BUFSIZE, new_imbufsize)
+ call imseti (im, IM_BUFFRAC, 0)
+ buf_size = new_imbufsize
+ return
+end
diff --git a/noao/digiphot/apphot/aputil/aprmwhite.x b/noao/digiphot/apphot/aputil/aprmwhite.x
new file mode 100644
index 00000000..071ea7ed
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aprmwhite.x
@@ -0,0 +1,22 @@
+include <ctype.h>
+
+# AP_RMWHITE -- Remove whitespace from a string.
+
+procedure ap_rmwhite (instr, outstr, maxch)
+
+char instr[ARB] # the input string
+char outstr[ARB] # the output string, may be the same as instr
+int maxch # maximum number of characters in outstr
+
+int ip, op
+
+begin
+ op = 1
+ for (ip = 1; (instr[ip] != EOS) && (op <= maxch); ip = ip + 1) {
+ if (IS_WHITE(instr[ip]))
+ next
+ outstr[op] = instr[ip]
+ op = op + 1
+ }
+ outstr[op] = EOS
+end
diff --git a/noao/digiphot/apphot/aputil/aptopt.x b/noao/digiphot/apphot/aputil/aptopt.x
new file mode 100644
index 00000000..0e5bc212
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aptopt.x
@@ -0,0 +1,232 @@
+include <mach.h>
+include <math.h>
+
+# APTOPT - One-dimensional centering routine using repeated convolutions to
+# locate image center.
+
+define MAX_SEARCH 3 # Max initial search steps
+
+int procedure aptopt (data, npix, center, sigma, tol, maxiter, ortho)
+
+real data[ARB] # initial data
+int npix # number of pixels
+real center # initial guess at center
+real sigma # sigma of Gaussian
+real tol # gap tolerance for sigma
+int maxiter # maximum number of iterations
+int ortho # orthogonalize weighting vector
+
+int i, iter
+pointer sp, wgt
+real newx, x[3], news, s[3], delx
+real adotr(), apqzero()
+
+begin
+ if (sigma <= 0.0)
+ return (-1)
+
+ # Allocate working space.
+ call smark (sp)
+ call salloc (wgt, npix, TY_REAL)
+
+ # Initialize.
+ x[1] = center
+ call mkt_prof_derv (Memr[wgt], npix, x[1], sigma, ortho)
+ s[1] = adotr (Memr[wgt], data, npix)
+ #if (abs (s[1]) <= EPSILONR) {
+ if (s[1] == 0.0) {
+ center = x[1]
+ call sfree (sp)
+ return (0)
+ } else
+ s[3] = s[1]
+
+ # Search for the correct interval.
+ for (i = 1; (s[3] * s[1] >= 0.0) && (i <= MAX_SEARCH); i = i + 1) {
+ s[3] = s[1]
+ x[3] = x[1]
+ x[1] = x[3] + sign (sigma, s[3])
+ call mkt_prof_derv (Memr[wgt], npix, x[1], sigma, ortho)
+ s[1] = adotr (Memr[wgt], data, npix)
+ #if (abs (s[1]) <= EPSILONR) {
+ if (s[1] == 0.0) {
+ center = x[1]
+ call sfree (sp)
+ return (0)
+ }
+ }
+
+ # Location not bracketed.
+ if (s[3] * s[1] > 0.0) {
+ call sfree (sp)
+ return (-1)
+ }
+
+ # Intialize the quadratic search.
+ delx = x[1] - x[3]
+ x[2] = x[3] - s[3] * delx / (s[1] - s[3])
+ call mkt_prof_derv (Memr[wgt], npix, x[2], sigma, ortho)
+ s[2] = adotr (Memr[wgt], data, npix)
+ #if (abs (s[2]) <= EPSILONR) {
+ if (s[2] == 0.0) {
+ center = x[2]
+ call sfree (sp)
+ return (1)
+ }
+
+ # Search quadratically.
+ for (iter = 2; iter <= maxiter; iter = iter + 1) {
+
+ # Check for completion.
+ #if (abs (s[2]) <= EPSILONR)
+ if (s[2] == 0.0)
+ break
+ if (abs (x[2] - x[1]) <= tol)
+ break
+ if (abs (x[3] - x[2]) <= tol)
+ break
+
+ # Compute new intermediate value.
+ newx = x[1] + apqzero (x, s)
+ call mkt_prof_derv (Memr[wgt], npix, newx, sigma, ortho)
+ news = adotr (Memr[wgt], data, npix)
+
+ if (s[1] * s[2] > 0.0) {
+ s[1] = s[2]
+ x[1] = x[2]
+ s[2] = news
+ x[2] = newx
+ } else {
+ s[3] = s[2]
+ x[3] = x[2]
+ s[2] = news
+ x[2] = newx
+ }
+ }
+
+ # Evaluate the center.
+ center = x[2]
+ call sfree (sp)
+ return (iter)
+end
+
+
+# AP_TPROFDER -- Procedure to estimate the approximating triangle function
+# and its derivatives.
+
+procedure ap_tprofder (data, der, npix, center, sigma, ampl)
+
+real data[ARB] # input data
+real der[ARB] # derivatives
+int npix # number of pixels
+real center # center of input Gaussian function
+real sigma # sigma of input Gaussian function
+real ampl # amplitude
+
+int i
+real x, xabs, width
+
+begin
+ width = sigma * 2.35
+ do i = 1, npix {
+ x = (i - center) / width
+ xabs = abs (x)
+ if (xabs <= 1.0) {
+ data[i] = ampl * (1.0 - xabs)
+ der[i] = x * data[i]
+ } else {
+ data[i] = 0.0
+ der[i] = 0.0
+ }
+ }
+end
+
+
+# MKT_PROF_DERV - Make orthogonal profile derivative vector.
+
+procedure mkt_prof_derv (weight, npix, center, sigma, norm)
+
+real weight[ARB] # input weight
+int npix # number of pixels
+real center # center
+real sigma # center
+int norm # orthogonalise weight
+
+pointer sp, der
+real coef
+real asumr(), adotr()
+
+begin
+ call smark (sp)
+ call salloc (der, npix, TY_REAL)
+
+ # Fetch the weighting function and derivatives.
+ call ap_tprofder (Memr[der], weight, npix, center, sigma, 1.0)
+
+ if (norm == YES) {
+
+ # Make orthogonal to level background.
+ coef = -asumr (weight, npix) / npix
+ call aaddkr (weight, coef, weight, npix)
+ coef = -asumr (Memr[der], npix) / npix
+ call aaddkr (Memr[der], coef, Memr[der], npix)
+
+ # Make orthogonal to profile vector.
+ coef = adotr (Memr[der], Memr[der], npix)
+ if (coef <= 0.0)
+ coef = 1.0
+ else
+ coef = adotr (weight, Memr[der], npix) / coef
+ call amulkr (Memr[der], coef, Memr[der], npix)
+ call asubr (weight, Memr[der], weight, npix)
+
+ # Normalize the final vector.
+ coef = adotr (weight, weight, npix)
+ if (coef <= 0.0)
+ coef = 1.0
+ else
+ coef = sqrt (1.0 / coef)
+ call amulkr (weight, coef, weight, npix)
+ }
+
+ call sfree (sp)
+end
+
+define QTOL .125
+
+# APQZERO - Solve for the root of a quadratic function defined by three
+# points.
+
+real procedure apqzero (x, y)
+
+real x[3]
+real y[3]
+
+real a, b, c, det, dx
+real x2, x3, y2, y3
+
+begin
+ # Compute the determinant.
+ x2 = x[2] - x[1]
+ x3 = x[3] - x[1]
+ y2 = y[2] - y[1]
+ y3 = y[3] - y[1]
+ det = x2 * x3 * (x2 - x3)
+
+ # Compute the shift in x.
+ #if (abs (det) > 100.0 * EPSILONR) {
+ if (abs (det) > 0.0) {
+ a = (x3 * y2 - x2 * y3) / det
+ b = - (x3 * x3 * y2 - x2 * x2 * y3) / det
+ c = a * y[1] / (b * b)
+ if (abs (c) > QTOL)
+ dx = (-b / (2.0 * a)) * (1.0 - sqrt (1.0 - 4.0 * c))
+ else
+ dx = - (y[1] / b) * (1.0 + c)
+ return (dx)
+ #} else if (abs (y3) > EPSILONR)
+ } else if (abs (y3) > 0.0)
+ return (-y[1] * x3 / y3)
+ else
+ return (0.0)
+end
diff --git a/noao/digiphot/apphot/aputil/apvectors.x b/noao/digiphot/apphot/aputil/apvectors.x
new file mode 100644
index 00000000..7c9dbcfd
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apvectors.x
@@ -0,0 +1,515 @@
+
+# AP_ABOXR -- Vector boxcar smooth. The input vector is convolved with a
+# uniform kernel of length knpix. Boundary extension is assumed to have
+# been performed on the input array before entering the routine.
+
+procedure ap_aboxr (in, out, npix, knpix)
+
+real in[npix+knpix-1] # input array with boundary exntension
+real out[npix] # output array
+int npix # number of pixels
+int knpix # length of the kernel
+
+int i
+real sum
+
+begin
+ sum = 0.0
+ do i = 1, knpix - 1
+ sum = sum + in[i]
+
+ do i = 1, npix {
+ sum = sum + in[i+knpix-1]
+ out[i] = sum / knpix
+ sum = sum - in[i]
+ }
+
+end
+
+
+# AP_SBOXR -- Boxcar smooth for a vector that has not been boundary
+# extended.
+
+procedure ap_sboxr (in, out, npix, nsmooth)
+
+real in[npix] # input array
+real out[npix] # output array
+int npix # number of pixels
+int nsmooth # half width of smoothing box
+
+int i, j, ib, ie, ns
+real sum
+
+begin
+ ns = 2 * nsmooth + 1
+ do i = 1, npix {
+ ib = max (i - nsmooth, 1)
+ ie = min (i + nsmooth, npix)
+ sum = 0.0
+ do j = ib, ie
+ sum = sum + in[j]
+ out[i] = sum / ns
+ }
+end
+
+
+# AP_ALIMR -- Compute the maximum and minimum data values and indices of a
+# 1D array.
+
+procedure ap_alimr (data, npts, mindat, maxdat, imin, imax)
+
+real data[npts] # data array
+int npts # number of points
+real mindat, maxdat # min and max data value
+int imin, imax # indices of min and max data values
+
+int i
+
+begin
+ imin = 1
+ imax = 1
+ mindat = data[1]
+ maxdat = data[1]
+
+ do i = 2, npts {
+ if (data[i] > maxdat) {
+ imax = i
+ maxdat = data[i]
+ }
+ if (data[i] < mindat) {
+ imin = i
+ mindat = data[i]
+ }
+ }
+end
+
+
+# AP_IALIMR -- Compute the maximum and minimum data values of a 1D indexed
+# array
+
+procedure ap_ialimr (data, index, npts, mindat, maxdat)
+
+real data[npts] # data array
+int index[npts] # index array
+int npts # number of points
+real mindat, maxdat # min and max data value
+
+int i
+
+begin
+ mindat = data[index[1]]
+ maxdat = data[index[1]]
+ do i = 2, npts {
+ if (data[index[i]] > maxdat)
+ maxdat = data[index[i]]
+ if (data[index[i]] < mindat)
+ mindat = data[index[i]]
+ }
+end
+
+
+# AP_ASUMR - Compute the sum of an index sorted array
+
+real procedure ap_asumr (data, index, npts)
+
+real data[npts] # data array
+int index[npts] # index array
+int npts # number of points
+
+double sum
+int i
+
+begin
+ sum = 0.0d0
+ do i = 1, npts
+ sum = sum + data[index[i]]
+ return (real (sum))
+end
+
+
+# AP_AMAXEL -- Find the maximum value and its index of a 1D array.
+
+procedure ap_amaxel (a, npts, maxdat, imax)
+
+real a[ARB] # the data array
+int npts # number of points
+real maxdat # maximum value
+int imax # imdex of max value
+
+int i
+
+begin
+ maxdat = a[1]
+ imax = 1
+ do i = 2, npts {
+ if (a[i] > maxdat) {
+ maxdat = a[i]
+ imax = i
+ }
+ }
+end
+
+
+# AHGMR -- Accumulate the histogram of the input vector. The output vector
+# hgm (the histogram) should be cleared prior to the first call. The procedure
+# returns the number of data values it could not include in the histogram.
+
+int procedure aphgmr (data, wgt, npix, hgm, nbins, z1, z2)
+
+real data[ARB] # data vector
+real wgt[ARB] # weights vector
+int npix # number of pixels
+real hgm[ARB] # output histogram
+int nbins # number of bins in histogram
+real z1, z2 # greyscale values of first and last bins
+
+real dz
+int bin, i, nreject
+
+begin
+ if (nbins < 2)
+ return (0)
+
+ nreject = 0
+ dz = real (nbins - 1) / real (z2 - z1)
+ do i = 1, npix {
+ if (data[i] < z1 || data[i] > z2) {
+ nreject = nreject + 1
+ wgt[i] = 0.0
+ next
+ }
+ bin = int ((data[i] - z1) * dz) + 1
+ hgm[bin] = hgm[bin] + 1.0
+ }
+
+ return (nreject)
+end
+
+
+# APHIGMR -- Accumulate the histogram of the input vector. The output vector
+# hgm (the histogram) should be cleared prior to the first call. The procedure
+# returns the number of data values it could not include in the histogram.
+
+int procedure aphigmr (data, wgt, index, npix, hgm, nbins, z1, z2)
+
+real data[ARB] # data vector
+real wgt[ARB] # weights vector
+int index[ARB] # index vector
+int npix # number of pixels
+real hgm[ARB] # output histogram
+int nbins # number of bins in histogram
+real z1, z2 # greyscale values of first and last bins
+
+real dz
+int bin, i, nreject
+
+begin
+ if (nbins < 2)
+ return (0)
+
+ nreject = 0
+ dz = real (nbins - 1) / real (z2 - z1)
+ do i = 1, npix {
+ if (data[index[i]] < z1 || data[index[i]] > z2) {
+ nreject = nreject + 1
+ wgt[index[i]] = 0.0
+ next
+ }
+ bin = int ((data[index[i]] - z1) * dz) + 1
+ hgm[bin] = hgm[bin] + 1.0
+ }
+
+ return (nreject)
+end
+
+
+# AP_IJTOR -- Compute radius values given the center coordinates, the line
+# number in a the subraster and the length of the subraster x axis.
+
+procedure ap_ijtor (r, nr, line, xc, yc)
+
+real r[nr] # array of output r values
+int nr # length of r array
+int line # line number
+real xc, yc # subraster center
+
+int i
+real temp
+
+begin
+ temp = (line - yc ) ** 2
+ do i = 1, nr
+ r[i] = sqrt ((i - xc) ** 2 + temp)
+end
+
+
+# AP_IJTOR2 -- Compute radius values given the center coordinates and the size
+# of the subraster.
+
+procedure ap_ijtor2 (r, nx, ny, xc, yc)
+
+real r[nx,ARB] # array of output r values
+int nx # x dimension of output array
+int ny # Y dimension of output array
+real xc, yc # subraster center
+
+int i, j
+real temp
+
+begin
+ do j = 1, ny {
+ temp = (j - yc) ** 2
+ do i = 1, nx
+ r[i,j] = sqrt ((i - xc) ** 2 + temp)
+ }
+end
+
+
+# AP_2DALIMR -- Compute min and max values of a 2D array along with
+# the indices of the minimum and maximum pixel.
+
+procedure ap_2dalimr (data, nx, ny, mindat, maxdat, imin, jmin, imax, jmax)
+
+real data[nx, ny] # data
+int nx, ny # array dimensions
+real mindat, maxdat # min, max data values
+int imin, jmin # indices of min element
+int imax, jmax # indices of max element
+
+int i, j
+
+begin
+ imin = 1
+ jmin = 1
+ imax = 1
+ jmax = 1
+ mindat = data[1,1]
+ maxdat = data[1,1]
+
+ do j = 2, ny {
+ do i = 1, nx {
+ if (data[i,j] < mindat) {
+ imin = i
+ jmin = j
+ mindat = data[i,j]
+ } else if (data[i,j] > maxdat) {
+ imax = i
+ jmax = j
+ maxdat = data[i,j]
+ }
+ }
+ }
+end
+
+
+define LOGPTR 20 # log2(maxpts) (1e6)
+
+# APQSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+procedure apqsort (data, a, b, npix)
+
+real data[ARB] # data array
+int a[ARB], b[ARB] # index array
+int npix # number of pixels
+
+int i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+real pivot
+
+begin
+ # Initialize the indices for an inplace sort.
+ do i = 1, npix
+ a[i] = i
+ call amovi (a, b, npix)
+
+ p = 1
+ lv[1] = 1
+ uv[1] = npix
+ while (p > 0) {
+
+ # If only one elem in subset pop stack otherwise pivot line.
+ if (lv[p] >= uv[p])
+ p = p - 1
+ else {
+ i = lv[p] - 1
+ j = uv[p]
+ pivot = data[b[j]]
+
+ while (i < j) {
+ for (i=i+1; data[b[i]] < pivot; i=i+1)
+ ;
+ for (j=j-1; j > i; j=j-1)
+ if (data[b[j]] <= pivot)
+ break
+ if (i < j) { # out of order pair
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+ }
+ }
+
+ j = uv[p] # move pivot to position i
+ temp = b[j] # interchange elements
+ b[j] = b[i]
+ b[i] = temp
+
+ if (i-lv[p] < uv[p] - i) { # stack so shorter done first
+ lv[p+1] = lv[p]
+ uv[p+1] = i - 1
+ lv[p] = i + 1
+ } else {
+ lv[p+1] = i + 1
+ uv[p+1] = uv[p]
+ uv[p] = i - 1
+ }
+
+ p = p + 1 # push onto stack
+ }
+ }
+end
+
+
+
+# APRECIPROCAL -- Compute the reciprocal of the absolute value of a vector.
+
+procedure apreciprocal (a, b, npts, value)
+
+real a[ARB] # the input vector
+real b[ARB] # the output vector
+int npts # the number of data points
+real value # the value to be assigned to b[i] if a[i] = 0
+
+int i
+
+begin
+ do i = 1, npts {
+ if (a[i] > 0.0)
+ b[i] = 1.0 / a[i]
+ else if (a[i] < 0.0)
+ b[i] = - 1.0 / a[i]
+ else
+ b[i] = value
+ }
+end
+
+
+# AP_WLIMR -- Set the weights of all data points outside a given minimum and
+# maximum values to zero. Compute the minimum and maximum values and their
+# indices of the remaining data.
+
+procedure ap_wlimr (pix, w, npts, datamin, datamax, dmin, dmax, imin, imax)
+
+real pix[ARB] # input pixel array
+real w[ARB] # weight array
+int npts # number of points
+real datamin # minimum good data point
+real datamax # maximum good data point
+real dmin # output data minimum
+real dmax # output data maximum
+int imin # index of the data minimum
+int imax # index of the data maximum
+
+int i
+real value
+
+begin
+ dmin = datamax
+ dmax = datamin
+ imin = 1
+ imax = 1
+
+ do i = 1, npts {
+ value = pix[i]
+ if ((value < datamin) || (value > datamax)) {
+ w[i] = 0.0
+ next
+ }
+ if (value < dmin) {
+ dmin = value
+ imin = i
+ } else if (value > dmax) {
+ dmax = value
+ imax = i
+ }
+ }
+end
+
+
+# APWSSQR -- Compute the weighted sum of the squares of a vector.
+
+real procedure apwssqr (a, wgt, npts)
+
+real a[ARB] # data array
+real wgt[ARB] # array of weights
+int npts # number of points
+
+int i
+real sum
+
+begin
+ sum = 0.0
+ do i = 1, npts
+ sum = sum + wgt[i] * a[i] * a[i]
+ return (sum)
+end
+
+
+# AP_XYTOR -- Change the single integer coord of the sky pixels to a radial
+# distance value. The integer coordinate is equal to coord = (i - xc + 1) +
+# blklen * (j - yc).
+
+procedure ap_xytor (coords, index, r, nskypix, xc, yc, blklen)
+
+int coords[ARB] # coordinate array
+int index[ARB] # the index array
+real r[ARB] # radial coordinates
+int nskypix # number of sky pixels
+real xc, yc # center of sky subraster
+int blklen # x dimension of sky subraster
+
+int i
+real x, y
+
+begin
+ do i = 1, nskypix {
+ x = real (mod (coords[index[i]], blklen))
+ if (x == 0)
+ x = real (blklen)
+ y = (coords[index[i]] - x) / blklen + 1
+ r[i] = sqrt ((x - xc) ** 2 + (y - yc) ** 2)
+ }
+end
+
+
+# AP_W1SUR -- Linearly combine two vectors where the weight for the first
+# vectors is 1.0 and is k2 for the second vector
+
+procedure ap_w1sur (a, b, c, npts, k2)
+
+real a[ARB] # the first input vector
+real b[ARB] # the second input vector
+real c[ARB] # the output vector
+int npts # number of points
+real k2 # the weigting factor for the second vector
+
+int i
+
+begin
+ do i = 1, npts
+ c[i] = a[i] + k2 * b[i]
+end
+
+
+# AP_INDEX -- Define an index array.
+
+procedure ap_index (index, npix)
+
+int index[ARB] # the index array
+int npix # the number of pixels
+
+int i
+
+begin
+ do i = 1, npix
+ index[i] = i
+end
diff --git a/noao/digiphot/apphot/aputil/mkpkg b/noao/digiphot/apphot/aputil/mkpkg
new file mode 100644
index 00000000..eefeed9d
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/mkpkg
@@ -0,0 +1,25 @@
+# Miscellaneous Routines used by the APPHOT Package
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+ apbsmooth.x
+ apclip.x
+ apdate.x <time.h>
+ apdiverr.x
+ apfnames.x
+ apgscur.x <gset.h> <fset.h>
+ apgtools.x <pkg/gtools.h> <gset.h>
+ aplucy.x
+ apnlfuncs.x <mach.h> <math.h>
+ apmapr.x
+ apmeds.x
+ apmoments.x
+ appcache.x <imhdr.h> <imset.h>
+ aprmwhite.x <ctype.h>
+ aptopt.x <math.h> <mach.h>
+ apvectors.x
+ ;
generated by cgit (git 2.34.1) at 2025-09-20 09:01:04 -0400