diff options
Diffstat (limited to 'noao/digiphot/daophot/group')
| -rw-r--r-- | noao/digiphot/daophot/group/dpgconfirm.x | 24 | ||||
| -rw-r--r-- | noao/digiphot/daophot/group/dpmkgroup.x | 484 | ||||
| -rw-r--r-- | noao/digiphot/daophot/group/dpsmpsf.x | 199 | ||||
| -rw-r--r-- | noao/digiphot/daophot/group/dpwrtgroup.x | 448 | ||||
| -rw-r--r-- | noao/digiphot/daophot/group/mkpkg | 18 | ||||
| -rw-r--r-- | noao/digiphot/daophot/group/t_group.x | 246 |
6 files changed, 1419 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/group/dpgconfirm.x b/noao/digiphot/daophot/group/dpgconfirm.x new file mode 100644 index 00000000..4dcb3003 --- /dev/null +++ b/noao/digiphot/daophot/group/dpgconfirm.x @@ -0,0 +1,24 @@ +# DP_GCONFIRM -- Confirm the critical GROUP task parameters. + +procedure dp_gconfirm (dao) + +pointer dao # pointer to the group structure + +begin + call printf ("\n") + + # Confirm the psf radius. + call dp_vpsfrad (dao) + + # Confirm the fitting radius. + call dp_vfitrad (dao) + + # Confirm the critical signal-to-noise ratio. + call dp_vcritsnratio (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/group/dpmkgroup.x b/noao/digiphot/daophot/group/dpmkgroup.x new file mode 100644 index 00000000..741fc8b0 --- /dev/null +++ b/noao/digiphot/daophot/group/dpmkgroup.x @@ -0,0 +1,484 @@ +include <mach.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +define EPS_OVERLAP (1.0e-10) + +# DP_MKGROUP -- Arrange the photometry results into natural groupings +# based on physical proximity and the signal-to-noise ratio. + +procedure dp_mkgroup (dao, im, grp) + +pointer dao # pointer to the daophot structure +pointer im # pointer to input image +int grp # the output file descriptor + +bool overlap +int i, maxgroup, curr_star, first_unknown, first_ingrp, nin_currgrp +int curr_point, crit_point, ier +pointer psffit, apsel, index, group_size, number +real fradius, critovlap, bright_mag, read_noise, fitradsq, critsep +real critsep_sq, xcurr, ycurr, dxcurr_frompsf, dycurr_frompsf, magcurr +real skycurr, magcrit, skycrit, dxcrit_frompsf, dycrit_frompsf +real deltax, deltay, radsq, rel_bright, radius, ratio, stand_err, dvdx, dvdy + +bool dp_gchkstar() +real dp_usepsf() + +begin + # Get the daophot pointers. + psffit = DP_PSFFIT(dao) + apsel = DP_APSEL(dao) + + # Return if the photometry file is empty. + if (DP_APNUM(apsel) <= 0) + return + + # Store the original fitting radius. + fradius = DP_FITRAD(dao) + DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao)) + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) + + # Get some working memory. + call malloc (index, DP_APNUM(apsel), TY_INT) + call malloc (group_size, DP_APNUM(apsel), TY_INT) + call malloc (number, DP_APNUM(apsel), TY_INT) + + # Initialize the group information. + call aclri (Memi[index], DP_APNUM(apsel)) + call aclri (Memi[group_size], DP_APNUM(apsel)) + call aclri (Memi[number], DP_APNUM(apsel)) + + # Check for INDEF results and fix them up.. + call dp_gpfix (dao, im, bright_mag) + + # Sort the stars into order of increasing y value. The star ID's, X, + # MAG and SKY values are still in the order in which they were read + # in. INDEX points to the proper value, i.e. Y (i) and X (index(i)). + + call quick (Memr[DP_APYCEN(apsel)], DP_APNUM(apsel), Memi[index], ier) + + # Bright_mag is the apparent magnitude of the brightest star + # in the input file. If this is INDEF then set to the PSFMAG. + + if (! IS_INDEFR(bright_mag)) + bright_mag = DAO_RELBRIGHT (psffit, bright_mag) + else + bright_mag = 1.0 + + # Smooth the PSF. + if ((DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit)) > 0) + call dp_smpsf (dao) + + # Define some important constants in terms of the daophot parameters. + + critovlap = DP_CRITSNRATIO(dao) * DP_CRITSNRATIO(dao) + read_noise = (DP_READNOISE (dao) / DP_PHOTADU(dao)) ** 2 + fitradsq = (DP_FITRAD(dao) + 1) * (DP_FITRAD(dao) + 1) + + ## Note that the definition of firadsq has been changed in this + ## version of group. This has been done so that the grouping + ## algorithm defaults to the one used by ALLSTAR in the case that + ## the critoverlap parameter is very large. + ## fitradsq = (2.0 * DP_FITRAD(dao)) * (2.0 * DP_FITRAD(dao)) + + # Define the critical separation. + critsep = DP_PSFRAD(dao) + DP_FITRAD(dao) + 1.0 + critsep_sq = critsep * critsep + + # Now we search the list for stars lying within one critical + # separation of each other. The stars are currently in a stack + # DP_APNUM(apsel) stars long. FIRST_INGRP points to the first star + # in the current group and starts out, of course, with the + # value CURR_STAR. CURR_STAR will point to the position in the stack + # occupied by the star which is currently the center of a circle + # of radius equal to the critical radius, within which we are + # looking for other stars. CURR_STAR starts out with a value of + # 1. FIRST_UNKNOWN points to the top position in the stack of the + # stars which have not yet been assigned to groups. FIRST_UNKNOWN + # starts out with the value 2. Each time through, the program goes + # down through the stack from FIRST_UNKNOWN to DP_APNUM(apsel) and + # looks for stars within the critical distance from the star at stack + # position CURR_STAR. When such a star is found, it changes places + # in the stack with the star at FIRST_UNKNOWN and FIRST_UNKNOWN is + # incremented by one. When the search has gotten to the last + # position in the stack (DP_APNUM(apsel)), the pointer CURR_STAR is + # incremented by one, and the search proceeds again from the new + # value of FIRST_UNKNOWN to DP_APNUM(apsel). If the pointer CURR_STAR + # catches up with the pointer FIRST_UNKNOWN, that means that the + # group currently being built up is complete. The number of stars in + # the newly-created group (the first star of which is at stack + # position FIRST_INGRP) is stored in array element GROUP_SIZE[FIRST_ + # INGRP]. Then a new group is started beginning with the star at the + # current position CURR_STAR, ( = FIRST_UNKNOWN for the moment), + # FIRST_UNKNOWN is incremented by 1, and the next group is built up + # as before. + + # Initialize. + maxgroup = 0 + curr_star = 1 + first_unknown = 1 + + # Loop through the list of stars. + while (curr_star <= DP_APNUM(apsel)) { + + # Initialize for the next group. + first_ingrp = curr_star + nin_currgrp = 1 + first_unknown = first_unknown + 1 + + # Begin defining a group. + while (curr_star < first_unknown) { + + # Get the parameters of the current star. + curr_point = Memi[index+curr_star-1] + xcurr = Memr[DP_APXCEN(apsel)+curr_point-1] + ycurr = Memr[DP_APYCEN(apsel)+curr_star-1] + call dp_wpsf (dao, im, xcurr, ycurr, dxcurr_frompsf, + dycurr_frompsf, 1) + dxcurr_frompsf = (dxcurr_frompsf - 1.0) / DP_PSFX(psffit) - 1.0 + dycurr_frompsf = (dycurr_frompsf - 1.0) / DP_PSFY(psffit) - 1.0 + magcurr = Memr[DP_APMAG(apsel)+curr_point-1] + skycurr = Memr[DP_APMSKY(apsel)+curr_point-1] + if (IS_INDEFR(magcurr)) + magcurr = bright_mag + else if (abs (DAO_MAGCHECK(psffit, magcurr)) > (MAX_EXPONENTR - + 1)) + magcurr = bright_mag + else + magcurr = DAO_RELBRIGHT (psffit, magcurr) + + # Go through the list of unassigned stars looking for stars + # within one critical distance of the current star. + + do i = first_unknown, DP_APNUM(apsel) { + + # Is this star within one critical separation of the + # current star ? + + overlap = dp_gchkstar (Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memi[index], i, xcurr, ycurr, + critsep_sq, deltax, deltay, radsq) + + # Break from the do loop if we are beyond the the critical + # separation. + + if (deltay > critsep) + break + if (! overlap) + next + + # Define the characteristics of the unknown star. + + crit_point = Memi[index+i-1] + magcrit = Memr[DP_APMAG(apsel)+crit_point-1] + skycrit = Memr[DP_APMSKY(apsel)+crit_point-1] + call dp_wpsf (dao, im, Memr[DP_APXCEN(apsel)+crit_point-1], + Memr[DP_APYCEN(apsel)+i-1], dxcrit_frompsf, + dycrit_frompsf, 1) + dxcrit_frompsf = (dxcrit_frompsf - 1.0) / DP_PSFX(psffit) - + 1.0 + dycrit_frompsf = (dycrit_frompsf - 1.0) / DP_PSFY(psffit) - + 1.0 + + # Check to see if stars overlap critically. + + if ((critovlap > EPS_OVERLAP) && (radsq > fitradsq)) { + + overlap = false + + # Rel_bright is the brightness of the star currently + # being tested relative to the star in the center. + + if (IS_INDEFR(magcrit)) + rel_bright = bright_mag + else if (abs (DAO_MAGCHECK (psffit, magcrit)) > + (MAX_EXPONENTR - 1)) + rel_bright = bright_mag + else + rel_bright = DAO_RELBRIGHT(psffit, magcrit) + + # Determine the point at which to evaluate the PSF. + + radius = sqrt (radsq) + ratio = (radius - (DP_FITRAD(dao) + 1)) / radius + deltax = ratio * deltax + deltay = ratio * deltay + + # Determine which star is brighter. + + if (magcurr > rel_bright) { + rel_bright = magcurr * + dp_usepsf (DP_PSFUNCTION(psffit), deltax, + deltay, DP_PSFHEIGHT(psffit), + Memr[DP_PSFPARS(psffit)], + Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), + dxcurr_frompsf, dycurr_frompsf, dvdx, dvdy) + } else { + rel_bright = rel_bright * + dp_usepsf (DP_PSFUNCTION(psffit), -deltax, + -deltay, DP_PSFHEIGHT(psffit), + Memr[DP_PSFPARS(psffit)], + Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), + dxcrit_frompsf, dycrit_frompsf, dvdx, dvdy) + } + + # Do the stars overlap ? + + if (! IS_INDEFR(skycurr) && ! IS_INDEFR(skycrit)) { + stand_err = read_noise + 0.5 * (skycurr + + skycrit) / DP_PHOTADU(dao) + if ((rel_bright * rel_bright) >= (stand_err * + critovlap)) + overlap = true + } + } + + # The two stars do overlap. Increment the pointers + # and move this star to the top of the stack. + + if (overlap) { + nin_currgrp = nin_currgrp + 1 + call dp_gswap2 (i, first_unknown, Memi[index], + Memr[DP_APYCEN(apsel)]) + first_unknown = first_unknown + 1 + } + } + + curr_star = curr_star + 1 + } + + # Check for maximum group size. + if (nin_currgrp > maxgroup) + maxgroup = nin_currgrp + + # Increment the number of groups versus size counter. Stars + # with undefined centers are always in a group with a membrship + # of 1 and are skipped. + + if (! IS_INDEFR(xcurr) && ! IS_INDEFR(ycurr)) + Memi[number+nin_currgrp-1] = Memi[number+nin_currgrp-1] + 1 + + # Define the group size. + Memi[group_size+first_ingrp-1] = nin_currgrp + } + + # Write out all the groups to the output file. + call dp_wrtgroup (dao, im, grp, Memi[number], Memi[index], + Memi[group_size], maxgroup) + + call mfree (number, TY_INT) + call mfree (index, TY_INT) + call mfree (group_size, TY_INT) + + # Restore the original fitting radius. + DP_FITRAD(dao) = fradius + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) +end + + +# DP_GSWAP2 -- Interchange two stars in the photometry list and then shift +# the list. + +procedure dp_gswap2 (star1, star2, index, y) + +int star1, star2 # the two star numbers to interchange +int index[ARB] # the index array +real y[ARB] # array of y positions + +int j, k, l, ihold +real yhold + +begin + yhold = y[star1] + ihold = index[star1] + + l = star1 + star2 + do j = star2, star1 - 1 { + k = l - j + y[k] = y[k-1] + index[k] = index [k-1] + } + + y[star2] = yhold + index[star2] = ihold +end + + +# DP_GCHKSTR -- Check to see if the unknown stars star is within one critical +# separation of the current star. + +bool procedure dp_gchkstar (x, y, index, i, xcurr, ycurr, crit, deltax, + deltay, radsq) + +real x[ARB] # array of x positions +real y[ARB] # array of y positions +int index[ARB] # index array from the quick sort +int i # position of test star in stack +real xcurr # x position of current star +real ycurr # y position of current star +real crit # critical radius squared +real deltax # output difference in x +real deltay # output difference in y +real radsq # separation squared + +begin + # Initialize the deltas. + deltax = MAX_REAL + deltay = MAX_REAL + + # Check to see if star is within a critical radius. Reject the + # star if any of the positions are INDEF. + if (IS_INDEFR(xcurr) || IS_INDEFR(ycurr)) { + return (false) + } else if (IS_INDEFR(y[i]) || IS_INDEFR(x[index[i]])) { + return (false) + } else { + deltay = y[i] - ycurr + deltax = x[index[i]] - xcurr + radsq = deltax * deltax + deltay * deltay + if (radsq <= crit) + return (true) + else + return (false) + } +end + + +# DP_GPFIX -- Check for INDEF photometry and get estimate of magnitude. + +procedure dp_gpfix (dao, im, bright_mag) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +real bright_mag # apparent magnitude of the brightest star + +int i, lowx, lowy, nxpix, nypix +pointer apsel, subim +real fitrad, fitradsq, mingdata, maxgdata, x, y, sky, mag + +bool fp_equalr() +pointer dp_gsubrast() +real dp_gaccum() + +begin + # Get pointers to the daophot structures. + apsel = DP_APSEL (dao) + + # Initialize. + fitrad = DP_FITRAD(dao) + fitradsq = fitrad * fitrad + 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) + bright_mag = MAX_REAL + + # Get a magnitude estimate for stars with INDEF magnitudes by summing + # all of the pixels within one fitting radius and scaling with respect + # to the PSFMAG. + + do i = 1, DP_APNUM(apsel) { + + # Check for undefined centers. + x = Memr[DP_APXCEN(apsel)+i-1] + y = Memr[DP_APYCEN(apsel)+i-1] + if (IS_INDEFR(x) || IS_INDEFR(y)) + next + + # Check for an undefined sky. + sky = Memr[DP_APMSKY(apsel)+i-1] + if (IS_INDEFR(sky)) + next + + # Correct the magnitudes. + mag = Memr[DP_APMAG(apsel)+i-1] + if (IS_INDEFR(mag)) { + + # Get subraster around star position. If the subraster + # cannot be extracted leave the magnitude as INDEF. + subim = dp_gsubrast (im, x, y, fitrad, lowx, lowy, nxpix, + nypix) + if (subim == NULL) + next + + # Estimate the value of the PSF. + mag = dp_gaccum (dao, Memr[subim], nxpix, nypix, lowx, lowy, + x, y, sky, fitradsq, mingdata, maxgdata) + + if (! IS_INDEFR(mag)) + Memr[DP_APMAG(apsel)+i-1] = mag + + + } else if (mag < bright_mag) + bright_mag = mag + } + + if (fp_equalr (bright_mag, MAX_REAL)) + bright_mag = INDEFR +end + + +# DP_GACCUM -- Accumulate the model counts. + +real procedure dp_gaccum (dao, subim, nxpix, nypix, lowx, lowy, x, y, sky, + fitradsq, mingdata, maxgdata) + +pointer dao # pointer to the daophot structure +real subim[nxpix,nypix] # the input data subraster +int nxpix, nypix # the dimensions of the data subraster +int lowx, lowy # the lower left corner of the subraster +real x, y # the coordinates of the star +real sky # the sky value of the star +real fitradsq # the fitting radius of the star +real mingdata, maxgdata # the minimum and maximum good data values + +int k, j +pointer psffit +real mag, dxfrom_psf, dyfrom_psf, numer, denom, dx, dy, radsq, dvdx, dvdy +real weight, value +real dp_usepsf() + +begin + psffit = DP_PSFFIT(dao) + + # Compute the distance from the psf star. + dxfrom_psf = (x - 1.0) / DP_PSFX(psffit) - 1.0 + dyfrom_psf = (y - 1.0) / DP_PSFY(psffit) - 1.0 + + denom = 0. + numer = 0. + mag = INDEFR + do k = 1, nypix { + do j = 1, nxpix { + + dx = real (lowx + j - 1) - x + dy = real (lowy + k - 1) - y + radsq = dx * dx + dy * dy + if ((radsq >= fitradsq) || (subim[j,k] < mingdata) || + (subim[j,k] > maxgdata)) + next + + value = 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) + radsq = radsq / fitradsq + weight = 5. / (5. + radsq / (1.0 - radsq)) + numer = numer + weight * value * (subim[j,k] - sky) + denom = denom + weight * value * value + } + } + + if (denom > 0.0 && numer > 0.0) + mag = DP_PSFMAG(psffit) - 2.5 * log10 (numer / denom) + + return (mag) +end diff --git a/noao/digiphot/daophot/group/dpsmpsf.x b/noao/digiphot/daophot/group/dpsmpsf.x new file mode 100644 index 00000000..96cdb324 --- /dev/null +++ b/noao/digiphot/daophot/group/dpsmpsf.x @@ -0,0 +1,199 @@ +include <mach.h> +include "../lib/daophotdef.h" + +define NSEC 4 +define IRMIN 4 + +# DP_SMPSF -- Smooth the psf before grouping. + +procedure dp_smpsf (dao) + +pointer dao # pointer to the daophot strucuture + +int k, icenter, irmax, nexpand +pointer psffit, sum, high, low, n +real rmax + +begin + # Get some pointers. + psffit = DP_PSFFIT(dao) + + # Get some constants. + icenter = (DP_PSFSIZE(psffit) + 1) / 2 + rmax = .7071068 * real (DP_PSFSIZE(psffit) - 1) + irmax = int (rmax + 1.0e-5) + nexpand = DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit) + + # Allocate working memory. + call malloc (sum, NSEC * irmax, TY_REAL) + call malloc (high, NSEC * irmax, TY_REAL) + call malloc (low, NSEC * irmax, TY_REAL) + call malloc (n, NSEC * irmax, TY_INT) + + # Do the smoothing. + do k = 1, nexpand { + + # Initialize. + call aclrr (Memr[sum], NSEC * irmax) + call amovkr (-MAX_REAL, Memr[high], NSEC * irmax) + call amovkr (MAX_REAL, Memr[low], NSEC * irmax) + call aclri (Memi[n], NSEC * irmax) + + # Acumulate. + call dp_smaccum (Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_PSFSIZE(psffit), Memr[sum], Memr[low], Memr[high], Memi[n], + NSEC, IRMIN, icenter, rmax, k) + + # Normalize. + call dp_smnorm (Memr[sum], Memr[low], Memr[high], Memi[n], NSEC, + IRMIN, irmax) + + # Smooth. + call dp_smo (Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_PSFSIZE(psffit), Memr[sum], NSEC, IRMIN, icenter, rmax, k) + } + + call mfree (sum, TY_REAL) + call mfree (low, TY_REAL) + call mfree (high, TY_REAL) + call mfree (n, TY_INT) +end + + +# DP_SMACCUM -- Accumulate the sums and limits + +procedure dp_smaccum (psflut, nxpsf, nypsf, sum, low, high, n, nsec, irmin, + icenter, rmax, k) + +real psflut[nxpsf,nypsf,ARB] # the psf lookup table +int nxpsf, nypsf # size of the psf lookup table +real sum[nsec,ARB] # array of sums +real low[nsec,ARB] # array of low values +real high[nsec,ARB] # array of high values +int n[nsec,ARB] # array of number of points +int nsec # dimension of sum arrays +int irmin # number of sums +int icenter # center of the array +real rmax # max radius +int k # third dimension array index + +int i, j, idx, idy, is, ir +real dxsq, dysq, r +int dp_isctr() + +begin + do j = 1, nypsf { + idy = j - icenter + dysq = idy ** 2 + do i = 1, nxpsf { + idx = i - icenter + dxsq = idx ** 2 + r = sqrt (dxsq + dysq) + if (r > rmax) + next + ir = int (r + 1.0e-5) + if (ir < irmin) + next + is = dp_isctr (idx, idy) + sum[is,ir] = sum[is,ir] + psflut[i,j,k] + if (psflut[i,j,k] > high[is,ir]) + high[is,ir] = psflut[i,j,k] + if (psflut[i,j,k] < low[is,ir]) + low[is,ir] = psflut[i,j,k] + n[is,ir] = n[is,ir] + 1 + } + } +end + + +# DP_SMNORM -- Normalize the sum + +procedure dp_smnorm (sum, low, high, n, nsec, irmin, irmax) + +real sum[nsec,ARB] # array of sums +real low[nsec,ARB] # array of low values +real high[nsec,ARB] # array of high values +int n[nsec,ARB] # array of counter +int nsec # array dimension +int irmin # radius index +int irmax # maximum radius index + +int ir, is + +begin + do ir = irmin, irmax { + do is = 1, nsec { + if (n[is,ir] > 2) + sum[is,ir] = (sum[is,ir] - high[is,ir] - low[is,ir]) / + (n[is,ir] - 2) + } + } +end + + +# DP_SMO -- Do the actual smoothing. + +procedure dp_smo (psflut, nxpsf, nypsf, sum, nrec, irmin, icenter, rmax, k) + +real psflut[nxpsf,nypsf,ARB] # the lookup table +int nxpsf, nypsf # size of the psf lookup table +real sum[nrec,ARB] # array of sums +int nrec # dimension of sum array +int irmin # min radius index +int icenter # index of center +real rmax # maximum radius +int k # index of third dimension + +int i, j, idx, idy, ir, is +real dysq, r +int dp_isctr() + +begin + do j = 1, nypsf { + idy = j - icenter + dysq = idy ** 2 + do i = 1, nxpsf { + idx = i - icenter + r = sqrt (real (idx ** 2) + dysq) + if (r > rmax) + next + ir = int (r + 1.0e-5) + if (ir < irmin) + next + is = dp_isctr (idx, idy) + psflut[i,j,k] = sum[is,ir] + } + } +end + + +# DP_ISCTR -- Convert an index pair into a numbered sector from 1 to 4. + +int procedure dp_isctr (i, j) + +int i # first index +int j # second index + +int isctr + +begin + if (i > 0) { + isctr = 1 + } else if (i < 0) { + isctr = 3 + } else { + if (j <= 0) + isctr = 1 + else + isctr = 3 + } + + if (j > 0) { + isctr = isctr + 1 + } else if (j == 0) { + if (i > 0) + isctr = 2 + } + + return (isctr) +end diff --git a/noao/digiphot/daophot/group/dpwrtgroup.x b/noao/digiphot/daophot/group/dpwrtgroup.x new file mode 100644 index 00000000..97b9a714 --- /dev/null +++ b/noao/digiphot/daophot/group/dpwrtgroup.x @@ -0,0 +1,448 @@ +include <mach.h> +include <time.h> +include <tbset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +define NCOLUMN 6 + +# DP_TGNEWGRP -- Create a new GROUP output ST table. + +procedure dp_tgnewgrp (tp, colpoint) + +pointer tp # pointer to outpu ST table +pointer colpoint[ARB] # array of pointers to columns + +pointer sp, colnames, colunits, colformat, col_dtype, col_len + +begin + # Allocate space for table definition. + call smark (sp) + call salloc (colnames, NCOLUMN * (SZ_COLNAME + 1), TY_CHAR) + call salloc (colunits, NCOLUMN * (SZ_COLUNITS + 1), TY_CHAR) + call salloc (colformat, NCOLUMN * (SZ_COLFMT + 1), TY_CHAR) + call salloc (col_dtype, NCOLUMN, TY_INT) + call salloc (col_len, NCOLUMN, TY_INT) + + # Set up the column definitions. + call strcpy (GROUP, Memc[colnames], SZ_COLNAME) + call strcpy (ID, Memc[colnames+SZ_COLNAME+1], SZ_COLNAME) + call strcpy (XCENTER, Memc[colnames+2*SZ_COLNAME+2], SZ_COLNAME) + call strcpy (YCENTER, Memc[colnames+3*SZ_COLNAME+3], SZ_COLNAME) + call strcpy (MAG, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME) + call strcpy (SKY, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME) + + # Set up the format definitions. + call strcpy ("%6d", Memc[colformat], SZ_COLFMT) + call strcpy ("%6d", Memc[colformat+SZ_COLFMT+1], SZ_COLFMT) + call strcpy ("%10.3f", Memc[colformat+2*SZ_COLFMT+2], SZ_COLFMT) + call strcpy ("%10.3f", Memc[colformat+3*SZ_COLFMT+3], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT) + call strcpy ("%15.7g", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT) + + # Set up the unit definitions. + call strcpy ("##", Memc[colunits], SZ_COLUNITS) + call strcpy ("##", Memc[colunits+SZ_COLUNITS+1], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+2*SZ_COLUNITS+2], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+3*SZ_COLUNITS+3], SZ_COLUNITS) + call strcpy ("MAGNITUDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS) + call strcpy ("ADC", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS) + + # Set up the data type definitions. + Memi[col_dtype] = TY_INT + Memi[col_dtype+1] = TY_INT + Memi[col_dtype+2] = TY_REAL + Memi[col_dtype+3] = TY_REAL + Memi[col_dtype+4] = TY_REAL + Memi[col_dtype+5] = TY_REAL + + # Define the column lengths. + Memi[col_len] = 1 + Memi[col_len+1] = 1 + Memi[col_len+2] = 1 + Memi[col_len+3] = 1 + Memi[col_len+4] = 1 + Memi[col_len+5] = 1 + + # Define and create the table. + call tbcdef (tp, colpoint, Memc[colnames], Memc[colunits], + Memc[colformat], Memi[col_dtype], Memi[col_len], NCOLUMN) + call tbtcre (tp) + + call sfree (sp) +end + + +# DP_XGGRPPARS -- Write out the parameters to the header of the GROUP text +# output file. + +procedure dp_xggrppars (dao, tp) + +pointer dao # pointer to the DAOPHOT structure +int tp # the output file descriptor + +pointer sp, outstr, date, time, psffit +int envfind() + +begin + # Allocate working space. + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + call salloc (date, SZ_DATE, TY_CHAR) + call salloc (time, SZ_DATE, TY_CHAR) + + psffit = DP_PSFFIT(dao) + + # Write the id. + if (envfind ("version", Memc[outstr], SZ_LINE) <= 0) + call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE) + call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE) + call dp_sparam (tp, "IRAF", Memc[outstr], "version", "") + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call dp_sparam (tp, "USER", Memc[outstr], "name", "") + call gethost (Memc[outstr], SZ_LINE) + call dp_sparam (tp, "HOST", Memc[outstr], "computer", "") + call dp_date (Memc[date], Memc[time], SZ_DATE) + call dp_sparam (tp, "DATE", Memc[date], "yyyy-mm-dd", "") + call dp_sparam (tp, "TIME", Memc[time], "hh:mm:ss", "") + call dp_sparam (tp, "PACKAGE", "daophot", "name", "") + call dp_sparam (tp, "TASK", "group", "name", "") + + # Write the file name parameters. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "IMAGE", Memc[outstr], "imagename", "") + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "PHOTFILE", Memc[outstr], "filename", "") + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "PSFIMAGE", Memc[outstr], "imagename", "") + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "GRPFILE", Memc[outstr], "filename", "") + + # Write out relevant data parameters. + call dp_rparam (tp, "SCALE", DP_SCALE(dao), "units/pix", "") + call dp_rparam (tp, "DATAMIN", DP_MINGDATA(dao), "counts", "") + call dp_rparam (tp, "DATAMAX", DP_MAXGDATA(dao), "counts", "") + call dp_rparam (tp, "GAIN", DP_PHOTADU(dao), "number", "") + call dp_rparam (tp, "READNOISE", DP_READNOISE(dao), "electrons", "") + + # Write out the observing parameters. + call dp_sparam (tp, "OTIME", DP_OTIME(dao), "timeunit", "") + call dp_rparam (tp, "XAIRMASS", DP_XAIRMASS(dao), "number", "") + call dp_sparam (tp, "IFILTER", DP_IFILTER(dao), "filter", "") + + # Write out the daophot parameters. + call dp_rparam (tp, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", "") + call dp_rparam (tp, "FITRAD", DP_SFITRAD(dao), "scaleunit", "") + call dp_rparam (tp, "PSFMAG", DP_PSFMAG(psffit), "magnitude", "") + call dp_rparam (tp, "CRITSNRATIO", DP_CRITSNRATIO(dao), "sigma", "") + call dp_iparam (tp, "MAXGROUP", DP_MAXGROUP(dao), "number", "") + #call dp_iparam (tp, "MAXNSTAR", DP_MAXNSTAR(dao), "number", "") + + # Write out the group size parameters. + # call dp_iparam (tp, "MINSZGROUP", 1, "number", "") + # call dp_iparam (tp, "MAXSZGROUP", MAX_INT, "number", "") + + call sfree(sp) +end + + +# DP_TGGRPPARS -- Write out the parameters to the header of the GROUP output +# ST table file. + +procedure dp_tggrppars (dao, tp) + +pointer dao # pointer to the DAOPHOT structure +pointer tp # pointer to the output table + +pointer sp, outstr, date, time, psffit +int envfind() + +begin + # Allocate working space. + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + call salloc (date, SZ_DATE, TY_CHAR) + call salloc (time, SZ_DATE, TY_CHAR) + + psffit = DP_PSFFIT(dao) + + # Write the id. + if (envfind ("version", Memc[outstr], SZ_LINE) <= 0) + call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE) + call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE) + call tbhadt (tp, "IRAF", Memc[outstr]) + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call tbhadt (tp, "USER", Memc[outstr]) + call gethost (Memc[outstr], SZ_LINE) + call tbhadt (tp, "HOST", Memc[outstr]) + call dp_date (Memc[date], Memc[time], SZ_DATE) + call tbhadt (tp, "DATE", Memc[date]) + call tbhadt (tp, "TIME", Memc[time]) + call tbhadt (tp, "PACKAGE", "daophot") + call tbhadt (tp, "TASK", "group") + + # Write the file name parameters. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "IMAGE", Memc[outstr]) + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PHOTFILE", Memc[outstr]) + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PSFIMAGE", Memc[outstr]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "GRPFILE", Memc[outstr]) + + # Write out relevant data parameters. + call tbhadr (tp, "SCALE", DP_SCALE(dao)) + call tbhadr (tp, "DATAMIN", DP_MINGDATA(dao)) + call tbhadr (tp, "DATAMAX", DP_MAXGDATA(dao)) + call tbhadr (tp, "GAIN", DP_PHOTADU(dao)) + call tbhadr (tp, "READNOISE", DP_READNOISE(dao)) + + # Write out the observing parameters. + call tbhadt (tp, "OTIME", DP_OTIME(dao)) + call tbhadr (tp, "XAIRMASS", DP_XAIRMASS(dao)) + call tbhadt (tp, "IFILTER", DP_IFILTER(dao)) + + # Write out the daophot parameters. + call tbhadr (tp, "PSFRAD", DP_SPSFRAD(dao)) + call tbhadr (tp, "FITRAD", DP_SFITRAD(dao)) + call tbhadr (tp, "PSFMAG", DP_PSFMAG(psffit)) + call tbhadr (tp, "CRITSNRATIO", DP_CRITSNRATIO(dao)) + call tbhadi (tp, "MAXGROUP", DP_MAXGROUP(dao)) + #call tbhadi (tp, "MAXNSTAR", DP_MAXNSTAR(dao)) + + # call tbhadi (tp, "MINSZGROUP", 1) + # call tbhadi (tp, "MAXSZGROUP", MAX_INT) + + call sfree(sp) +end + + +# DP_WRTGROUP -- Write out each group into a text file or an ST table. + +procedure dp_wrtgroup (dao, im, grp, number, index, group_size, maxgroup) + +pointer dao # pointer to daophot structure +pointer im # the input image descriptor +int grp # the output file descriptor +int number[ARB] # number in group of each size +int index[ARB] # index to results +int group_size[ARB] # size of groups +int maxgroup # maximum group size + +begin + if (DP_TEXT(dao) == YES) + call dp_xwrtgroup (dao, im, grp, number, index, group_size, + maxgroup) + else + call dp_twrtgroup (dao, im, grp, number, index, group_size, + maxgroup) +end + + +define GR_NAMESTR "#N%4tGROUP%10tID%16tXCENTER%26tYCENTER%36tMAG%48tMSKY\ +%80t\\\n" +define GR_UNITSTR "#U%4t##%10t##%16tpixels%26tpixels%36tmagnitudes%48tcounts\ +%80t\\\n" +define GR_FORMATSTR "#F%4t%%-9d%10t%%-6d%16t%%-10.3f%26t%%-10.3f%36t%%-12.3f\ +%48t%%-15.7g%80t \n" +define GR_DATASTR "%-9d%10t%-6d%16t%-10.3f%26t%-10.3f%36t%-12.3f%48t%-15.7g\ +%80t \n" + + +# DP_XWRTGROUP -- Write each group into the GROUP output ST table. + +procedure dp_xwrtgroup (dao, im, grp, number, index, group_size, maxgroup) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int grp # the output file descriptor +int number[ARB] # number in group of each size +int index[ARB] # index to results +int group_size[ARB] # size of groups +int maxgroup # maximum group size + +int i, j, k, id, ngroup, nstars, first_ingrp +pointer apsel +real x, y, mag, sky + +begin + # Get the daophot pointer. + apsel = DP_APSEL(dao) + + # Print results to the standard output. + if (DP_VERBOSE(dao) == YES) { + call printf (" Size of Number of\n") + call printf (" group groups\n\n") + do i = 1, maxgroup { + if (number[i] <= 0) + next + call printf ("%8d %9d\n") + call pargi (i) + call pargi (number[i]) + } + } + + # Add header parameters to the table. + call dp_xggrppars (dao, grp) + + # Write out the banner. + call fprintf (grp, "#\n") + call fprintf (grp, GR_NAMESTR) + call fprintf (grp, GR_UNITSTR) + call fprintf (grp, GR_FORMATSTR) + call fprintf (grp, "#\n") + + # Write out each group. + ngroup = 1 + first_ingrp = 1 + while (first_ingrp <= DP_APNUM(apsel)) { + + do j = first_ingrp, first_ingrp + group_size[first_ingrp] - 1 { + + # Test the center. + k = index[j] + x = Memr[DP_APXCEN(apsel)+k-1] + y = Memr[DP_APYCEN(apsel)+j-1] + if (IS_INDEFR(x) || IS_INDEFR(y)) + break + call dp_wout (dao, im, x, y, x, y, 1) + + # Get the rest of the numbers. + id = Memi[DP_APID(apsel)+k-1] + mag = Memr[DP_APMAG(apsel)+k-1] + sky = Memr[DP_APMSKY(apsel)+k-1] + + # Write the results. + call fprintf (grp, GR_DATASTR) + call pargi (ngroup) + call pargi (id) + call pargr (x) + call pargr (y) + call pargr (mag) + call pargr (sky) + } + + ngroup = ngroup + 1 + first_ingrp = first_ingrp + group_size[first_ingrp] + } + + # Compute the number of groups and the number of stars. + ngroup = 0 + nstars = 0 + do i = 1, maxgroup { + if (number[i] <= 0) + next + nstars = nstars + i * number[i] + ngroup = ngroup + number[i] + } + + if (DP_VERBOSE(dao) == YES) { + call printf ("\nTotal of %d stars in %d groups\n") + call pargi (nstars) + call pargi (ngroup) + } +end + + +# DP_TWRTGROUP -- Write each group into the GROUP output text file. + +procedure dp_twrtgroup (dao, im, grp, number, index, group_size, maxgroup) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +pointer grp # pointer to group output file +int number[ARB] # number in group of each size +int index[ARB] # index to results +int group_size[ARB] # size of groups +int maxgroup # maximum group size + +int row, first_ingrp, ngroup, nstars, i, j, k, id +pointer apsel, sp, colpoint +real x, y, mag, sky + +begin + # Get the daophot pointer. + apsel = DP_APSEL(dao) + + # Allocate space for the column pointers. + call smark( sp) + call salloc (colpoint, NCOLUMN, TY_INT) + + # Get the necessary info to create the ST table. + call dp_tgnewgrp (grp, Memi[colpoint]) + + # Add header parameters to the ST table. + call dp_tggrppars (dao, grp) + + # Optionally print results to the standard output. + if (DP_VERBOSE(dao) == YES) { + call printf (" Size of Number of\n") + call printf (" group groups\n\n") + do i = 1, maxgroup { + if (number[i] <= 0) + next + call printf (" %d %d\n") + call pargi (i) + call pargi (number[i]) + } + } + + # Initialize for writing. + ngroup = 1 + row = 0 + + # Initialize for reading. + first_ingrp = 1 + + # Write out the data for all the groups. + while (first_ingrp <= DP_APNUM(apsel)) { + + do j = first_ingrp, first_ingrp + group_size[first_ingrp] - 1 { + + # Test the center. + k = index[j] + x = Memr[DP_APXCEN(apsel)+k-1] + y = Memr[DP_APYCEN(apsel)+j-1] + if (IS_INDEFR(x) || IS_INDEFR(y)) + break + call dp_wout (dao, im, x, y, x, y, 1) + + # Get the rest of the values. + id = Memi[DP_APID(apsel)+k-1] + mag = Memr[DP_APMAG(apsel)+k-1] + sky = Memr[DP_APMSKY(apsel)+k-1] + + # Copy the values to the correct table row. + row = row + 1 + call tbrpti (grp, Memi[colpoint], ngroup, 1, row) + call tbrpti (grp, Memi[colpoint+1], id, 1, row) + call tbrptr (grp, Memi[colpoint+2], x, 1, row) + call tbrptr (grp, Memi[colpoint+3], y, 1, row) + call tbrptr (grp, Memi[colpoint+4], mag, 1, row) + call tbrptr (grp, Memi[colpoint+5], sky, 1, row) + } + + ngroup = ngroup + 1 + first_ingrp = first_ingrp + group_size[first_ingrp] + } + + # Compute the number of groups and the number of stars. + ngroup = 0 + nstars = 0 + do i = 1, maxgroup { + if (number[i] <= 0) + next + nstars = nstars + i * number[i] + ngroup = ngroup + number[i] + } + + # Optionally print out a summary of the results. + if (DP_VERBOSE(dao) == YES) { + call printf ("\nTotal of %d stars in %d groups\n") + call pargi (nstars) + call pargi (ngroup) + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/group/mkpkg b/noao/digiphot/daophot/group/mkpkg new file mode 100644 index 00000000..35173772 --- /dev/null +++ b/noao/digiphot/daophot/group/mkpkg @@ -0,0 +1,18 @@ +# GROUP task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + dpmkgroup.x <mach.h> ../lib/apseldef.h \ + ../lib/daophotdef.h + dpsmpsf.x <mach.h> ../lib/daophotdef.h + dpgconfirm.x + dpwrtgroup.x <mach.h> <time.h> \ + <tbset.h> ../lib/daophotdef.h \ + ../lib/apseldef.h + t_group.x <fset.h> <imhdr.h> \ + ../lib/daophotdef.h + ; diff --git a/noao/digiphot/daophot/group/t_group.x b/noao/digiphot/daophot/group/t_group.x new file mode 100644 index 00000000..5596c0b3 --- /dev/null +++ b/noao/digiphot/daophot/group/t_group.x @@ -0,0 +1,246 @@ +include <fset.h> +include <imhdr.h> +include "../lib/daophotdef.h" + +# T_GROUP -- Procedure to divide the stars in a photometry table into +# natural groups based on the magnitude level at which they overlap. + +procedure t_group () + +pointer image # name of the image +pointer apfile # aperture photometry file +pointer psfimage # name of the output PSF +pointer groupfile # output group table + +pointer sp, im, dao, outfname, str +int apd, root, cache, verbose, verify, update, grp, tp, wcs +int imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist +int req_size, old_size, buf_size, memstat +bool ap_text + +pointer immap(), tbtopn() +int access(), fnldir(), strlen(), strncmp(), fstati(), btoi() +int imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb(), fntgfnb() +int open(), 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 memory. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (apfile, SZ_FNAME, TY_CHAR) + call salloc (psfimage, SZ_FNAME, TY_CHAR) + call salloc (groupfile, 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[apfile], SZ_FNAME) + call clgstr ("psfimage", Memc[psfimage], SZ_FNAME) + call clgstr ("groupfile", Memc[groupfile], SZ_FNAME) + verbose = btoi (clgetb ("verbose")) + 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[apfile], NO) + lalist = fntlenb (alist) + pimlist = imtopen (Memc[psfimage]) + lpimlist = imtlen (pimlist) + olist = fntopnb (Memc[groupfile], NO) + lolist = fntlenb (olist) + + # Test that the lengths of the photometry file, psf image, and output + # file lists are the same as the length of the input image list. + + if ((limlist != lalist) && (strncmp (Memc[apfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call sfree (sp) + call error (0, + "Incompatible image and photometry file list lengths") + } + + if ((limlist != lpimlist) && (strncmp (Memc[psfimage], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call sfree (sp) + call error (0, + "Incompatible image and psf file list lengths") + } + + if ((limlist != lolist) && (strncmp (Memc[groupfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call sfree (sp) + call error (0, + "Incompatible image and group file list lengths") + } + + # Initialize the daophot structure and get the pset parameters. + call dp_gppars (dao) + call dp_seti (dao, VERBOSE, verbose) + + # Optionally verify and update the parameters. + if (verify == YES) { + call dp_gconfirm (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) + + # Open the PSF structure. + call dp_fitsetup (dao) + + # Open the photometry list structure. + call dp_apselsetup (dao) + + # Loop over the image list. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open input image and grab some header parameters. + im = immap (Memc[image], READ_ONLY, 0) + call dp_imkeys (dao, im) + call dp_sets (dao, INIMAGE, Memc[image]) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = dp_memstat (cache, req_size, old_size) + if (memstat == YES) + call dp_pcache (im, INDEFI, buf_size) + + # Open input photometry list and read in the photometry. + if (fntgfnb (alist, Memc[apfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[apfile], SZ_FNAME) + root = fnldir (Memc[apfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[apfile+root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[apfile])) + call dp_inname (Memc[image], Memc[outfname], "mag", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[apfile], Memc[outfname], SZ_FNAME) + ap_text = itob (access (Memc[outfname], 0, TEXT_FILE)) + if (ap_text) + apd = open (Memc[outfname], READ_ONLY, TEXT_FILE) + else + apd = tbtopn (Memc[outfname], READ_ONLY, 0) + call dp_wgetapert (dao, im, apd, DP_MAXNSTAR(dao), ap_text) + call dp_sets (dao, INPHOTFILE, Memc[outfname]) + + # Read 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) + tp = immap (Memc[outfname], READ_ONLY, 0) + call dp_readpsf (dao, tp) + call dp_sets (dao, PSFIMAGE, Memc[outfname]) + + # Open output GROUP file. If the output is "default", dir$default + # or a directory specification then the extension "grp" is added to + # the image name and a suitable version number is appended to the + # output name. + + if (fntgfnb (olist, Memc[groupfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[groupfile], SZ_FNAME) + root = fnldir (Memc[groupfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[groupfile + root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[groupfile])) + call dp_outname (Memc[image], Memc[outfname], "grp", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[groupfile], Memc[outfname], SZ_FNAME) + if (DP_TEXT(dao) == YES) + grp = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + grp = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTPHOTFILE, Memc[outfname]) + + # Now go and group the stars. + call dp_mkgroup (dao, im, grp) + + # Close up the input image. + call imunmap (im) + + # Close up the photometry file. + if (ap_text) + call close (apd) + else + call tbtclo (apd) + + # Close up the PSF image. + call imunmap (tp) + + # Close up the group table. + if (DP_TEXT(dao) == YES) + call close (grp) + else + call tbtclo (grp) + + # Uncache memory. + call fixmem (old_size) + + } + + # Close the image/file lists. + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + + # Free the photometry structure + call dp_apclose (dao) + + # Free the PSF structure. + call dp_fitclose (dao) + + # Free the daophot structure. + call dp_free (dao) + + call sfree(sp) +end |