include include include 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