aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/daophot/allstar/dpalinit.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/daophot/allstar/dpalinit.x')
-rw-r--r--noao/digiphot/daophot/allstar/dpalinit.x665
1 files changed, 665 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/allstar/dpalinit.x b/noao/digiphot/daophot/allstar/dpalinit.x
new file mode 100644
index 00000000..4f1afb89
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalinit.x
@@ -0,0 +1,665 @@
+include <mach.h>
+include <imhdr.h>
+include <math.h>
+include "../lib/daophotdef.h"
+include "../lib/allstardef.h"
+
+# DP_SETWT -- Initialize the subtracted image and compute the initial weights
+# image from the input image.
+
+procedure dp_setwt (dao, im)
+
+pointer dao # pointer to daophot structure
+pointer im # input image decriptor
+
+int i, ncol, nline
+pointer sp, v1, v2, v3, v4, line1, line2, line3, line4
+pointer allstar, wtim, dataim, subtim
+real rdnoise, mingdata, maxgdata
+int imgnlr(), impnlr()
+
+begin
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (v3, IM_MAXDIM, TY_LONG)
+ call salloc (v4, IM_MAXDIM, TY_LONG)
+
+ # Define the allstar pointer.
+ allstar = DP_ALLSTAR (dao)
+
+ # Define some useful constants.
+ rdnoise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+ if (IS_INDEFR(DP_MINGDATA(dao)))
+ mingdata = -MAX_REAL
+ else
+ mingdata = DP_MINGDATA(dao)
+ if (IS_INDEFR(DP_MAXGDATA(dao)))
+ maxgdata = MAX_REAL
+ else
+ maxgdata = DP_MAXGDATA(dao)
+
+ wtim = DP_WEIGHTS (allstar)
+ dataim = DP_DATA(allstar)
+ subtim = DP_SUBT(allstar)
+ ncol = IM_LEN (im,1)
+ nline = IM_LEN(im,2)
+
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovkl (long(1), Meml[v2], IM_MAXDIM)
+ call amovkl (long(1), Meml[v3], IM_MAXDIM)
+ call amovkl (long(1), Meml[v4], IM_MAXDIM)
+ do i = 1, nline {
+
+ # Get the next line
+ if (imgnlr (im, line1, Meml[v1]) == EOF)
+ break
+
+ # Initialize the subtracted image.
+ if (DP_CACHE(allstar, A_DCOPY) == YES) {
+ call amovr (Memr[line1], Memr[dataim], ncol)
+ dataim = dataim + ncol
+ } else if (impnlr (dataim, line2, Meml[v2]) != EOF) {
+ call amovr (Memr[line1], Memr[line2], ncol)
+ }
+
+ # Initialize the weights image.
+ if (DP_CACHE(allstar, A_WEIGHT) == YES) {
+ call dp_wtvector (Memr[line1], Memr[wtim], ncol, mingdata,
+ maxgdata, rdnoise)
+ wtim = wtim + ncol
+ } else if (impnlr (wtim, line3, Meml[v3]) != EOF) {
+ call dp_wtvector (Memr[line1], Memr[line3], ncol, mingdata,
+ maxgdata, rdnoise)
+ }
+
+ # Initilize the subtracted image.
+ if (DP_CACHE(allstar, A_WEIGHT) == YES) {
+ ;
+ } else if (impnlr (subtim, line4, Meml[v4]) != EOF) {
+ call amovkr (0.0, Memr[line4], ncol)
+ }
+
+ }
+
+ # Make sure all the changes are written to disk.
+ if (DP_CACHE(allstar, A_WEIGHT) == NO)
+ call imflush (DP_WEIGHTS(allstar))
+ if (DP_CACHE(allstar, A_DCOPY) == NO)
+ call imflush (DP_DATA(allstar))
+ if (DP_CACHE(allstar, A_SUBT) == NO)
+ call imflush (DP_SUBT(allstar))
+
+ call sfree (sp)
+end
+
+
+# DP_WTVECTOR -- Compute the initial weights for a vector of input data.
+
+procedure dp_wtvector (a, b, ncols, mingdata, maxgdata, rnoisesq)
+
+real a[ARB] # input array
+real b[ARB] # output array
+int ncols # number of points
+real mingdata # minimum good data value
+real maxgdata # maximum good data value
+real rnoisesq # read noise squared in adu
+
+int i
+
+begin
+ do i = 1, ncols {
+ if (a[i] < mingdata || a[i] > maxgdata)
+ b[i] = -MAX_REAL
+ else
+ b[i] = rnoisesq
+ }
+end
+
+
+define DELTA_MAG 12.5
+define INIT_REL_BRIGHT 0.003
+
+# DP_ALZERO -- Convert from magnitudes to relative brightnesses before
+# fitting the stars. If the star's magnitude is undefined or more than
+# 12.5 magnitudes fainter than the PSF magnitude, a default relative
+# brightness is defined.
+
+procedure dp_alzero (dao, mag, nstars)
+
+pointer dao # pointer to the daophot strucuture
+real mag[ARB] # the magnitude array
+int nstars # number of stars
+
+int i
+pointer psffit
+real faint
+
+begin
+ psffit = DP_PSFFIT (dao)
+ faint = DP_PSFMAG(psffit) + DELTA_MAG
+
+ do i = 1, nstars {
+ if (IS_INDEFR(mag[i])) {
+ mag[i] = INIT_REL_BRIGHT
+ } else if (mag[i] >= faint) {
+ mag[i] = INIT_REL_BRIGHT
+ } else {
+ mag[i] = DAO_RELBRIGHT (psffit, mag[i])
+ }
+ }
+end
+
+
+# DP_STRIP -- Remove stars with undefined centers, stars with centers that
+# are off the input image, and duplicate stars from the input photometry
+# list, where duplicate stars are defined as those within a distance of
+# radius pixels of each other. The duplicate stars are moved to the end of
+# the star list and the number of stars is recomputed.
+
+procedure dp_strip (id, x, y, mag, sky, skip, aier, nstar, sepradsq, ncol,
+ nline, fitradius, verbose)
+
+int id[ARB] # array of star ids
+real x[ARB] # array of star x values
+real y[ARB] # array of star y values
+real mag[ARB] # array of star magnitudes
+real sky[ARB] # array of star sky values
+int skip[ARB] # array of fit/nofit indicators
+int aier[ARB] # array of error codes
+int nstar # number of stars (may change)
+real sepradsq # separation radius squared
+int ncol # number of columns in the image
+int nline # number of columns in the image
+real fitradius # the fitting radius
+int verbose # print error messages ?
+
+int i, j, ier, ahold
+pointer sp, index, idhold
+real seprad, dx, dy, xhold, yhold, maghold, skyhold
+
+begin
+ # Duplicate stars are impossible.
+ if (nstar <= 1)
+ return
+
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (index, nstar, TY_INT)
+
+ # Sort the data in y.
+ call quick (y, nstar, Memi[index], ier)
+
+ # Rectify the remaining arrays.
+ call dp_irectify (id, Memi[index], nstar)
+ call dp_rectify (x, Memi[index], nstar)
+ call dp_rectify (mag, Memi[index], nstar)
+ call dp_rectify (sky, Memi[index], nstar)
+
+ # Determine whether any stars are close to star i.
+ seprad = sqrt (sepradsq)
+ do i = 1, nstar - 1 {
+
+ # This star was rejected on a previous loop.
+ if (skip[i] == YES)
+ next
+
+ # Reject if star has an INDEF valued position or is off image.
+ if (IS_INDEFR(x[i]) || IS_INDEFR(y[i])) {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d has an undefined x or y value\n")
+ call pargi (id[i])
+ } else if ((int (x[i] - fitradius) + 1) > ncol) {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf ("REJECTING: Star %d is outside the image\n")
+ call pargi (id[i])
+ } else if (int (x[i] + fitradius) < 1) {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf ("REJECTING: Star %d is outside the image\n")
+ call pargi (id[i])
+ } else if ((int (y[i] - fitradius) + 1) > nline) {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf ("REJECTING: Star %d is outside the image\n")
+ call pargi (id[i])
+ } else if (int (y[i] + fitradius) < 1) {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf ("REJECTING: Star %d is outside the image\n")
+ call pargi (id[i])
+ }
+
+ # This star was rejected on this loop.
+ if (skip[i] == YES)
+ next
+
+ # Loop over the remaining stars.
+ do j = i + 1, nstar {
+
+ # Star has already been rejected on previous loop.
+ if (skip[j] == YES)
+ next
+
+ # Test for INDEF.
+ if (IS_INDEFR(x[j]) || IS_INDEFR(y[j])) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d has an undefined x or y value\n")
+ call pargi (id[j])
+ } else if ((int (x[j] - fitradius) + 1) > ncol) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d is outside the image\n")
+ call pargi (id[j])
+ } else if (int (x[j] + fitradius) < 1) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d is outside the image\n")
+ call pargi (id[j])
+ } else if ((int (y[j] - fitradius) + 1) > nline) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d is outside the image\n")
+ call pargi (id[j])
+ } else if (int (y[j] + fitradius) < 1) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_OFFIMAGE
+ if (verbose == YES)
+ call printf (
+ "REJECTING: Star %d is outside the image\n")
+ call pargi (id[j])
+ }
+
+ # Star was rejected.
+ if (skip[j] == YES)
+ next
+
+ # Test for proximity.
+ dy = y[j] - y[i]
+ if (dy > seprad)
+ break
+ dx = x[j] - x[i]
+ if (abs (dx) > seprad)
+ next
+ if ((dx * dx + dy * dy) > sepradsq)
+ next
+
+ # Set the magnitude of the star to INDEF and skip.
+ if (mag[j] <= mag[i]) {
+ mag[j] = INDEFR
+ skip[j] = YES
+ aier[j] = ALLERR_MERGE
+ if (verbose == YES) {
+ call printf (
+ "REJECTING: Star %d has merged with star %d\n")
+ call pargi (id[j])
+ call pargi (id[i])
+ }
+ } else {
+ mag[i] = INDEFR
+ skip[i] = YES
+ aier[i] = ALLERR_MERGE
+ if (verbose == YES) {
+ call printf (
+ "REJECTING: Star %d has merged with star %d\n")
+ call pargi (id[i])
+ call pargi (id[j])
+ }
+ break
+ }
+ }
+ }
+
+ # Remove the duplicate stars.
+ for (i = 1; i <= nstar; i = i + 1) {
+
+ # Redefine the number of stars by removing skipped stars from
+ # the end of the star list.
+ while (nstar >= 1) {
+ if (skip[nstar] == NO)
+ break
+ nstar = nstar - 1
+ }
+ if (i > nstar)
+ break
+ if (skip[i] == NO)
+ next
+
+ # Switch the rejected star with the one at the end of the list.
+ idhold = id[i]
+ xhold = x[i]
+ yhold = y[i]
+ maghold = mag[i]
+ skyhold = sky[i]
+ ahold = aier[i]
+
+ id[i] = id[nstar]
+ x[i] = x[nstar]
+ y[i] = y[nstar]
+ mag[i] = mag[nstar]
+ sky[i] = sky[nstar]
+ skip[i] = NO
+ aier[i] = aier[nstar]
+
+ id[nstar] = idhold
+ x[nstar] = xhold
+ y[nstar] = yhold
+ mag[nstar] = maghold
+ sky[nstar] = skyhold
+ skip[nstar] = YES
+ aier[nstar] = ahold
+
+ nstar = nstar - 1
+ }
+
+ call sfree (sp)
+end
+
+
+# DP_WSTINIT -- Initialize the weight and scratch arrays / images.
+
+procedure dp_wstinit (dao, im, xcen, ycen, mag, nstar, radius, x1, x2, y1, y2)
+
+pointer dao # pointer to the daophot structure
+pointer im # pointer to the input image
+real xcen[ARB] # the x centers array
+real ycen[ARB] # the y centers array
+real mag[ARB] # the magnitude array
+int nstar # the number of stars
+real radius # radius for subraction
+int x1, x2 # column limits
+int y1, y2 # line limits
+
+int j, l, ncol, npix, nl, ns, lx, mx
+pointer psffit, allstar, databuf, wtbuf, owtbuf, subtbuf
+real rsq, psfradius, psfradsq, x, y, dy, dysq, deltax, deltay
+pointer imgs2r(), imps2r()
+
+begin
+ # Set up some pointers.
+ psffit = DP_PSFFIT(dao)
+ allstar = DP_ALLSTAR(dao)
+
+ # Set up some constants.
+ ncol = IM_LEN(im,1)
+ npix = x2 - x1 + 1
+ rsq = radius ** 2
+ if (DP_PSFSIZE(psffit) == 0)
+ psfradius = DP_PSFRAD(dao)
+ else
+ psfradius = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0
+ psfradsq = psfradius * psfradius
+
+ # Begin the search of the groups at the first group.
+ ns = 0
+ nl = 0
+
+ ns = 0
+ do j = y1, y2 {
+
+ # Get the data.
+ if (DP_CACHE(allstar,A_DCOPY) == YES)
+ databuf = DP_DATA(allstar) + (j - 1) * ncol + x1 - 1
+ else
+ databuf = imgs2r (DP_DATA(allstar), x1, x2, j, j)
+ if (DP_CACHE(allstar,A_WEIGHT) == YES) {
+ wtbuf = DP_WEIGHTS(allstar) + (j - 1) * ncol + x1 - 1
+ owtbuf = wtbuf
+ } else {
+ owtbuf = imps2r (DP_WEIGHTS(allstar), x1, x2, j, j)
+ wtbuf = imgs2r (DP_WEIGHTS(allstar), x1, x2, j, j)
+ }
+ if (DP_CACHE(allstar,A_SUBT) == YES)
+ subtbuf = DP_SUBT(allstar) + (j - 1) * ncol + x1 - 1
+ else
+ subtbuf = imps2r (DP_SUBT(allstar), x1, x2, j, j)
+
+ # Set all the weights in the working array negative.
+ call aabsr (Memr[wtbuf], Memr[owtbuf], npix)
+ call anegr (Memr[owtbuf], Memr[owtbuf], npix)
+ call amovkr (0.0, Memr[subtbuf], npix)
+
+ # Set the y coordinate.
+ y = real (j)
+
+ # Initialize the weight and scratch arrays.
+
+ # Find all the points within one working radius of each star.
+ # Set all the sigmas positive again and copy them from DATA
+ # to SUBT.
+ do l = ns + 1, nstar {
+ dy = y - ycen[l]
+ if (dy > radius) {
+ ns = l
+ next
+ } else if (dy < -radius)
+ break
+ dysq = dy ** 2
+ lx = max (x1, min (x2, int (xcen[l] - radius) + 1))
+ mx = max (x1, min (x2, int (xcen[l] + radius)))
+ x = xcen[l] - lx + 1.0
+ lx = lx - x1 + 1
+ mx = mx - x1 + 1
+ call dp_wstvector (x, Memr[databuf+lx-1], Memr[owtbuf+lx-1],
+ Memr[subtbuf+lx-1], mx - lx + 1, dysq, rsq, -MAX_REAL)
+ }
+
+ do l = nl + 1, nstar {
+ dy = y - ycen[l]
+ if (dy > psfradius) {
+ nl = l
+ next
+ } else if (dy < -psfradius)
+ break
+ dysq = dy ** 2
+ lx = max (x1, min (x2, int (xcen[l] - psfradius) + 1))
+ mx = max (x1, min (x2, int (xcen[l] + psfradius)))
+ x = xcen[l] - lx + 1.0
+ lx = lx - x1 + 1
+ mx = mx - x1 + 1
+ call dp_wpsf (dao, im, xcen[l], ycen[l], deltax, deltay, 1)
+ deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+ deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+ call dp_alsubstar (psffit, x, dy, deltax, deltay, mag[l],
+ Memr[subtbuf+lx-1], Memr[owtbuf+lx-1], mx - lx + 1, dysq,
+ psfradsq)
+ }
+ }
+
+ if (DP_CACHE(allstar,A_WEIGHT) == NO)
+ call imflush (DP_WEIGHTS(allstar))
+ if (DP_CACHE(allstar,A_SUBT) == NO)
+ call imflush (DP_SUBT(allstar))
+end
+
+
+# DP_WSTVECTOR -- Set the weight and scratch vector
+
+procedure dp_wstvector (xcen, data, weight, subt, npix, dysq, rsq, badwt)
+
+real xcen # x coordinate of center
+real data[ARB] # the data array
+real weight[ARB] # the weight array
+real subt[ARB] # the subtracted array array
+int npix # the number of pixels
+real dysq # the y distance squared
+real rsq # the inclusion radius squared
+real badwt # badwt value
+
+int i
+real dx
+
+begin
+ do i = 1, npix {
+ dx = (real (i) - xcen)
+ if ((dx ** 2 + dysq) <= rsq) {
+ if (weight[i] <= badwt)
+ next
+ weight[i] = abs (weight[i])
+ subt[i] = data[i]
+ } else if (dx >= 0.0)
+ break
+ }
+end
+
+
+# DP_ALSUBSTAR -- The subtraction vector.
+
+procedure dp_alsubstar (psffit, xcen, dy, deltax, deltay, mag, subt, weight,
+ npix, dysq, psfradsq)
+
+pointer psffit # pointer to the psf fitting structure
+real xcen # x coordinate of center
+real dy # y offset from psf center
+real deltax # x distance from psf position
+real deltay # y distance from psf position
+real mag # the magnitude of the star
+real subt[ARB] # the subtracted array array
+real weight[ARB] # the weight array
+int npix # the number of pixels
+real dysq # the y distance squared
+real psfradsq # the inclusion radius squared
+
+int i
+real dx, dvdxc, dvdyc
+real dp_usepsf()
+
+begin
+ do i = 1, npix {
+ if (weight[i] < 0.0)
+ next
+ dx = real (i) - xcen
+ if ((dx ** 2 + dysq) < psfradsq) {
+ subt[i] = subt[i] - mag * dp_usepsf (DP_PSFUNCTION(psffit),
+ dx, dy, DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+ Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit),
+ DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), deltax, deltay,
+ dvdxc, dvdyc)
+ } else if (dx > 0.0)
+ break
+ }
+end
+
+
+# DP_ALSKY -- Recompute the sky values.
+
+procedure dp_alsky (dao, im, xcen, ycen, sky, nstar, x1, x2, y1, y2, rinsky,
+ routsky, minnsky, badwt)
+
+pointer dao # pointer to the daophot structure
+int im # the input image pointer
+real xcen[ARB] # x coordinates
+real ycen[ARB] # y coordinates
+real sky[ARB] # sky values
+int nstar # number of stars
+int x1, x2 # coordinate limits
+int y1, y2 # coordinate limits
+real rinsky # inner radius of the sky annulus
+real routsky # outer radius of the sky annulus
+int minnsky # minimum number of sky pixels
+real badwt # the bad weight value
+
+int istar, i, j, ncols, lenskybuf, lx, mx, ly, my, line1, line2, nsky, ier
+pointer allstar, sub, psub, wgt, pwgt, skyvals, index
+real risq, rosq, dannulus, dysq, dx, rsq
+pointer dp_gst(), dp_gwt()
+
+begin
+ # Get the allstar pointer.
+ allstar = DP_ALLSTAR(dao)
+
+ # Define some constants
+ ncols = IM_LEN(im,1)
+ risq = rinsky ** 2
+ rosq = routsky ** 2
+ dannulus = routsky - rinsky
+
+ # Allcoate some working memory.
+ lenskybuf = PI * (2.0 * rinsky + dannulus + 1.0) * (dannulus + 0.5)
+ call malloc (skyvals, lenskybuf, TY_REAL)
+ call malloc (index, lenskybuf, TY_INT)
+
+ # Accumulate the sky buffer.
+ do istar = 1, nstar {
+
+ # Get the data.
+ lx = max (x1, min (x2, int (xcen[istar] - routsky) + 1))
+ mx = max (x1, min (x2, int (xcen[istar] + routsky)))
+ ly = max (y1, min (y2, int (ycen[istar] - routsky) + 1))
+ my = max (y1, min (y2, int (ycen[istar] + routsky)))
+ line1 = ly
+ line2 = my
+ sub = dp_gst (dao, im, line1, line2, READ_ONLY, NO)
+ wgt = dp_gwt (dao, im, line1, line2, READ_ONLY, NO)
+
+ nsky = 0
+ psub = sub + (ly - DP_SYOFF(allstar)) * ncols
+ pwgt = wgt + (ly - DP_WYOFF(allstar)) * ncols
+ do j = ly, my {
+ dysq = (real (j) - ycen[istar]) ** 2
+ do i = lx, mx {
+ dx = (real (i) - xcen[istar])
+ rsq = dx ** 2 + dysq
+ if (rsq > rosq) {
+ if (dx > 0.0)
+ break
+ else
+ next
+ }
+ if (rsq < risq)
+ next
+ if (Memr[pwgt+i-1] <= badwt)
+ next
+ Memr[skyvals+nsky] = Memr[psub+i-1]
+ nsky = nsky + 1
+ }
+ psub = psub + ncols
+ pwgt = pwgt + ncols
+ }
+
+ # Compute the new sky value.
+ if (nsky > minnsky) {
+ call quick (Memr[skyvals], nsky, Memi[index], ier)
+ j = nint (0.2 * nsky)
+ dx = 0.0
+ do i = (nsky + 1) / 2 - j, (nsky / 2) + j + 1
+ dx = dx + Memr[skyvals+i-1]
+ sky[istar] = dx / real ((nsky / 2) + 2 * j + 2 -
+ (nsky + 1) / 2)
+ }
+ }
+
+
+ call mfree (skyvals, TY_REAL)
+ call mfree (index, TY_INT)
+ sub = dp_gst (dao, im, line1, line2, READ_ONLY, YES)
+ wgt = dp_gwt (dao, im, line1, line2, READ_ONLY, YES)
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:28:47 -0400