diff options
Diffstat (limited to 'noao/digiphot/daophot/nstar')
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpggroup.x | 386 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpmemnstar.x | 158 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpnconfirm.x | 34 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpnstar.x | 355 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpnstarfit.x | 1383 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/dpntwrite.x | 600 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/mkpkg | 24 | ||||
| -rw-r--r-- | noao/digiphot/daophot/nstar/t_nstar.x | 308 |
8 files changed, 3248 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/nstar/dpggroup.x b/noao/digiphot/daophot/nstar/dpggroup.x new file mode 100644 index 00000000..fa180c17 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpggroup.x @@ -0,0 +1,386 @@ +include "../../lib/ptkeysdef.h" +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +# DP_GNSTPSF -- Procedure to initialize for reading the group file fields from +# a photometry text file . The group file fields are ID, GROUP, X, Y, MAG, ERR, +# and SKY. + +procedure dp_gnstpsf (fields, sel_fields, max_nfields) + +int fields[ARB] # array of selected fields +char sel_fields[ARB] # names of selected containing fields +int max_nfields # maximum number of fields selected + +int i +int strlen() + +begin + # Initialize the fields string. + sel_fields[1] = EOS + + # Encode the fields. + do i = 1, max_nfields { + switch (fields[i]) { + case DP_PAPID: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (ID) + case DP_PAPXCEN: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (XCENTER) + case DP_PAPYCEN: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (YCENTER) + case DP_PAPSKY: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (SKY) + case DP_PAPMAG1: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (MAG) + case DP_PAPGROUP: + call sprintf (sel_fields[strlen(sel_fields)+1], SZ_LINE, "%s ") + call pargstr (GROUP) + } + } + + # Backspace over the terminating blank character. + if (sel_fields[1] != EOS) + sel_fields[strlen(sel_fields)] = EOS +end + + +# DP_TNSTINIT -- Procedure to initialize for reading the group file fields from +# a photometry table. The group file fields are ID, GROUP, X, Y, MAG, ERR, +# and SKY. + +procedure dp_tnstinit (tp, colpoint) + +pointer tp # the table descriptor +pointer colpoint[ARB] # the column descriptor + +begin + # Get the id. + # First the ID + call tbcfnd (tp, ID, colpoint[1], 1) + if (colpoint[1] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (ID) + } + + # Get the x position. + call tbcfnd (tp, XCENTER, colpoint[2], 1) + if (colpoint[2] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (XCENTER) + } + + # Get the y position. + call tbcfnd (tp, YCENTER, colpoint[3], 1) + if (colpoint[3] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (YCENTER) + } + + # Get the magnitude. + call tbcfnd (tp, MAG, colpoint[4], 1) + if (colpoint[4] == NULL) # No column + call tbcfnd (tp, APMAG, colpoint[4], 1) + if (colpoint[4] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (APMAG) + } + + # Get the sky. + call tbcfnd (tp, SKY, colpoint[5], 1) + if (colpoint[5] == NULL) + call tbcfnd (tp, SKY, colpoint[5], 1) + if (colpoint[5] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (SKY) + } + + # Get the group number. + call tbcfnd (tp, GROUP, colpoint[6], 1) + if (colpoint[6] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (GROUP) + } +end + + +# DP_GGROUP -- Read in a single group. + +int procedure dp_ggroup (dao, tp, key, fields, indices, colpoint, max_row, + max_group, in_record, curr_group) + +pointer dao # pointer to daophot structure +int tp # input file/table descriptor +pointer key # pointer to text database structure +char fields[ARB] # nstar fields to be read +int indices[ARB] # array of text file field pointers +int colpoint[ARB] # array of column pointers +int max_row # number of rows in table +int max_group # maximum group size +int in_record # pointer to current input record +int curr_group # current group number + +bool nullflag +int istar, group, buf_size +pointer apsel +int dp_nstsel() + +begin + # If the current input record is set to zero we are at EOF. In_record + # is initialized to one on entry to this routine. + + if (in_record == 0) + return (0) + + # Get the next group number. Note that the last star read on the + # previous call is the first star in the new group. + + apsel = DP_APSEL(dao) + if (in_record == 1) { + buf_size = max_group + group = 0 + istar = 0 + } else { + Memi[DP_APID(apsel)] = Memi[DP_APID(apsel)+istar-1] + Memr[DP_APXCEN(apsel)] = Memr[DP_APXCEN(apsel)+istar-1] + Memr[DP_APYCEN(apsel)] = Memr[DP_APYCEN(apsel)+istar-1] + Memr[DP_APMAG(apsel)] = Memr[DP_APMAG(apsel)+istar-1] + Memr[DP_APMSKY(apsel)] = Memr[DP_APMSKY(apsel)+istar-1] + istar = 1 + } + + # Loop over the stars in a single group. + repeat { + + # Set the current group. + curr_group = group + + # Read in the photometry for a single star. + + # In this case we have a text database file. + if (key != NULL) { + + if (dp_nstsel (key, tp, fields, indices, Memi[DP_APID(apsel)+ + istar], group, Memr[DP_APXCEN(apsel)+istar], + Memr[DP_APYCEN(apsel)+istar], Memr[DP_APMSKY(apsel)+ + istar], Memr[DP_APMAG(apsel)+ istar]) == EOF) { + in_record = 0 + break + } + + # In this case we have a table. + } else { + + if (in_record > max_row) { + in_record = 0 + break + } else { + call tbrgti (tp, colpoint[1], Memi[DP_APID(apsel)+istar], + nullflag, 1, in_record) + call tbrgtr (tp, colpoint[2], Memr[DP_APXCEN(apsel)+istar], + nullflag, 1, in_record) + call tbrgtr (tp, colpoint[3], Memr[DP_APYCEN(apsel)+istar], + nullflag, 1, in_record) + call tbrgtr (tp, colpoint[4], Memr[DP_APMAG(apsel)+istar], + nullflag, 1, in_record) + call tbrgtr (tp, colpoint[5], Memr[DP_APMSKY(apsel)+istar], + nullflag, 1, in_record) + call tbrgti (tp, colpoint[6], group, nullflag, 1, in_record) + } + } + + # Increment the record and star counters. + in_record = in_record + 1 + istar = istar + 1 + + # Allocate more memory as needed. + if (istar == buf_size) { + buf_size = buf_size + max_group + call dp_rmemapsel (dao, indices, NAPPAR, buf_size + 1) + } + + } until ((group != curr_group) && (curr_group != 0)) + + # Return the number of stars in the group. + if (in_record == 0) { + if (curr_group == 0) + curr_group = group + return (istar) + } else + return (istar - 1) +end + + +# DP_NSTSEL -- Read in the required photometry records from a text file. + +int procedure dp_nstsel (key, fd, fields, indices, id, group, x, y, sky, mag) + +pointer key # pointer to key structure +int fd # text file descriptor +char fields[ARB] # fields to be output +int indices[ARB] # indices array +int id # star id number +int group # group number +real x # x center +real y # y center +real sky # sky value +real mag # magnitude + +int nchars, nunique, uunique, funique, ncontinue, recptr +int first_rec, nselect, record +pointer line +int getline(), strncmp(), pt_choose() +data first_rec /YES/ + +begin + # Initialize the file read. + if (first_rec == YES) { + nunique = 0 + uunique = 0 + funique = 0 + nselect = 0 + record = 0 + call malloc (line, SZ_LINE, TY_CHAR) + } + + # Initialize the record read. + ncontinue = 0 + recptr = 1 + + # Loop over the text file records. + repeat { + + # Read in a line of the text file. + nchars = getline (fd, Memc[line]) + if (nchars == EOF) + break + + # Determine the type of record. + if (Memc[line] == KY_CHAR_POUND) { + + if (strncmp (Memc[line], KY_CHAR_KEYWORD, KY_LEN_STR) == 0) { + call pt_kyadd (key, Memc[line], nchars) + } else if (strncmp (Memc[line], KY_CHAR_NAME, + KY_LEN_STR) == 0) { + nunique = nunique + 1 + call pt_kname (key, Memc[line], nchars, nunique) + } else if (strncmp (Memc[line], KY_CHAR_UNITS, + KY_LEN_STR) == 0) { + uunique = uunique + 1 + call pt_knunits (key, Memc[line], nchars, uunique) + } else if (strncmp (Memc[line], KY_CHAR_FORMAT, + KY_LEN_STR) == 0) { + funique = funique + 1 + call pt_knformats (key, Memc[line], nchars, funique) + } + + } else if (Memc[line] == KY_CHAR_NEWLINE) { + # skip blank lines + + } else { + + # Construct the text file record. + call pt_mkrec (key, Memc[line], nchars, first_rec, recptr, + ncontinue) + + # Construct output record when there is no continuation char. + if (Memc[line+nchars-2] != KY_CHAR_CONT) { + + # Select the appropriate records. + if (nselect <= 0) { + nselect = pt_choose (key, fields) + if (nselect < NAPGROUP) { + call eprintf ( + "The group file does not have the correct format\n") + break + } + } + + # Construct the output record by moving the selected fields + # into the data structures. + + call dp_gnst (key, indices, id, group, x, y, sky, mag) + record = record + 1 + first_rec = NO + + # Record is complete so exit the loop. + break + } + } + + } + + # Return EOF or the record number. + if (nchars == EOF || (nselect < NAPGROUP)) { + first_rec = YES + nunique = 0 + uunique = 0 + funique = 0 + nselect = 0 + call mfree (line, TY_CHAR) + return (EOF) + } else + return (record) +end + + +# DP_GNST -- Decode the standard GROUP text file fields into the appropriate +# arrays. + +procedure dp_gnst (key, fields, id, group, x, y, sky, mag) + +pointer key # pointer to keys strucuture +int fields[ARB] # fields array +int id # star id +int group # group id +real x # x position +real y # y position +real sky # sky value +real mag # magnitude + +int i, index, elem, maxch, kip, ip +int ctoi(), ctor() +char buffer[SZ_LINE] + +begin + do i = 1, KY_NSELECT(key) { + + # Find the key. + index = Memi[KY_SELECT(key)+i-1] + elem = Memi[KY_ELEM_SELECT(key)+i-1] + maxch = Memi[KY_LEN_SELECT(key)+i-1] + kip = Memi[KY_PTRS(key)+index-1] + (elem - 1) * maxch + call amovc (Memc[kip], buffer, maxch) + buffer[maxch+1] = EOS + + # Decode the output value. + ip = 1 + switch (fields[i]) { + case DP_PAPID: + if (ctoi (buffer, ip, id) <= 0) + call error (0, "ERROR: Error reading ID field.") + case DP_PAPGROUP: + if (ctoi (buffer, ip, group) <= 0) + call error (0, "ERROR: Error reading GROUP field.") + case DP_PAPXCEN: + if (ctor (buffer, ip, x) <= 0) + call error (0, "ERROR: Error reading XCENTER field.") + case DP_PAPYCEN: + if (ctor (buffer, ip, y) <= 0) + call error (0, "ERROR: Error reading YCENTER field.") + case DP_PAPSKY: + if (ctor (buffer, ip, sky) <= 0) + call error (0, "ERROR: Error reading MSKY field.") + case DP_PAPMAG1: + if (ctor (buffer, ip, mag) <= 0) + call error (0, "ERROR: Error reading MAG field.") + default: + call printf ("Error reading the photometry file.\n") + } + + } +end diff --git a/noao/digiphot/daophot/nstar/dpmemnstar.x b/noao/digiphot/daophot/nstar/dpmemnstar.x new file mode 100644 index 00000000..04bd25f4 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpmemnstar.x @@ -0,0 +1,158 @@ +include "../lib/daophotdef.h" +include "../lib/nstardef.h" + +# DP_NSTARSETUP -- Procedure to set up the NSTAR parameters. + +procedure dp_nstarsetup (dp) + +pointer dp # pointer to daophot structure + +pointer nstar + +begin + # Allocate Memory + call malloc (DP_NSTAR(dp), LEN_NSTARSTRUCT, TY_STRUCT) + nstar = DP_NSTAR(dp) + + DP_NNPIX(nstar) = NULL + DP_NNUMER(nstar) = NULL + DP_NDENOM(nstar) = NULL + DP_NRPIXSQ(nstar) = NULL + DP_NSKIP(nstar) = NULL + DP_NXCLAMP(nstar) = NULL + DP_NXOLD(nstar) = NULL + DP_NX(nstar) = NULL + DP_NV(nstar) = NULL + DP_NSUMWT(nstar) = NULL + DP_NC(nstar) = NULL + DP_NIER(nstar) = NULL +end + + +# DP_MEMNSTAR -- Procedure to allocate sufficient memory for NSTAR. + +procedure dp_memnstar (dao, max_star) + +pointer dao # pointer to daophot structure +int max_star # maximum number of stars + +pointer nstar + +begin + nstar = DP_NSTAR(dao) + + if (DP_NNPIX (nstar) != NULL) + call mfree (DP_NNPIX (nstar), TY_INT) + call malloc (DP_NNPIX (nstar), max_star + 1, TY_INT) + + if (DP_NNUMER (nstar) != NULL) + call mfree (DP_NNUMER (nstar), TY_REAL) + call malloc (DP_NNUMER (nstar), max_star + 1, TY_REAL) + + if (DP_NDENOM (nstar) != NULL) + call mfree (DP_NDENOM (nstar), TY_REAL) + call malloc (DP_NDENOM (nstar), max_star + 1, TY_REAL) + + if (DP_NRPIXSQ (nstar) != NULL) + call mfree (DP_NRPIXSQ (nstar), TY_REAL) + call malloc (DP_NRPIXSQ (nstar), max_star + 1, TY_REAL) + + if (DP_NSKIP (nstar) != NULL) + call mfree (DP_NSKIP (nstar), TY_INT) + call malloc (DP_NSKIP (nstar), max_star + 1, TY_INT) + + if (DP_NIER (nstar) != NULL) + call mfree (DP_NIER (nstar), TY_INT) + call malloc (DP_NIER(nstar), max_star + 1, TY_INT) + + if (DP_NSUMWT (nstar) != NULL) + call mfree (DP_NSUMWT (nstar), TY_REAL) + call malloc (DP_NSUMWT (nstar), max_star + 1, TY_REAL) + + if (DP_RECENTER(dao) == YES) { + + if (DP_NXCLAMP (nstar) != NULL) + call mfree (DP_NXCLAMP (nstar), TY_REAL) + call malloc (DP_NXCLAMP (nstar), 3 * max_star + 1, TY_REAL) + + if (DP_NXOLD (nstar) != NULL) + call mfree (DP_NXOLD (nstar), TY_REAL) + call malloc (DP_NXOLD (nstar), 3 * max_star + 1, TY_REAL) + + if (DP_NX (nstar) != NULL) + call mfree (DP_NX (nstar), TY_REAL) + call malloc (DP_NX (nstar), 3 * max_star + 1, TY_REAL) + + if (DP_NC (nstar) != NULL) + call mfree (DP_NC (nstar), TY_REAL) + call malloc (DP_NC (nstar), (3 * max_star + 1) * (3 * max_star + 1), + TY_REAL) + + if (DP_NV (nstar) != NULL) + call mfree (DP_NV (nstar), TY_REAL) + call malloc (DP_NV (nstar), 3 * max_star + 1, TY_REAL) + + } else { + + if (DP_NXCLAMP (nstar) != NULL) + call mfree (DP_NXCLAMP (nstar), TY_REAL) + call malloc (DP_NXCLAMP (nstar), max_star + 1, TY_REAL) + + if (DP_NXOLD (nstar) != NULL) + call mfree (DP_NXOLD (nstar), TY_REAL) + call malloc (DP_NXOLD (nstar), max_star + 1, TY_REAL) + + if (DP_NX (nstar) != NULL) + call mfree (DP_NX (nstar), TY_REAL) + call malloc (DP_NX (nstar), max_star + 1, TY_REAL) + + if (DP_NC (nstar) != NULL) + call mfree (DP_NC (nstar), TY_REAL) + call malloc (DP_NC (nstar), (max_star + 1) * (max_star + 1), + TY_REAL) + + if (DP_NV (nstar) != NULL) + call mfree (DP_NV (nstar), TY_REAL) + call malloc (DP_NV (nstar), max_star + 1, TY_REAL) + } +end + + +# DP_NSCLOSE -- Procedure to close up the NSTAR parameters. + +procedure dp_nsclose (dp) + +pointer dp # pointer to daophot structure + +pointer nstar + +begin + nstar = DP_NSTAR(dp) + + if (DP_NNPIX (nstar) != NULL) + call mfree (DP_NNPIX (nstar), TY_INT) + if (DP_NNUMER (nstar) != NULL) + call mfree (DP_NNUMER (nstar), TY_REAL) + if (DP_NDENOM (nstar) != NULL) + call mfree (DP_NDENOM (nstar), TY_REAL) + if (DP_NRPIXSQ (nstar) != NULL) + call mfree (DP_NRPIXSQ (nstar), TY_REAL) + if (DP_NSKIP (nstar) != NULL) + call mfree (DP_NSKIP (nstar), TY_INT) + if (DP_NXCLAMP (nstar) != NULL) + call mfree (DP_NXCLAMP (nstar), TY_REAL) + if (DP_NXOLD (nstar) != NULL) + call mfree (DP_NXOLD (nstar), TY_REAL) + if (DP_NX (nstar) != NULL) + call mfree (DP_NX (nstar), TY_REAL) + if (DP_NV (nstar) != NULL) + call mfree (DP_NV (nstar), TY_REAL) + if (DP_NSUMWT (nstar) != NULL) + call mfree (DP_NSUMWT (nstar), TY_REAL) + if (DP_NC (nstar) != NULL) + call mfree (DP_NC (nstar), TY_REAL) + if (DP_NIER (nstar) != NULL) + call mfree (DP_NIER (nstar), TY_INT) + + call mfree (nstar, TY_STRUCT) +end diff --git a/noao/digiphot/daophot/nstar/dpnconfirm.x b/noao/digiphot/daophot/nstar/dpnconfirm.x new file mode 100644 index 00000000..ecbd3056 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpnconfirm.x @@ -0,0 +1,34 @@ +include "../lib/daophotdef.h" + +# DP_NCONFIRM -- Procedure to confirm the critical nstar parameters. + +procedure dp_nconfirm (dao) + +pointer dao # pointer to the group structure + +int dp_stati() + +begin + call printf ("\n") + + # Confirm recentering and sky fitting. + call dp_vrecenter (dao) + call dp_vfitsky (dao) + if (dp_stati (dao, FITSKY) == NO) + call dp_vgroupsky (dao) + + # Confirm the psf radius. + call dp_vpsfrad (dao) + + # Confirm the fitting radius. + call dp_vfitrad (dao) + + # Confirm the maximum group size. + call dp_vmaxgroup (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/nstar/dpnstar.x b/noao/digiphot/daophot/nstar/dpnstar.x new file mode 100644 index 00000000..66ba10a4 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpnstar.x @@ -0,0 +1,355 @@ +include <imhdr.h> +include <tbset.h> +include <mach.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/nstardef.h" + + +# DP_NPHOT -- Read in the groups and fit all the stars in each group +# simultaneously. + +procedure dp_nphot (dao, im, grp, nst, rej, ap_text) + +pointer dao # pointer to the daophot structure +pointer im # input image descriptor +int grp # input group file descriptor +int nst # output photometry file descriptor +int rej # rejections file descriptor +bool ap_text # photometry text file + +bool converge +int out_record, rout_record, in_record, nrow_in_table, old_size, niter +int cdimen, stat +pointer sp, indices, fields, icolpoint, ocolpoint, key, nstar, apsel +real radius, mean_sky + +int dp_ggroup(), tbpsta(), dp_nstarfit(), dp_nxycheck() +real dp_nmsky() + +begin + # Store the original value of the fitting radius. + radius = DP_FITRAD(dao) + DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao)) + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) + + # Allocate working memory. + call smark (sp) + call salloc (icolpoint, NAPGROUP, TY_INT) + call salloc (indices, NAPPAR, TY_INT) + call salloc (fields, SZ_LINE, TY_CHAR) + call salloc (ocolpoint, NST_NOUTCOL, TY_INT) + + # Get some daophot pointers. + apsel = DP_APSEL(dao) + nstar = DP_NSTAR (dao) + + # Allocate space for and define the output table. + if (DP_TEXT(dao) == YES) { + call dp_xpnewnstar (dao, nst) + if (rej != NULL) + call dp_xpnewnstar (dao, rej) + } else { + call dp_tpnewnstar (dao, nst, Memi[ocolpoint]) + if (rej != NULL) + call dp_tpnewnstar (dao, rej, Memi[ocolpoint]) + } + + # Set up the input file. + call dp_gnindices (Memi[indices]) + if (ap_text) { + call pt_kyinit (key) + call dp_gnstpsf (Memi[indices], Memc[fields], NAPGROUP) + nrow_in_table = 0 + } else { + key = NULL + call dp_tnstinit (grp, Memi[icolpoint]) + nrow_in_table = tbpsta (grp, TBL_NROWS) + } + + # Allocate some memory for fitting. + call dp_memapsel (dao, Memi[indices], NAPPAR, DP_MAXGROUP(dao) + 1) + call dp_memnstar (dao, DP_MAXGROUP(dao) + 1) + + # Initialize the input/output files for reading/writing. + in_record = 1 + out_record = 0 + rout_record = 0 + + repeat { + + # Read in the next group of stars. + old_size = dp_ggroup (dao, grp, key, Memc[fields], Memi[indices], + Memi[icolpoint], nrow_in_table, DP_MAXGROUP(dao), in_record, + DP_NGNUM(nstar)) + if (old_size <= 0) + break + DP_NNUM(nstar) = old_size + + # Convert coordinates if necessary. + call dp_win (dao, im, Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], DP_NNUM(nstar)) + + # Print out the group number and number of stars. + if (DP_VERBOSE (dao) == YES) { + call printf ("Group: %4d contains %2d stars\n") + call pargi (DP_NGNUM (nstar)) + call pargi (DP_NNUM (nstar)) + } + + # If the group center is undefined or off image skip to next + # group. + if (dp_nxycheck (im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + DP_NNUM(nstar), DP_FITRAD(dao), stat) <= 0) { + + if (DP_VERBOSE(dao) == YES) { + if (stat < 0) + call printf ( + "\tGroup %d center is undefined: skipping\n") + else if (stat == 0) + call printf ("\tGroup %d is off the image: skipping\n") + call pargi (DP_NGNUM(nstar)) + } + next + + # If too many stars skip to the next group. + } else if (DP_NNUM (nstar) > DP_MAXGROUP(dao)) { + + if (DP_VERBOSE(dao) == YES) { + call printf ("\tGroup %d larger than %d stars: skipping\n") + call pargi (DP_NGNUM(nstar)) + call pargi (DP_MAXGROUP(dao)) + } + niter = 0 + call realloc (DP_NIER(nstar), DP_NNUM(nstar), TY_INT) + call dp_nstindef (dao, DP_NNUM(nstar), NSTERR_BIGGROUP) + + # If the mean sky value of group is undefined skip to next group. + } else if (IS_INDEFR (dp_nmsky (Memr[DP_APMSKY(apsel)], + DP_NNUM(nstar), mean_sky))) { + + if (DP_VERBOSE(dao) == YES) { + call printf ( + "\tGroup %d sky value is undefined: skipping\n") + call pargi (DP_NGNUM(nstar)) + } + niter = 0 + call dp_nstindef (dao, DP_NNUM(nstar), NSTERR_INDEFSKY) + + # Do the fit. + } else { + + # Estimate values of INDEF magnitudes. + call dp_nmaginit (dao, Memr[DP_APMAG(apsel)], + Memr[DP_APERR(apsel)], DP_NNUM(nstar)) + + # Set up for the fit. + if (DP_RECENTER(dao) == YES) + cdimen = 3 * DP_NNUM(nstar) + 1 + else + cdimen = DP_NNUM(nstar) + 1 + + # Iterate. + for (niter = 1; niter <= DP_MAXITER(dao); niter = niter + 1) { + DP_NNUM(nstar) = dp_nstarfit (dao, im, DP_NNUM(nstar), + mean_sky, cdimen, niter, converge) + if (DP_NNUM(nstar) <= 0) + break + if (converge) + break + } + + # Solution did not converge. + niter = min (niter, DP_MAXITER(dao)) + } + + # Now write out the results. + if (DP_TEXT(dao) == YES) + call dp_xntwrite (dao, im, nst, rej, niter, old_size) + else + call dp_tntwrite (dao, im, nst, rej, niter, old_size, + out_record, rout_record, Memi[ocolpoint]) + } + + if (ap_text) + call pt_kyfree (key) + + call sfree (sp) + + # Restore the original value of the fitting radius. + DP_FITRAD(dao) = radius + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) +end + + +# DP_NSTINDEF -- Set magnitudes and fitting parameters of unfit stars to indef. + +procedure dp_nstindef (dao, max_nstars, pier) + +pointer dao # pointer to the daophot structure +int max_nstars # number of stars +int pier # the photometry error code + +int i +pointer apsel, nstar + +begin + apsel = DP_APSEL(dao) + nstar = DP_NSTAR(dao) + do i = 1, max_nstars { + Memr[DP_APMAG(apsel)+i-1] = INDEFR + Memr[DP_APERR(apsel)+i-1] = INDEFR + Memr[DP_APCHI(apsel)+i-1] = INDEFR + Memr[DP_APSHARP(apsel)+i-1] = INDEFR + Memi[DP_NIER(nstar)+i-1] = pier + } +end + + +# DP_GNINDICES -- Get the memory allocation fields. + +procedure dp_gnindices (indices) + +int indices[ARB] # index array + +begin + indices[1] = DP_PAPID + indices[2] = DP_PAPXCEN + indices[3] = DP_PAPYCEN + indices[4] = DP_PAPMAG1 + indices[5] = DP_PAPSKY + indices[6] = DP_PAPGROUP + indices[7] = DP_PAPMERR1 + indices[8] = DP_PAPNITER + indices[9] = DP_PAPSHARP + indices[10] = DP_PAPCHI +end + + +define GRP_CTRINDEF -1 +define GRP_CTROFFIMAGE 0 +define GRP_CTROK 1 + +# DP_NXYCHECK -- Check that the center of the group is defined. -1 is +# returned if the center of the group is undefined, 0, is returned if the +# group is entirely off the image, 1 is returned if the group is at least +# partially on the input image. + +int procedure dp_nxycheck (im, x, y, group_size, radius, stat) + +pointer im # pointer to the input image +real x[ARB] # array of x values +real y[ARB] # array of y values +int group_size # the size of the group +real radius # the fitting radius +int stat # the return status + +int i, nxy +real xmin, xmax, ymin, ymax, xx, yy + +begin + # Initialize. + stat = GRP_CTRINDEF + xmin = MAX_REAL + xmax = -MAX_REAL + ymin = MAX_REAL + ymax = -MAX_REAL + + # Compute the minimum and maximum x and y values. + nxy = 0 + do i = 1, group_size { + xx = x[i] + yy = y[i] + if (IS_INDEFR(xx) || IS_INDEFR(yy)) + next + nxy = nxy + 1 + if (xx < xmin) + xmin = xx + if (xx > xmax) + xmax = xx + if (yy < ymin) + ymin = yy + if (yy > ymax) + ymax = yy + } + if (nxy <= 0) + return (stat) + + # Test the min and max values. + stat = GRP_CTROFFIMAGE + if ((int (xmin - radius) + 1) > IM_LEN(im,1)) + return (stat) + if (int (xmax + radius) < 1) + return (stat) + if ((int (ymin - radius) + 1) > IM_LEN(im,2)) + return (stat) + if (int (ymax + radius) < 1) + return (stat) + + # The group is on the image. + stat = GRP_CTROK + return (stat) +end + + +# DP_NMSKY -- Compute the mean sky value for the group of stars. + +real procedure dp_nmsky (sky, group_size, msky) + +real sky[ARB] # the array of sky values +int group_size # the size of the group of stars +real msky # the mean sky value + +int i, nsky +real sky_sum + +begin + sky_sum = 0.0 + nsky = 0 + + do i = 1, group_size { + if (IS_INDEFR(sky[i])) + next + sky_sum = sky_sum + sky[i] + nsky = nsky + 1 + } + + if (nsky <= 0) + msky = INDEFR + else + msky = sky_sum / nsky + + return (msky) +end + + +define MIN_REL_BRIGHT 1.0E-04 # minimum relative brightness +define INIT_REL_BRIGHT 0.01 # initial relative brightness + +# DP_NMAGINIT -- Initialize the magnitude and magnitude error arrays before +# fitting the group. + +procedure dp_nmaginit (dao, mag, magerr, group_size) + +pointer dao # pointer to the daophot strucuture +real mag[ARB] # the magnitude array +real magerr[ARB] # the magnitude error array +int group_size # size of the group + +int i +pointer psffit + +begin + psffit = DP_PSFFIT (dao) + do i = 1, group_size { + if (IS_INDEFR(mag[i])) + mag[i] = INIT_REL_BRIGHT + else { + mag[i] = DAO_RELBRIGHT (psffit, mag[i]) + if (mag[i] <= MIN_REL_BRIGHT) + mag[i] = INIT_REL_BRIGHT + } + magerr[i] = 0.0 + } +end diff --git a/noao/digiphot/daophot/nstar/dpnstarfit.x b/noao/digiphot/daophot/nstar/dpnstarfit.x new file mode 100644 index 00000000..8b864c86 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpnstarfit.x @@ -0,0 +1,1383 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/nstardef.h" + + +# DP_NSTARFIT -- Fit the stellar group. + +int procedure dp_nstarfit (dao, im, nin_group, mean_sky, cdimen, iter, + converge) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +int nin_group # original group size +real mean_sky # the mean sky for the group +int cdimen # dimensions of the coefficient matrix +int iter # the current iteration +bool converge # did the fit converge ? + +bool clip, refit +int i, j, k, xpix, ypix, group_size, nterm, ncols, mindex, ifaint, tifaint +int ntmin, ixmin, ixmax, iymin, iymax, ifxmin, ifxmax, ifymin, ifymax, flag +pointer apsel, psffit, nstar, subim, ypixel, pixel +real fitradsq, cutoff, psfradsq, sepcrit, sepmin, perr, peakerr, sky_value +real read_noise, mingdata, maxgdata, chigrp, wcrit, xmin, xmax, ymin, ymax +real datum, sumres, grpwt, xtemp, ytemp, weight, ds, pred_pixval, sigmasq +real relerr, dswt, faint, tfaint, noise + +bool dp_nstmerge(), dp_ntomit(), dp_ntmin(), dp_ncheckc() +pointer imgs2r() +real dp_ntskyval(), dp_ntsubtract() + +begin + # Define the daophot pointers. + psffit = DP_PSFFIT (dao) + apsel = DP_APSEL(dao) + nstar = DP_NSTAR (dao) + + # Set up some daophot constants. At some point these will be computed + # when the NSTAR task is started up instead of at the beginning of + # each group fit. For the moment it is convenient and not too + # costly to compute them here. Also initialize the fit. + + if (iter == 1) { + + # Compute the fitting and psf radii. + fitradsq = DP_FITRAD (dao) ** 2 + psfradsq = DP_PSFRAD(dao) ** 2 + cutoff = CUT_FACTOR * fitradsq + + # Compute the merging parameters. + if (IS_INDEFR(DP_MERGERAD(dao))) { + sepcrit = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 + + Memr[DP_PSFPARS(psffit)+1] ** 2) + sepmin = min (1.0, FRACTION_MINSEP * sepcrit) + } else { + sepcrit = DP_MERGERAD(dao) ** 2 + sepmin = min (1.0, FRACTION_MINSEP * sepcrit) + } + + # Compute the noise parameters. + read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2 + perr = 0.01 * DP_FLATERR(dao) + peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] * + Memr[DP_PSFPARS(psffit)+1]) + + # Compute the minimum and maximum good pixel values. + 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) + + # Initialize the fit. + chigrp = 1.0 + clip = false + if (DP_RECENTER(dao) == YES) { + nterm = 3 * nin_group + ntmin = 3 + } else { + nterm = nin_group + ntmin = 1 + } + if (DP_FITSKY(dao) == YES) { + nterm = nterm + 1 + ntmin = ntmin + 1 + } + call aclrr (Memr[DP_NXOLD(nstar)], nterm) + call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) + } + + # Start a new iteration. + group_size = nin_group + if (DP_RECENTER(dao) == YES) + nterm = 3 * group_size + else + nterm = group_size + if (DP_FITSKY(dao) == YES) + nterm = nterm + 1 + + # Begin fitting the current group of stars. + repeat { + + # Initialize the convergence criteria. + converge = false + + # Check that there is at least 1 star in the group. + if (group_size < 1) + return (group_size) + + # Set up the critical error for star rejection. + if (iter >= NITER_MAX) + wcrit = WCRIT_MAX + else if (iter >= NITER_MED) + wcrit = WCRIT_MED + else if (iter >= NITER_MIN) + wcrit = WCRIT_MIN + else + wcrit = MAX_REAL + + # Set the sky fitting derivative. + if (DP_FITSKY(dao) == YES) + Memr[DP_NX(nstar)+nterm-1] = -1.0 + + # Initialize arrays. + call aclrr (Memr[DP_APCHI(apsel)], group_size) + call aclrr (Memr[DP_NSUMWT(nstar)], group_size) + call aclrr (Memr[DP_NNUMER(nstar)], group_size) + call aclrr (Memr[DP_NDENOM(nstar)], group_size) + call amovki (NSTERR_OK, Memi[DP_NIER(nstar)], group_size) + + # Compute the minimum and maximum x and y values. + call alimr (Memr[DP_APXCEN(apsel)], group_size, xmin, xmax) + call alimr (Memr[DP_APYCEN(apsel)], group_size, ymin, ymax) + + # Check to see whether any two stars are within the critical + # difference from each other. + + if ((group_size > 1) && dp_nstmerge (Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APERR(apsel)], group_size, sepcrit, sepmin, wcrit, + i, j, k)) { + + # Compute the new centroid and brightness. + call dp_nstcen (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], i, j, k) + + # Print out verbose comments. + if (DP_VERBOSE (dao) == YES) { + call printf ("\tStar %-5d has merged with star %-5d\n") + call pargi (Memi[DP_APID(apsel)+k-1]) + call pargi (Memi[DP_APID(apsel)+i-1]) + } + + # Recompute the mean sky. + if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES)) + mean_sky = (mean_sky * group_size - + Memr[DP_APMSKY(apsel)+k-1]) / (group_size - 1) + + # Now remove the k-th star from the group. + call dp_remove (k, group_size, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], + Memi[DP_NIER(nstar)], NSTERR_MERGE) + + # After deleting a star, resize the matrix, release all of + # the clamps and back up the iteration counter. + + if (DP_RECENTER(dao) == YES) + nterm = 3 * group_size + else + nterm = group_size + if (DP_FITSKY(dao) == YES) + nterm = nterm + 1 + clip = false + call aclrr (Memr[DP_NXOLD(nstar)], nterm) + call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) + iter = max (1, iter - 1) + next + } + + # Now we can proceed with the iteration. If this is the first + # iteration read in the subraster. Determine the size of the + # subraster we need to extract from the image for this group. + + if (iter == 1) { + ixmin = max (1, int (xmin - DP_PSFRAD(dao) - + DP_FITRAD(dao)) + 1) + iymin = max (1, int (ymin - DP_PSFRAD(dao) - + DP_FITRAD(dao)) + 1) + ixmax = min (IM_LEN(im,1), int (xmax + DP_PSFRAD(dao) + + DP_FITRAD(dao))) + iymax = min (IM_LEN(im,2), int (ymax + DP_PSFRAD(dao) + + DP_FITRAD(dao))) + subim = imgs2r (im, ixmin, ixmax, iymin, iymax) + ncols = ixmax - ixmin + 1 + } + + # Compute the area on the subraster that is off interest to + # the current iteration. + + ifxmin = max (ixmin, int (xmin - DP_FITRAD(dao)) + 1) + ifymin = max (iymin, int (ymin - DP_FITRAD(dao)) + 1) + ifxmax = min (ixmax, int (xmax + DP_FITRAD(dao))) + ifymax = min (iymax, int (ymax + DP_FITRAD(dao))) + + # Zero the normal matrix and the vector of residuals. + + call aclrr (Memr[DP_NV(nstar)], nterm) + call aclrr (Memr[DP_NC(nstar)], cdimen * cdimen) + call aclri (Memi[DP_NNPIX(nstar)], group_size) + + sumres = 0.0 + grpwt = 0.0 + + # Loop over the pixels in the subraster. + + ypixel = subim + (ifymin - iymin) * ncols + ifxmin - ixmin - 1 + do ypix = ifymin, ifymax { + + ytemp = ypix + pixel = ypixel + + do xpix = ifxmin, ifxmax { + + xtemp = xpix + pixel = pixel + 1 + datum = Memr[pixel] + + # Reject pixel if not in good data range. + if ((datum < mingdata) || (datum > maxgdata)) + next + + # If this pixel is within one fitting radius of at + # least one star then include it in the solution. + # While figuring this out, compute the squared distance + # of this pixel from the centroid of each star in the + # group, and sum its contribution to the number + # of pixels within one fitting radius of each star. + + if (dp_ntomit (Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)], + Memi[DP_NSKIP(nstar)], group_size, xtemp, ytemp, + cutoff)) + next + + if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == NO)) { + sky_value = dp_ntskyval (Memr[DP_APMSKY(apsel)], + Memi[DP_NSKIP(nstar)], group_size) + if (IS_INDEFR(sky_value)) + sky_value = mean_sky + } else + sky_value = mean_sky + + # Subtract the mean sky from the pixel. + #ds = datum - mean_sky + ds = datum - sky_value + + # Now loop over the stars and subtract from this pixel + # the light contribution from each star within one psf + # radius. + # + # The condition equation for pixel[xpix,ypix] is + # + # residual = data[xpix,ypix] - sky - sum (scale * + # psf[xpix-xc,ypix-yc]) + # + # The scale, xc, and yc are iterated until + # + # sum (weight * residual ** 2) + # + # is minimized. Weight will be a function of 1) the + # distance of the pixel from the center of the nearest + # star, 2) the model-predicted brightness of the pixel + # (taking into consideration the readout noise, the + # photons/ADU, and the interpolating error of the PSF, + # and, 3) the size of the residual itself. One is + # necessary to prevent the non-linear least squares + # solution from oscillating: oftimes it will come + # to pass that if you include a pixel in the solution + # then the predicted shift of the centroid will cause + # that pixel to be excluded in the next iteration, and the + # new predicted shift of the centroid will cause that + # pixel to be included again. This could go on ad + # infinitum. The cure is to have the weight of the pixel + # go asymptotically to zero as its distance from the + # stellar centroid approaches the fitting radius. In + # the case like that just described, the solution can + # then find a real minimum of the sum of the weighted + # squared residuals with that pixel at some low-weight + # position just inside the fitting radius. Two is + # just sensible weighting. Three is a crude attempt + # at making the solution more robust against bad pixels. + + weight = dp_ntsubtract (dao, im, Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)], + Memr[DP_APMAG(apsel)], Memi[DP_NSKIP(nstar)], + Memr[DP_NX(nstar)], group_size, xtemp, ytemp, ds, + psfradsq, fitradsq, DP_RECENTER(dao)) + + # At this point the vector X contains the first + # derivative of the condition equation, for the pixel + # under consideration, with respect to each of the + # fitting parameters for all of these stars. We now need + # to add these values into the normal matrix and the vector + # of residuals. + + # The expected random error in the pixel is the quadratic + # sum of the Poisson statistics, plus the readout noise, + # plus an estimated error of 0.75% of the total brightness + # of the total brightness for the difficulty of flat- + # fielding and bias subtraction, plus an estimated error + # of the same fraction of the fourth derivative at the + # peak of the profile, to account for the difficulty + # of accurately interpolating within the point-spread + # function. The fourth derivative of the PSF is + # proportional to H / sigma ** 4 (sigma is the Gaussian + # width parameter for the stellar core); using the geometric + # mean of sigma(x) and sigma(y), this becomes H / [sigma(x) + # * sigma(y)] ** 2. The ratio of the fitting error to this + # quantity is estimated to be approximately 0.027 from a + # good-seeing CTIO frame. (This is probably a function of + # seeing, sampling etc.) + + # Get the residual from the PSF fit and the pixel + # intensity as predicted by the fit. Pred_pixval = raw data + # minus residual = model predicted value of the intensity at + # this point. + + pred_pixval = max (0.0, datum - ds) + if ((pred_pixval > maxgdata) && (iter >= 4)) + next + + #sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + + #(perr * pred_pixval) ** 2 + (peakerr * + #(pred_pixval - mean_sky)) ** 2 + sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + + (perr * pred_pixval) ** 2 + (peakerr * + (pred_pixval - sky_value)) ** 2 + if (sigmasq <= 0.0) + next + relerr = abs (ds) / sqrt (sigmasq) + + # Add this residual into the weighted sum of the + # absolute relative residuals. + + sumres = sumres + relerr * weight + grpwt = grpwt + weight + + # Add into the accumulating sums of the weighted + # absolute relative residuals and of the image sharpness + # parameter for each of the stars. Include in the + # sharpness index only those pixels within NCORE_SIGMA + # sigma of the centroid of the object. This saves time + # and floating underflows by excluding pixels + # which contribute less than about one part in a + # million to the index. + + call dp_acsharp (Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memi[DP_NSKIP(nstar)], + Memi[DP_NNPIX(nstar)], Memr[DP_NNUMER(nstar)], + Memr[DP_NDENOM(nstar)], Memr[DP_NSUMWT(nstar)], + Memr[DP_APCHI(apsel)], group_size, xtemp, ytemp, + Memr[DP_PSFPARS(psffit)], Memr[DP_PSFPARS(psffit)+1], + ds, sigmasq, relerr, weight) + + # If clipping is in effect, reduce the weight of a bad + # pixel. A pixel having a residual of 2.5 sigma gets + # reduced to half weight and one with a rersidual of 5 + # sigma gets weight of 1 / 257. + + weight = weight / sigmasq + if (clip) + weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) * + chigrp)) ** DP_CLIPEXP(dao)) + dswt = ds * weight + + # Work this pixel into the normal matrix. + call dp_mataccum (Memr[DP_NX(nstar)], Memi[DP_NSKIP(nstar)], + group_size, Memr[DP_NC(nstar)], Memr[DP_NV(nstar)], + cdimen, nterm, weight, dswt, DP_RECENTER(dao), + DP_FITSKY(dao)) + + } + + ypixel = ypixel + ncols + } + + # Make sure that every star in the group has at least MIN_FIT_PIXEL + # pixels within one fitting radius. + + refit = dp_ntmin (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], Memi[DP_NNPIX(nstar)], + Memi[DP_NIER(nstar)], group_size, nterm, DP_RECENTER(dao), + DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, DP_VERBOSE(dao)) + if (group_size < 1) + return (group_size) + if (refit) + next + + # Reflect the normal matrix across the diagonal. + call dp_mreflect (Memr[DP_NC(nstar)], cdimen, nterm) + + # Compute the robust estimate of the standard deviation of the + # residuals for the group as a whole, and for each star. This + # estimate is sqrt (PI/2) * weighted mean absolute relative + # residual. Do you like that "absolute relative" stuff? (PBS) + # NO! (LED) + # + # CHI = CHI_NORM * SUM (weight * resid) / (# of pixels) + # + # This gets corrected for bias by being multiplied by + # + # SQRT (# of pixels) / (# of pixels - 3) + + # Then the value is driven towards unity, depending on + # exactly how many pixels were involved: if CHIOLD is based + # on a total weight of 3, then it is extremely poorly determined + # and we just want to keep CHIOLD = 1. The larger GRPWT is, the + # better determined CHIOLD is, and the less we want to force it + # toward unity. So, just take the weighted average of CHIOLD and + # unity, with weights GRPWT - 3 and 1, respectively. + + if (grpwt > ntmin) { + chigrp = CHI_NORM * sumres * sqrt (1.0 / (grpwt * (grpwt - + ntmin))) + chigrp = ((grpwt - ntmin) * chigrp + ntmin) / grpwt + } + + # CHIOLD has been pulled toward its expected value of unity to + # keep the statistics of a small number of pixels from completely + # dominating the error analysis. Similarly, the photometric + # errors for the individual stars will be pulled toward unity + # now. Later on, if the number of stars in the group is + # greated than one, CHI will be nudged toward the group average. + # In order to work optimally, of course, this requires that + # the # of photons / ADU, the READNOISE and the other noise + # contributors are properly specified. + + call dp_ntchi (Memr[DP_NSUMWT(nstar)], Memr[DP_APCHI(apsel)], + group_size, ntmin, chigrp, grpwt) + + # Invert the matrix. + call invers (Memr[DP_NC(nstar)], cdimen, nterm, flag) + + # Check for a singular matrix. + refit = dp_ncheckc (Memr[DP_NC(nstar)], cdimen, nterm, + Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size, + DP_RECENTER(dao), DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, + DP_VERBOSE(dao)) + if (group_size < 1) + return (group_size) + if (refit) + next + + # Solve for position and scale factor increments. + call mvmul (Memr[DP_NC(nstar)], cdimen, nterm, Memr[DP_NV(nstar)], + Memr[DP_NX(nstar)]) + + if (iter <= 1) + refit = true + else + refit = false + + # Fit the sky. + if (DP_FITSKY(dao) == YES) { + noise = sqrt (abs (mean_sky / DP_PHOTADU(dao)) + read_noise) + mean_sky = mean_sky - max (-3.0 * noise, + min (Memr[DP_NX(nstar)+nterm-1], 3.0 * noise)) + if (abs (Memr[DP_NX(nstar)+nterm-1]) > (1.0e-4 * mean_sky)) + refit = true + } + + # In the beginning, the brightness of each star will be permitted + # to change by no more than 2 magnitudes per iteration, and the x,y + # coordinates of each centroid will be permitted to change by + # no more than 0.4 pixels per iteration. Any time that the + # parameter correction changes sign from one iteration to the + # next, the maximum permissible change will be reduced by a factor + # of two. These clamps are released any time a star disappears. + + call dp_ntclamp (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)], + Memr[DP_NSUMWT(nstar)], group_size, Memr[DP_NC(nstar)], cdimen, + Memr[DP_NX(nstar)], Memr[DP_NXOLD(nstar)], + Memr[DP_NXCLAMP(nstar)], DP_RECENTER(dao), clip, refit) + + # Check whether the estimated centroid of the any star has + # moved so far out of the limits of the picture that it has fewer + # than 4 or 5 pixels within one fitting radius. + + call dp_ntcentroid (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size, + ixmin, ixmax, iymin, iymax, fitradsq, DP_FITSKY(dao), + DP_GROUPSKY(dao), mean_sky, refit, DP_VERBOSE(dao)) + + if (group_size < 1) + return (group_size) + + # Update matrix dimensions. + if (DP_RECENTER(dao) == YES) + nterm = 3 * group_size + else + nterm = group_size + if (DP_FITSKY(dao) == YES) + nterm = nterm + 1 + + # Now check whether any of the stars is too faint (more than 12.5 + # magnitudes fainter than the PSF star). If several stars are too + # faint, delete the faintest one, and set the brightness of the + # other faint ones to 12.5 magnitudes below the PSF star. That way + # on the next iteration we will see whether these stars want to + # grow or to disappear. + + faint = 0.0 + ifaint = 0 + call dp_ntfmag (Memr[DP_APMAG(apsel)], group_size, tfaint, tifaint) + + # If at least one star is more than 12.5 magnitudes fainter + # than the PSF then ifaint is the relative index of the faintest + # of them, and faint is the relative brightness of the + # faintest of them. + + # No very faint star was detected. + if (tifaint <= 0) { + + # If the solution has not converged and if the number of + # iterations is less than MIN_ITER, perform another iteration + # with no questions asked. + + if ((refit) && (iter < MIN_ITER)) + return (group_size) + + # If the solution doesn't think it has converged, after the + # fourth iteration delete the least certain star if it is less + # less than a one-sigma detection; after the eighth iteration + # delete the least certain star if it is less than a 1.5 sigma + # detection; after the twelfth iteration OR if the solution + # thinks it has converged, delete the least certain star if it + # is less than a two-sigma detection. + + call dp_fsnoise (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)], + group_size, faint, ifaint) + + if ((refit) && (iter < DP_MAXITER (dao)) && (faint < wcrit)) + return (group_size) + } + + # Reject the appropriate star. + if ((tifaint > 0) || (faint >= MIN_FAINT)) { + + # Either (a) the solution thinks it has not converged + # and the faintest star is more uncertain than sqrt(wcrit) + # or (b) the solution thinks it has converged and the + # faintest star is more uncertain than two-sigma. + + if (DP_VERBOSE (dao) == YES) { + mindex = max (tifaint, ifaint) + call printf ( + "\tStar %-5d has been deleted because it is too faint\n") + call pargi (Memi[DP_APID(apsel)+mindex-1]) + } + + if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES) && + (group_size > 1)) + mean_sky = (mean_sky * group_size - + Memr[DP_APMSKY(apsel)+max(tifaint,ifaint)-1]) / + (group_size - 1) + + call dp_remove (max (tifaint, ifaint), group_size, + Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], + NSTERR_FAINT) + + if (group_size < 1) + return (group_size) + + if (DP_RECENTER(dao) == YES) + nterm = 3 * group_size + else + nterm = group_size + if (DP_FITSKY(dao) == YES) + nterm = nterm + 1 + call aclrr (Memr[DP_NXOLD(nstar)], nterm) + call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) + clip = false + iter = max (1, iter - 1) + next + } + + # Solution has either converged or gone to the maximum number + # of iterations. + + if ((iter < DP_MAXITER(dao)) && (! clip)) { + + # The first convergence milestone has been reached. Turn on the + # clipper, loosen the clamps and keep on going. + + if (DP_CLIPEXP(dao) > 0) + clip = true + converge = false + call aclrr (Memr[DP_NXOLD(nstar)], nterm) + call amaxkr (Memr[DP_NXCLAMP(nstar)], 0.25, + Memr[DP_NXCLAMP(nstar)], nterm) + return (group_size) + } + + converge = true + + } until (converge) + + return (group_size) +end + + +# DP_NSTMERGE -- Decide whether two stars in a group should merge. + +bool procedure dp_nstmerge (xcen, ycen, mag, magerr, group_size, sepcrit, + sepmin, wcrit, i, j, k) + +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real mag[ARB] # array of magnitudes +real magerr[ARB] # array of magnitude errors +int group_size # group size +real sepcrit # the critical separation squared +real sepmin # the minimum separation +real wcrit # critical error for rejection +int i, j, k # output indices + +real separation + +begin + do i = 1, group_size { + do j = 1, i - 1 { + + # Compute the separation. + separation = (xcen[j] - xcen[i]) ** 2 + + (ycen[j] - ycen[i]) ** 2 + if (separation > sepcrit) + next + + # Find the fainter of the two stars. + k = j + if (mag[i] < mag[j]) + k = i + if ((separation < sepmin) || ((magerr[k] / mag[k]) > wcrit)) + return (true) + } + } + + return (false) +end + + +# DP_NSTCEN -- Recompute the centroid and brightness of the i-th star. + +procedure dp_nstcen (xcen, ycen, mag, i, j, k) + +real xcen[ARB] # the x centers +real ycen[ARB] # the y centers +real mag[ARB] # the magnitudes +int i, j, k # array indices + +begin + # Now eliminate the fainter of the two stars. The k-th + # star is now the fainter of the two, the i-th the + # brighter. + + if (mag[i] < mag[j]) + i = j + + # Now replace the centroid of the i-th star with the + # weighted mean of the most recent estimates of the + # centroids of the i-th and the k-th stars, and the + # brightness of i-th with the sum of the brightnesses. + + xcen[i] = xcen[i] * mag[i] + xcen[k] * mag[k] + ycen[i] = ycen[i] * mag[i] + ycen[k] * mag[k] + mag[i] = mag[i] + mag[k] + xcen[i] = xcen[i] / mag[i] + ycen[i] = ycen[i] / mag[i] +end + + +# DP_NTOMIT -- Check whether a pixel is within one fitting radius of another +# star. + +bool procedure dp_ntomit (xcen, ycen, rpixsq, skip, group_size, fx, fy, cutoff) + +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real rpixsq[ARB] # array of distances squared +int skip[ARB] # array of skip values +int group_size # the group size +real fx, fy # pixel position in image +real cutoff # radius cuttoff for pixel inclusion + +bool omit +int i + +begin + omit = true + do i = 1, group_size { + skip[i] = YES + rpixsq[i] = (fx - xcen[i]) ** 2 + (fy - ycen[i]) ** 2 + if (rpixsq[i] > cutoff) + next + skip[i] = NO + omit = false + } + + return (omit) +end + + +# DP_NTSKYVAL -- Compute the average sky value to use for a particular +# pixel by averaging the sky values of all stars for which the +# pixel is within one fitting radius. + +real procedure dp_ntskyval (sky, skip, nstar) + +real sky[ARB] # array of sky values +int skip[ARB] # array of skip values +int nstar # the number of stars + +int i, npts +real sum + +begin + sum = 0.0 + npts = 0 + do i = 1, nstar { + if (skip[i] == YES) + next + sum = sum + sky[i] + npts = npts + 1 + } + if (npts <= 0) + return (INDEFR) + else + return (sum / npts) +end + + +# DP_NTSUBTRACT -- Procedure to subtract the contribution of a particular +# pixel from a particular star. + +real procedure dp_ntsubtract (dao, im, xcen, ycen, rpixsq, mag, skip, x, + group_size, fx, fy, ds, psfradsq, fitradsq, recenter) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real rpixsq[ARB] # array of distances squared +real mag[ARB] # array of magnitudes +int skip[ARB] # array of skip values +real x[ARB] # x accumulation vector array +int group_size # size of the group +real fx, fy # center of pixel in image +real ds # pixel value +real psfradsq # psf radius squared +real fitradsq # fit radius squared +int recenter # recenter the coordinates + +int i, i3, k +pointer psffit +real weight, dx, dy, deltax, deltay, val, dvdx, dvdy, rsq +real dp_usepsf() + +begin + psffit = DP_PSFFIT(dao) + + weight = 0.0 + do i = 1, group_size { + if (rpixsq[i] >= psfradsq) + next + dx = fx - xcen[i] + dy = fy - ycen[i] + call dp_wpsf (dao, im, xcen[i], ycen[i], deltax, deltay, 1) + deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0 + deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0 + val = 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, + dvdx, dvdy) + ds = ds - mag[i] * val + if (skip[i] == YES) + next + rsq = rpixsq[i] / fitradsq + weight = max (weight, 5.0 / (5.0 + rsq / (1.0 - rsq))) + if (recenter == YES) { + i3 = 3 * i + k = i3 - 2 + x[k] = -val + k = i3 - 1 + x[k] = -mag[i] * dvdx + x[i3] = -mag[i] * dvdy + } else + x[i] = -val + } + + return (weight) +end + + +# DP_ACSHARP -- Procedure to accumulate sums of the weighted absolute +# relative residuals and the image sharpness parameter for each of the +# stars. + +procedure dp_acsharp (xcen, ycen, skip, npix, numer, denom, sumwt, chi, + group_size, fx, fy, fwhmx, fwhmy, ds, sigmasq, relerr, weight) + +real xcen[ARB] # array of object x centers +real ycen[ARB] # array of object y centers +int skip[ARB] # array of skip values +int npix[ARB] # array of numbers of pixels +real numer[ARB] # numerator array +real denom[ARB] # denominator array +real sumwt[ARB] # array of summed weights +real chi[ARB] # array of chis +int group_size # group size paramter. +real fx, fy # position of data in image +real fwhmx, fwhmy # gaussian core widths in x and y +real ds # the data residual +real sigmasq # the sigma squared +real relerr # the relative error +real weight # the weight + +int i +real rhosq, dfdsig + +begin + do i = 1, group_size { + + if (skip[i] == YES) + next + + # Accumulate the number of pixels, chi and sum of the weights. + npix[i] = npix[i] + 1 + chi[i] = chi[i] + relerr * weight + sumwt[i] = sumwt[i] + weight + + # Include in the sharpness index only those pixels + # within NCORE_SIGMASQ of the centroid of the + # object. (This saves time and floating underflows + # by excluding pixels which contribute very little + # to the index. + + rhosq = ((xcen[i] - fx) / fwhmx) ** 2 + ((ycen[i] - fy) / + fwhmy) ** 2 + if (rhosq > NCORE_SIGMASQ) + next + rhosq = 0.6931472 * rhosq + dfdsig = exp (-rhosq) * (rhosq - 1.0) + numer[i] = numer[i] + dfdsig * ds / sigmasq + denom[i] = denom[i] + (dfdsig ** 2) / sigmasq + } +end + + +# DP_MATACCUM -- Procedure to accumulate the data into the matrices. + +procedure dp_mataccum (x, skip, group_size, c, v, cdimen, nterm, weight, dswt, + recenter, fitsky) + +real x[ARB] # x array +int skip[ARB] # skip vector +int group_size # size of the group +real c[cdimen,ARB] # coefficient matrix +real v[ARB] # vector array +int cdimen # dimensions of the coefficient matrix +int nterm # the number of terms +real weight # weight +real dswt # data weight +int recenter # recenter the coordinates +int fitsky # fit the sky value + +int i, i3, i3m2, k, j, l + +begin + if (fitsky == YES) { + c[nterm,nterm] = c[nterm,nterm] + weight + v[nterm] = v[nterm] - dswt + } + + do i = 1, group_size { + if (skip[i] == YES) + next + if (recenter == YES) { + i3 = i * 3 + i3m2 = i3 - 2 + do k = i3m2, i3 { + if (fitsky == YES) + c[nterm,k] = c[nterm,k] - x[k] * weight + v[k] = v[k] + x[k] * dswt + } + do j = 1, i { + if (skip[j] == YES) + next + do k = i3m2, i3 { + do l = 3 * j - 2, min (k, 3 * j) + c[k,l] = c[k,l] + x[k] * x[l] * weight + } + } + } else { + v[i] = v[i] + x[i] * dswt + if (fitsky == YES) + c[nterm,i] = c[nterm,i] - x[i] * weight + do j = 1, i { + if (skip[j] == YES) + next + c[i,j] = c[i,j] + x[i] * x[j] * weight + } + } + } +end + + +# DP_NTMIN -- Make sure that every star in the group has at least +# MIN_NPIX pixels within one fitting radius. + +bool procedure dp_ntmin (ids, xcen, ycen, mag, sky, npix, nier, group_size, + nterm, recenter, fitsky, groupsky, mean_sky, verbose) + +int ids[ARB] # array of ids +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real mag[ARB] # array of magnitudes +real sky[ARB] # array of sky values +int npix[ARB] # array of pixel numbers +int nier[ARB] # array of error codes +int group_size # size of the group +int nterm # number of terms +int recenter # recenter the objects +int fitsky # fit the sky value +int groupsky # use group sky value +real mean_sky # the current mean sky value +int verbose # verbose flag + +bool redo +int i + +begin + redo = false + do i = 1, group_size { + if (npix[i] >= MIN_NPIX) + next + redo = true + if (verbose == YES) { + call printf ( + "\tStar %-5d has been deleted: too few good pixels\n") + call pargi (ids[i]) + } + if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) + mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1) + call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, + nier, NSTERR_NOPIX) + if (group_size <= 0) + return (redo) + if (recenter == YES) + nterm = 3 * group_size + else + nterm = group_size + if (fitsky == YES) + nterm = nterm + 1 + } + + return (redo) +end + + +# DP_MREFLECT -- Reflect the normal matrix around the diagonal. + +procedure dp_mreflect (c, cdimen, nterm) + +real c[cdimen,ARB] # coefficient matrix +int cdimen # dimension of the c matrix +int nterm # number of terms + +int l, k + +begin + # Reflect the normal matrix across the diagonal. + do l = 2, nterm { + do k = 1, l - 1 + c[k,l] = c[l,k] + } +end + + +# DP_NTCHI -- Compute the chi value for each star. + +procedure dp_ntchi (sumwt, chi, group_size, ntmin, chigrp, grpwt) + +real sumwt[ARB] # sum of the weights +real chi[ARB] # the chis:wq +int group_size # size of the group +int ntmin # minimum number of points +real chigrp # the group chi +real grpwt # the group weight + +int i + +begin + do i = 1, group_size { + if (sumwt[i] > ntmin) { + chi[i] = CHI_NORM * chi[i] / sqrt ((sumwt[i] - ntmin) * + sumwt[i]) + sumwt[i] = ((sumwt[i] - ntmin) * chi[i] + MIN_SUMWT) / + sumwt[i] + } else { + chi[i] = chigrp + sumwt[i] = grpwt + } + } +end + + +# DP_NCHECKC -- Check the inverted matrix for singularity. + +bool procedure dp_ncheckc (c, cdimen, nterm, ids, xcen, ycen, mag, sky, + nier, group_size, recenter, fitsky, groupsky, mean_sky, verbose) + +real c[cdimen,ARB] # coefficient matrix +int cdimen # dimension of the c matrix +int nterm # number of terms +int ids[ARB] # array of ids +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real mag[ARB] # array of magnitudes +real sky[ARB] # array of sky values +int nier[ARB] # array of error codes +int group_size # size of the group +int recenter # recenter the objects +int fitsky # fit the sky value +int groupsky # use group sky value +real mean_sky # the current mean sky value +int verbose # verbose flag + +bool redo +int j, starno, i + +begin + redo = false + do j = 1, nterm { + if (c[j,j] > 0.0) + next + redo = true + if ((j == nterm) && (fitsky == YES)) + starno = 0 + else if (recenter == YES) + starno = (j + 2) / 3 + else + starno = j + if (starno == 0) { + if (verbose == YES) { + do i = 1, group_size { + call printf ( + "\tStar %-5d has been deleted: singular matrix\n") + call pargi (ids[i]) + } + } + call amovki (NSTERR_SINGULAR, nier, group_size) + group_size = 0 + } else { + if (verbose == YES) { + call printf ( + "\tStar %-5d has been deleted: singular matrix\n") + call pargi (ids[starno]) + } + if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) + mean_sky = (mean_sky * group_size - sky[starno]) / + (group_size - 1) + call dp_remove (starno, group_size, ids, xcen, ycen, mag, + sky, nier, NSTERR_SINGULAR) + } + if (group_size <= 0) + return (redo) + if (recenter == YES) + nterm = 3 * group_size + else + nterm = group_size + if (fitsky == YES) + nterm = nterm + 1 + } + + return (redo) +end + + +# DP_NTCLAMP -- Restrict the amount the solution can vary on each iteration. + +procedure dp_ntclamp (xcen, ycen, mag, magerr, sumwt, group_size, c, cdimen, + x, xold, clamp, recenter, clip, redo) + +real xcen[ARB] # x centers array +real ycen[ARB] # y centers array +real mag[ARB] # magnitude array +real magerr[ARB] # magnitude errors array +real sumwt[ARB] # array of weight sums +int group_size # size of the group +real c[cdimen, ARB] # coefficient matrix +int cdimen # dimensions of c +real x[ARB] # x vector +real xold[ARB] # old x vector +real clamp[ARB] # clamp on solution matrix +int recenter # recenter the objects +bool clip # clip the matrix +bool redo # redo the solution + +int i, l, j, k +real df + +begin + do i = 1, group_size { + + + # If any correction has changed sign since the last + # iteration, reduce the maximum permissible change by + # a factor of two. + + # Note that the sign of the correction is such that it + # must be SUBTRACTED from the current value of the + # parameter to get the improved parameter value. This being + # the case, if the correction to the brightness is + # negative (the least-squares thinks that the star should + # be brighter) a change of 1 magnitude is a change of a + # factor of 2.5; if the brightness correction is positive + # (the star should be fainter) a change of 1 magnitude + # is a change of 60%. + + if (recenter == YES) { + + l = 3 * i + k = l - 1 + j = l - 2 + + if ((xold[j] * x[j]) < 0.0) + clamp[j] = 0.5 * clamp[j] + mag[i] = mag[i] - x[j] / (1.0 + max (x[j] / + (MAX_DELTA_FAINTER * mag[i]), -x[j] / (MAX_DELTA_BRIGHTER * + mag[i])) / clamp[j]) + xold[j] = x[j] + + if ((xold[k] * x[k]) < 0.0) + clamp[k] = 0.5 * clamp[k] + if ((xold[l] * x[l]) < 0.0) + clamp[l] = 0.5 * clamp[l] + xcen[i] = xcen[i] - x[k] / (1.0 + abs(x[k]) / (clamp[k] * + MAX_DELTA_PIX)) + ycen[i] = ycen[i] - x[l] / (1.0 + abs(x[l]) / (clamp[l] * + MAX_DELTA_PIX)) + xold[k] = x[k] + xold[l] = x[l] + magerr[i] =sumwt[i] * sqrt (c[j,j]) + + } else { + + if ((xold[i] * x[i]) < 0.0) + clamp[i] = 0.5 * clamp[i] + mag[i] = mag[i] - x[i] / (1.0 + max (x[i] / + (MAX_DELTA_FAINTER * mag[i]), -x[i] / (MAX_DELTA_BRIGHTER * + mag[i])) / clamp[i]) + xold[i] = x[i] + magerr[i] =sumwt[i] * sqrt (c[i,i]) + } + + + # There are two milestones in the convergence process: the fits + # proceed normally until each star's magnitude changes by less + # than its standard error or MAX_NEW_ERRMAG magnitudes, whichever + # is greater, and its x and y centroids change by less than 0.02 + # pixel. At this point the least squares begins to apply + # down-weighting of pixels with large residuals as described + # above. The fits then continue until each star's + # magnitude changes by less than MAX (MAX_NEW_ERRMAG * std. error, + # MAX_NEW_RELBRIGHT2 magnitude), ad its centroids change by + # less than 0.002 pixel. + + if (redo) + next + + if (clip) { + if (abs (x[j]) > max (MAX_NEW_ERRMAG * magerr[i], + MAX_NEW_RELBRIGHT2 * mag[i])) { + redo = true + } else if (recenter == YES) { + df = (MAX_NEW_ERRMAG * sumwt[i]) ** 2 + if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR2)) + redo = true + else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR2)) + redo = true + } + } else { + if (abs (x[j]) > max (magerr[i], MAX_NEW_RELBRIGHT1 * + mag[i])) { + redo = true + } else if (recenter == YES) { + df = sumwt[i] ** 2 + if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR1)) + redo = true + else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR1)) + redo = true + } + } + } +end + + +# DP_NTCENTROID -- Check the new centroids to see if they have moved too +# far off the edge of the image. + +procedure dp_ntcentroid (ids, xcen, ycen, mag, sky, nier, group_size, ixmin, + ixmax, iymin, iymax, fitradsq, fitsky, groupsky, mean_sky, redo, + verbose) + +int ids[ARB] # array of ids +real xcen[ARB] # array of x centers +real ycen[ARB] # array of y centers +real mag[ARB] # array of magnitudes +real sky[ARB] # array of sky values +int nier[ARB] # array of error codes +int group_size # size of the group +int ixmin,ixmax # subraster x limits +int iymin,iymax # subraster y limits +real fitradsq # fit radius squared +int fitsky # fit the sky value +int groupsky # use the group sky value +real mean_sky # the mean sky value +bool redo # redo fit +int verbose # verbose mode + +int i +real dx, dy + +begin + # Check whether the centroid of any star has moved so far outside + # the picture that it has fewer than four or five pixels within + # one fitting radius. + + do i = 1, group_size { + + # If the centroid of the star is outside the picture in x or + # y, then DX or DY is its distance from the center of the edge + # pixel; otherwise DX and DY are zero. + + dx = max (ixmin - xcen[i], xcen[i] - ixmax, 0.0) + dy = max (iymin - ycen[i], ycen[i] - iymax, 0.0) + if ((dx <= MAX_PIX_INCREMENT) && (dy <= MAX_PIX_INCREMENT)) + next + if (((dx + 1.0) ** 2 + (dy + 1.0) ** 2) < fitradsq) + next + + # Print a warning message about the star. + if (verbose == YES) { + call printf ( + "\tStar %-5d has been deleted: new center too far off image\n") + call pargi (i) + } + + # Adjust the sky. + if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) + mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1) + + # Delete it. + call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, nier, + NSTERR_OFFIMAGE) + if (group_size < 1) + break + redo = true + } +end + + +# DP_NTFMAG -- Check for faint stars. + +procedure dp_ntfmag (mag, group_size, faint, ifaint) + +real mag[ARB] # array of magnitudes +int group_size # size of the group +real faint # faintest magnitude +int ifaint # index of faintest magnitude + +int i + +begin + faint = 1.0 + ifaint = 0 + do i = 1, group_size { + if (mag[i] > MIN_REL_BRIGHT) + next + if (mag[i] <= faint) { + faint = mag[i] + ifaint = i + } + mag[i] = MIN_REL_BRIGHT + } +end + + +# DP_FSNOISE -- Compute the smallest signal to noise ratio. + +procedure dp_fsnoise (mag, magerr, group_size, faint, ifaint) + +real mag[ARB] # array of magnitudes +real magerr[ARB] # array of magnitude errors +int group_size # size of group +real faint # faint value +int ifaint # faint index + +int i +real weight + +begin + faint = 0.0 + ifaint = 0 + do i = 1, group_size { + weight = magerr[i] / mag[i] + if (weight < faint) + next + faint = weight + ifaint = i + } +end + + +# DP_REMOVE -- Remove the i-th star from the list of stars in the current +# group by moving it to the end of the group. + +procedure dp_remove (i, nstar, ids, xcen, ycen, mag, sky, nier, pier) + +int i # the star to be removed +int nstar # the number of stars in the group +int ids[ARB] # the array of star ids +real xcen[ARB] # the array of star x positions +real ycen[ARB] # the array of star y positions +real mag[ARB] # the array of magnitudes +real sky[ARB] # the array of sky values. +int nier[ARB] # array of error codes +int pier # error code for deleted star + +int ihold, phold +real xhold, yhold, shold, mhold + +begin + nier[i] = pier + if (i != nstar) { + + ihold = ids[i] + xhold = xcen[i] + yhold = ycen[i] + shold = sky[i] + mhold = mag[i] + phold = nier[i] + + ids[i] = ids[nstar] + xcen[i] = xcen[nstar] + ycen[i] = ycen[nstar] + sky[i] = sky[nstar] + mag[i] = mag[nstar] + nier[i] = nier[nstar] + + ids[nstar] = ihold + xcen[nstar] = xhold + ycen[nstar] = yhold + sky[nstar] = shold + mag[nstar] = mhold + nier[nstar] = phold + } + nstar = nstar - 1 +end diff --git a/noao/digiphot/daophot/nstar/dpntwrite.x b/noao/digiphot/daophot/nstar/dpntwrite.x new file mode 100644 index 00000000..248a9f28 --- /dev/null +++ b/noao/digiphot/daophot/nstar/dpntwrite.x @@ -0,0 +1,600 @@ +include <tbset.h> +include <time.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/nstardef.h" + + +# DP_TPNEWNSTAR -- Create a new NSTAR output ST table. + +procedure dp_tpnewnstar (dao, nst, colpoint) + +pointer dao # pointer to the daophot structure +pointer nst # pointer to output photometry file +pointer colpoint[ARB] # array of column pointers + + +int i +pointer sp, colnames, colunits, colformat, col_dtype, col_len + +begin + # Allocate space for table definition. + call smark (sp) + call salloc (colnames, NST_NOUTCOL * (SZ_COLNAME + 1), TY_CHAR) + call salloc (colunits, NST_NOUTCOL * (SZ_COLUNITS + 1), TY_CHAR) + call salloc (colformat, NST_NOUTCOL * (SZ_COLFMT + 1), TY_CHAR) + call salloc (col_dtype, NST_NOUTCOL, TY_INT) + call salloc (col_len, NST_NOUTCOL, TY_INT) + + # Set up the column definitions. + call strcpy (ID, Memc[colnames], SZ_COLNAME) + call strcpy (GROUP, 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 (MAGERR, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME) + call strcpy (SKY, Memc[colnames+6*SZ_COLNAME+6], SZ_COLNAME) + call strcpy (NITER, Memc[colnames+7*SZ_COLNAME+7], SZ_COLNAME) + call strcpy (SHARP, Memc[colnames+8*SZ_COLNAME+8], SZ_COLNAME) + call strcpy (CHI, Memc[colnames+9*SZ_COLNAME+9], SZ_COLNAME) + call strcpy (PIER, Memc[colnames+10*SZ_COLNAME+10], SZ_COLNAME) + call strcpy (PERROR, Memc[colnames+11*SZ_COLNAME+11], SZ_COLNAME) + + # Se up the column 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 ("%14.3f", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT) + call strcpy ("%15.7g", Memc[colformat+6*SZ_COLFMT + 6], SZ_COLFMT) + call strcpy ("%6d", Memc[colformat+7*SZ_COLFMT+7], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+8*SZ_COLFMT+8], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+9*SZ_COLFMT+9], SZ_COLFMT) + call strcpy ("%6d", Memc[colformat+10*SZ_COLFMT+10], SZ_COLFMT) + call strcpy ("%13s", Memc[colformat+11*SZ_COLFMT+11], SZ_COLFMT) + + # Set up the column unit definitions. + call strcpy ("NUMBER", Memc[colunits], SZ_COLUNITS) + call strcpy ("NUMBER", 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 ("MAGNITIDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS) + call strcpy ("MAGNITIDES", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS) + call strcpy ("COUNTS", Memc[colunits+6*SZ_COLUNITS+6], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+7*SZ_COLUNITS+7], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+8*SZ_COLUNITS+8], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+9*SZ_COLUNITS+9], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+10*SZ_COLUNITS+10], SZ_COLUNITS) + call strcpy ("PERRORS", Memc[colunits+11*SZ_COLUNITS+11], SZ_COLUNITS) + + # Define the column datatypes. + 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 + Memi[col_dtype+6] = TY_REAL + Memi[col_dtype+7] = TY_INT + Memi[col_dtype+8] = TY_REAL + Memi[col_dtype+9] = TY_REAL + Memi[col_dtype+10] = TY_INT + Memi[col_dtype+11] = -13 + + # Define columnlengths. + do i = 1, NST_NOUTCOL + Memi[col_len+i-1] = 1 + + # Define and create the table. + call tbcdef (nst, colpoint, Memc[colnames], Memc[colunits], + Memc[colformat], Memi[col_dtype], Memi[col_len], NST_NOUTCOL) + call tbtcre (nst) + + # Write out some header parameters. + call dp_tnstarpars (dao, nst) + + call sfree (sp) + +end + + +define NST_NAME1STR "#N%4tID%10tGROUP%16tXCENTER%26tYCENTER%36tMAG%48t\ +MERR%62tMSKY%80t\\\n" +define NST_UNIT1STR "#U%4t##%10t##%16tpixels%26tpixels%36tmagnitudes%48t\ +magnitudes%62tcounts%80t\\\n" +define NST_FORMAT1STR "#F%4t%%-9d%10t%%-6d%16t%%-10.3f%26t%%-10.3f%36t\ +%%-12.3f%48t%%-14.3f%62t%%-15.7g%80t \n" + +define NST_NAME2STR "#N%12tNITER%18tSHARPNESS%30tCHI%42tPIER%48tPERROR%80t\\\n" +define NST_UNIT2STR "#U%12t##%18t##%30t##%42t##%48tperrors%80t\\\n" +define NST_FORMAT2STR "#F%12t%%-17d%18t%%-12.3f%30t%%-12.3f%42t%%-6d\ +%48t%%-13s%80t \n" + +# DP_XPNEWNSTAR -- Create a new NSTAR output text file. + +procedure dp_xpnewnstar (dao, nst) + +pointer dao # pointer to the daophot structure +int nst # the output file descriptor + + +begin + # Write out some header parameters. + call dp_xnstarpars (dao, nst) + + # Write out the header banner. + call fprintf (nst, "#\n") + call fprintf (nst, NST_NAME1STR) + call fprintf (nst, NST_UNIT1STR) + call fprintf (nst, NST_FORMAT1STR) + call fprintf (nst, "#\n") + call fprintf (nst, NST_NAME2STR) + call fprintf (nst, NST_UNIT2STR) + call fprintf (nst, NST_FORMAT2STR) + call fprintf (nst, "#\n") +end + + +# DP_XNSTARPARS -- Add parameters to the header of the output NSTAR text file. + +procedure dp_xnstarpars (dao, nst) + +pointer dao # pointer to the daophot structure +int nst # the output file descriptor + +pointer psffit, sp, outstr, date, time, comment +bool itob() +int envfind() + +begin + # Get pointers + psffit = DP_PSFFIT(dao) + + # 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) + call salloc (comment, SZ_LINE, TY_CHAR) + Memc[comment] = EOS + + # 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 (nst, "IRAF", Memc[outstr], "version", Memc[comment]) + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call dp_sparam (nst, "USER", Memc[outstr], "name", Memc[comment]) + call gethost (Memc[outstr], SZ_LINE) + call dp_sparam (nst, "HOST", Memc[outstr], "computer", Memc[comment]) + call dp_date (Memc[date], Memc[time], SZ_DATE) + call dp_sparam (nst, "DATE", Memc[date], "yyyy-mm-dd", Memc[comment]) + call dp_sparam (nst, "TIME", Memc[time], "hh:mm:ss", Memc[comment]) + call dp_sparam (nst, "PACKAGE", "daophot", "name", Memc[comment]) + call dp_sparam (nst, "TASK", "nstar", "name", Memc[comment]) + + # Write out the file names. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_FNAME) + call dp_sparam (nst, "IMAGE", Memc[outstr], "imagename", + Memc[comment]) + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_FNAME) + call dp_sparam (nst, "GRPFILE", Memc[outstr], "filename", Memc[comment]) + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_FNAME) + call dp_sparam (nst, "PSFIMAGE", Memc[outstr], "imagename", + Memc[comment]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME) + call dp_sparam (nst, "NSTARFILE", Memc[outstr], "filename", + Memc[comment]) + if (DP_OUTREJFILE(dao) == EOS) + call dp_sparam (nst, "REJFILE", "\"\"", "filename", Memc[comment]) + else { + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_FNAME) + call dp_sparam (nst, "REJFILE", Memc[outstr], "filename", + Memc[comment]) + } + + # Write out the data dependent parameters. + call dp_rparam (nst, "SCALE", DP_SCALE(dao), "units/pix", + Memc[comment]) + call dp_rparam (nst, "DATAMIN", DP_MINGDATA(dao), "counts", + Memc[comment]) + call dp_rparam (nst, "DATAMAX", DP_MAXGDATA(dao), "counts", + Memc[comment]) + call dp_rparam (nst, "GAIN", DP_PHOTADU(dao), "number", Memc[comment]) + call dp_rparam (nst, "READNOISE", DP_READNOISE(dao), "electrons", + Memc[comment]) + + # Write out the observing parameters. + call dp_sparam (nst, "OTIME", DP_OTIME(dao), "timeunit", Memc[comment]) + call dp_rparam (nst, "XAIRMASS", DP_XAIRMASS(dao), "number", + Memc[comment]) + call dp_sparam (nst, "IFILTER", DP_IFILTER(dao), "filter", + Memc[comment]) + + # Write out the fitting parameters. + call dp_bparam (nst, "RECENTER", itob (DP_RECENTER(dao)), "switch", + Memc[comment]) + call dp_bparam (nst, "FITSKY", itob (DP_FITSKY(dao)), "switch", + Memc[comment]) + call dp_bparam (nst, "GRPSKY", itob (DP_GROUPSKY(dao)), "switch", + Memc[comment]) + call dp_rparam (nst, "PSFMAG", DP_PSFMAG(psffit), "magnitude", + Memc[comment]) + call dp_rparam (nst, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", + Memc[comment]) + call dp_rparam (nst, "FITRAD", DP_SFITRAD(dao), "scaleunit", + Memc[comment]) + call dp_iparam (nst, "MAXITER", DP_MAXITER(dao), "number", + Memc[comment]) + call dp_iparam (nst, "MAXGROUP", DP_MAXGROUP(dao), "number", + Memc[comment]) + call dp_rparam (nst, "FLATERROR", DP_FLATERR(dao), "percentage", + Memc[comment]) + call dp_rparam (nst, "PROFERROR", DP_PROFERR(dao), "percentage", + Memc[comment]) + call dp_iparam (nst, "CLIPEXP", DP_CLIPEXP(dao), "number", + Memc[comment]) + call dp_rparam (nst, "CLIPRANGE", DP_CLIPRANGE(dao), "sigma", + Memc[comment]) + call dp_rparam (nst, "MERGERAD", DP_SMERGERAD(dao), "scaleunit", + Memc[comment]) + + call sfree(sp) +end + + +# DP_TNSTARPARS -- Add parameters to the header of the output NSTAR ST table. + +procedure dp_tnstarpars (dao, nst) + +pointer dao # pointer to the daophot structure +pointer nst # pointer to the output photometry table + +pointer psffit, sp, outstr, date, time +bool itob() +int envfind() + +begin + # Get pointers + psffit = DP_PSFFIT(dao) + + # 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) + + # 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 (nst, "IRAF", Memc[outstr]) + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call tbhadt (nst, "USER", Memc[outstr]) + call gethost (Memc[outstr], SZ_LINE) + call tbhadt (nst, "HOST", Memc[outstr]) + call dp_date (Memc[date], Memc[time], SZ_DATE) + call tbhadt (nst, "DATE", Memc[date]) + call tbhadt (nst, "TIME", Memc[time]) + call tbhadt (nst, "PACKAGE", "daophot") + call tbhadt (nst, "TASK", "nstar") + + # Write out the file names. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_FNAME) + call tbhadt (nst, "IMAGE", Memc[outstr]) + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_FNAME) + call tbhadt (nst, "GRPFILE", Memc[outstr]) + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_FNAME) + call tbhadt (nst, "PSFIMAGE", Memc[outstr]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME) + call tbhadt (nst, "NSTARFILE", Memc[outstr]) + if (DP_OUTREJFILE(dao) == EOS) + call tbhadt (nst, "REJFILE", "\"\"") + else { + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_FNAME) + call tbhadt (nst, "REJFILE", Memc[outstr]) + } + + # Write out the data dependent parameters. + call tbhadr (nst, "SCALE", DP_SCALE(dao)) + call tbhadr (nst, "DATAMIN", DP_MINGDATA(dao)) + call tbhadr (nst, "DATAMAX", DP_MAXGDATA(dao)) + call tbhadr (nst, "GAIN", DP_PHOTADU(dao)) + call tbhadr (nst, "READNOISE", DP_READNOISE(dao)) + + # Write out the observing parameters. + call tbhadt (nst, "OTIME", DP_OTIME(dao)) + call tbhadr (nst, "XAIRMASS", DP_XAIRMASS(dao)) + call tbhadt (nst, "IFILTER", DP_IFILTER(dao)) + + # Write out the fitting parameters. + call tbhadb (nst, "RECENTER", itob (DP_RECENTER(dao))) + call tbhadb (nst, "FITSKY", itob (DP_FITSKY(dao))) + call tbhadb (nst, "GRPSKY", itob (DP_GROUPSKY(dao))) + call tbhadr (nst, "PSFMAG", DP_PSFMAG(psffit)) + call tbhadr (nst, "PSFRAD", DP_SPSFRAD(dao)) + call tbhadr (nst, "FITRAD", DP_SFITRAD(dao)) + call tbhadi (nst, "MAXITER", DP_MAXITER(dao)) + call tbhadi (nst, "MAXGROUP", DP_MAXGROUP(dao)) + call tbhadr (nst, "FLATERROR", DP_FLATERR(dao)) + call tbhadr (nst, "PROFERROR", DP_PROFERR(dao)) + call tbhadi (nst, "CLIPEXP", DP_CLIPEXP(dao)) + call tbhadr (nst, "CLIPRANGE", DP_CLIPRANGE(dao)) + call tbhadr (nst, "MERGERAD", DP_SMERGERAD(dao)) + + call sfree(sp) +end + + +# DP_TNTWRITE -- Write out the NSTAR results to an ST table. + +procedure dp_tntwrite (dao, im, nst, rej, niter, old_size, output_row, + routput_row, colpoint) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +pointer nst # output photometry file descriptor +int rej # output rejections file descriptor +int niter # number of iterations +int old_size # original size of group +int output_row # output photometry file row number +int routput_row # output rejections file row number +pointer colpoint[ARB] # column pointer array + +int i, id, nkeep, nreject, pier, plen, iter +pointer psffit, nstar, apsel, sp, perror +real xcen, ycen, mag, errmag, sharp +int dp_gnsterr() + +begin + # Get some daophot pointers. + nstar = DP_NSTAR(dao) + apsel = DP_APSEL(dao) + psffit = DP_PSFFIT(dao) + + # Fill in the INDEFS. + nkeep = DP_NNUM(nstar) + nreject = old_size - nkeep + if (nreject > 0) { + call amovkr (INDEFR, Memr[DP_APMAG(apsel)+nkeep], nreject) + call amovkr (INDEFR, Memr[DP_APERR(apsel)+nkeep], nreject) + call amovkr (INDEFR, Memr[DP_APCHI(apsel)+nkeep], nreject) + } + + call smark (sp) + call salloc (perror, SZ_FNAME, TY_CHAR) + + # Now write out the results. + do i = 1, old_size { + + # Get the results. + id = Memi[DP_APID (apsel)+i-1] + xcen = Memr[DP_APXCEN (apsel)+i-1] + ycen = Memr[DP_APYCEN (apsel)+i-1] + if (IS_INDEFR(xcen) || IS_INDEFR(ycen)) + next + call dp_wout (dao, im, xcen, ycen, xcen, ycen, 1) + mag = Memr[DP_APMAG(apsel)+i-1] + errmag = Memr[DP_APERR(apsel)+i-1] + if (! IS_INDEFR(mag)) { + if (! IS_INDEFR(errmag)) + errmag = 1.085736 * errmag / mag + sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] * + Memr[DP_PSFPARS(psffit)+1] * Memr[DP_NNUMER(nstar)+i-1] / + (mag * DP_PSFHEIGHT(psffit) * Memr[DP_NDENOM(nstar)+i-1]) + if ((sharp < MIN_SHARP) || (sharp > MAX_SHARP)) + sharp = INDEFR + mag = DP_PSFMAG (psffit) - 2.5 * log10 (mag) + } else + sharp = INDEFR + pier = Memi[DP_NIER(nstar)+i-1] + if (pier == NSTERR_OK) + iter = niter + else + iter = 0 + plen = dp_gnsterr (pier, Memc[perror], SZ_FNAME) + + # Write the results to the standard output. + if (DP_VERBOSE (dao) == YES) { + call printf ( + "\tID: %5d XCEN: %8.2f YCEN: %8.2f MAG: %8.2f\n") + call pargi (id) + call pargr (xcen) + call pargr (ycen) + call pargr (mag) + } + + # Write the output row to the proper table. + if ((rej != NULL) && (pier != NSTERR_OK)) { + routput_row = routput_row + 1 + call tbrpti (rej, colpoint[1], id, 1, routput_row) + call tbrpti (rej, colpoint[2], DP_NGNUM(nstar), 1, routput_row) + call tbrptr (rej, colpoint[3], xcen, 1, routput_row) + call tbrptr (rej, colpoint[4], ycen, 1, routput_row) + call tbrptr (rej, colpoint[5], mag, 1, routput_row) + call tbrptr (rej, colpoint[6], errmag, 1, routput_row) + call tbrptr (rej, colpoint[7], Memr[DP_APMSKY(apsel)+i-1], + 1, routput_row) + call tbrpti (rej, colpoint[8], iter, 1, routput_row) + call tbrptr (rej, colpoint[9], sharp, 1, routput_row) + call tbrptr (rej, colpoint[10], Memr[DP_APCHI(apsel)+i-1], + 1, routput_row) + call tbrpti (rej, colpoint[11], pier, 1, routput_row) + call tbrptt (rej, colpoint[12], Memc[perror], plen, 1, + routput_row) + } else { + output_row = output_row + 1 + call tbrpti (nst, colpoint[1], id, 1, output_row) + call tbrpti (nst, colpoint[2], DP_NGNUM(nstar), 1, output_row) + call tbrptr (nst, colpoint[3], xcen, 1, output_row) + call tbrptr (nst, colpoint[4], ycen, 1, output_row) + call tbrptr (nst, colpoint[5], mag, 1, output_row) + call tbrptr (nst, colpoint[6], errmag, 1, output_row) + call tbrptr (nst, colpoint[7], Memr[DP_APMSKY(apsel)+i-1], + 1, output_row) + call tbrpti (nst, colpoint[8], iter, 1, output_row) + call tbrptr (nst, colpoint[9], sharp, 1, output_row) + call tbrptr (nst, colpoint[10], Memr[DP_APCHI(apsel)+i-1], + 1, output_row) + call tbrpti (nst, colpoint[11], pier, 1, output_row) + call tbrptt (nst, colpoint[12], Memc[perror], plen, 1, + output_row) + } + } + + call sfree (sp) +end + + +define NST_DATA1STR "%-9d%10t%-6d%16t%-10.3f%26t%-10.3f%36t%-12.3f%48t\ +%-14.3f%62t%-15.7g%80t\\\n" +define NST_DATA2STR "%12t%-6d%18t%-12.3f%30t%-12.3f%42t%-6d%48t%-13.13s%80t \n" + +# DP_XNTWRITE -- Write out the NSTAR results to a text file. + +procedure dp_xntwrite (dao, im, nst, rej, niter, old_size) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int nst # the output photometry file descriptor +int rej # the output rejections file descriptor +int niter # the number of the iteration +int old_size # old size of group + +int i, id, nkeep, nreject, pier, plen, iter +pointer nstar, psffit, apsel, sp, perror +real xcen, ycen, mag, errmag, sharp +int dp_gnsterr() + +begin + # Get some daophot pointers. + nstar = DP_NSTAR(dao) + psffit = DP_PSFFIT(dao) + apsel = DP_APSEL(dao) + + # Fill in the results for the rejected stars with INDEFS. + nkeep = DP_NNUM(nstar) + nreject = old_size - nkeep + if (nreject > 0) { + call amovkr (INDEFR, Memr[DP_APMAG(apsel)+nkeep], nreject) + call amovkr (INDEFR, Memr[DP_APERR(apsel)+nkeep], nreject) + call amovkr (INDEFR, Memr[DP_APCHI(apsel)+nkeep], nreject) + } + + call smark (sp) + call salloc (perror, SZ_FNAME, TY_CHAR) + + # Now write out the results. + do i = 1, old_size { + + # Get the results. + id = Memi[DP_APID (apsel)+i-1] + xcen = Memr[DP_APXCEN (apsel)+i-1] + ycen = Memr[DP_APYCEN (apsel)+i-1] + if (IS_INDEFR(xcen) || IS_INDEFR(ycen)) + next + call dp_wout (dao, im, xcen, ycen, xcen, ycen, 1) + mag = Memr[DP_APMAG(apsel)+i-1] + errmag = Memr[DP_APERR(apsel)+i-1] + if (! IS_INDEFR(mag)) { + if (! IS_INDEFR(errmag)) + errmag = 1.085736 * errmag / mag + sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] * + Memr[DP_PSFPARS(psffit)+1] * Memr[DP_NNUMER(nstar)+i-1] / + (mag * DP_PSFHEIGHT(psffit) * Memr[DP_NDENOM(nstar)+i-1]) + if ((sharp < MIN_SHARP) || (sharp > MAX_SHARP)) + sharp = INDEFR + mag = DP_PSFMAG (psffit) - 2.5 * log10 (mag) + } else + sharp = INDEFR + pier = Memi[DP_NIER(nstar)+i-1] + if (pier == NSTERR_OK) + iter = niter + else + iter = 0 + plen = dp_gnsterr (pier, Memc[perror], SZ_FNAME) + + # Write the results to the standard output. + if (DP_VERBOSE (dao) == YES) { + call printf ( + "\tID: %5d XCEN: %8.2f YCEN: %8.2f MAG: %8.2f\n") + call pargi (id) + call pargr (xcen) + call pargr (ycen) + call pargr (mag) + } + + # Write the results to the output file. + if ((rej != NULL) && (pier != NSTERR_OK)) { + call fprintf (rej, NST_DATA1STR) + call pargi (id) + call pargi (DP_NGNUM(nstar)) + call pargr (xcen) + call pargr (ycen) + call pargr (mag) + call pargr (errmag) + call pargr (Memr[DP_APMSKY(apsel)+i-1]) + call fprintf (rej, NST_DATA2STR) + call pargi (iter) + call pargr (sharp) + call pargr (Memr[DP_APCHI(apsel)+i-1]) + call pargi (pier) + call pargstr (Memc[perror]) + } else { + call fprintf (nst, NST_DATA1STR) + call pargi (id) + call pargi (DP_NGNUM(nstar)) + call pargr (xcen) + call pargr (ycen) + call pargr (mag) + call pargr (errmag) + call pargr (Memr[DP_APMSKY(apsel)+i-1]) + call fprintf (nst, NST_DATA2STR) + call pargi (iter) + call pargr (sharp) + call pargr (Memr[DP_APCHI(apsel)+i-1]) + call pargi (pier) + call pargstr (Memc[perror]) + } + } + + call sfree (sp) +end + + +# DP_GNSTERR -- Set the NSTAR task error code string. + +int procedure dp_gnsterr (ier, perror, maxch) + +int ier # the integer error code +char perror # the output error code string +int maxch # the maximum size of the error code string + +int plen +int gstrcpy() + +begin + switch (ier) { + case NSTERR_OK: + plen = gstrcpy ("No_error", perror, maxch) + case NSTERR_BIGGROUP: + plen = gstrcpy ("Big_group", perror, maxch) + case NSTERR_INDEFSKY: + plen = gstrcpy ("Bad_sky", perror, maxch) + case NSTERR_NOPIX: + plen = gstrcpy ("Npix_too_few", perror, maxch) + case NSTERR_SINGULAR: + plen = gstrcpy ("Singular", perror, maxch) + case NSTERR_FAINT: + plen = gstrcpy ("Too_faint", perror, maxch) + case NSTERR_MERGE: + plen = gstrcpy ("Merged", perror, maxch) + case NSTERR_OFFIMAGE: + plen = gstrcpy ("Off_image", perror, maxch) + default: + plen = gstrcpy ("No_error", perror, maxch) + } + + return (plen) +end diff --git a/noao/digiphot/daophot/nstar/mkpkg b/noao/digiphot/daophot/nstar/mkpkg new file mode 100644 index 00000000..327f3487 --- /dev/null +++ b/noao/digiphot/daophot/nstar/mkpkg @@ -0,0 +1,24 @@ +# NSTAR task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + dpggroup.x "../../lib/ptkeysdef.h" ../lib/daophotdef.h \ + ../lib/apseldef.h + dpmemnstar.x ../lib/daophotdef.h ../lib/nstardef.h + dpnconfirm.x ../lib/daophotdef.h + dpnstar.x <imhdr.h> <tbset.h> \ + <mach.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/nstardef.h + dpnstarfit.x <imhdr.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/nstardef.h \ + <mach.h> + dpntwrite.x <tbset.h> <time.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/nstardef.h + t_nstar.x <fset.h> <imhdr.h> \ + ../lib/daophotdef.h + ; diff --git a/noao/digiphot/daophot/nstar/t_nstar.x b/noao/digiphot/daophot/nstar/t_nstar.x new file mode 100644 index 00000000..8f6c542e --- /dev/null +++ b/noao/digiphot/daophot/nstar/t_nstar.x @@ -0,0 +1,308 @@ +include <fset.h> +include <imhdr.h> +include "../lib/daophotdef.h" + +# T_NSTAR -- Procedure to fit the PSF to multiple stars. + +procedure t_nstar () + +pointer image # input image descriptor +pointer groupfile # input group file descriptor +pointer psfimage # input PSF image descriptor +pointer nstarfile # output photometry file descriptor +pointer rejfile # output rejections file descriptor +int verbose # print messages +int verify # verify the critical parameters +int update # update the parameter set +int cache # cache the input image pixels + +pointer sp, outfname, im, psfim, dao, str +int imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist +int rlist, lrlist, root, grp, nst, rejfd, wcs, req_size, old_size +int buf_size, memstat +bool ap_text + +pointer immap(), tbtopn() +int strlen(), strncmp(), fnldir(), fstati(), open(), btoi() +int access(), imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb() +int fntgfnb(), 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) + + # Allocate working memory. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (groupfile, SZ_FNAME, TY_CHAR) + call salloc (psfimage, SZ_FNAME, TY_CHAR) + call salloc (nstarfile, SZ_FNAME, TY_CHAR) + call salloc (rejfile, SZ_FNAME, TY_CHAR) + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the input and output file names. + call clgstr ("image", Memc[image], SZ_FNAME) + call clgstr ("groupfile", Memc[groupfile], SZ_FNAME) + call clgstr ("psfimage", Memc[psfimage], SZ_FNAME) + call clgstr ("nstarfile", Memc[nstarfile], SZ_FNAME) + call clgstr ("rejfile", Memc[rejfile], SZ_FNAME) + + # Get the task mode parameters. + 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[groupfile], NO) + lalist = fntlenb (alist) + pimlist = imtopen (Memc[psfimage]) + lpimlist = imtlen (pimlist) + olist = fntopnb (Memc[nstarfile], NO) + lolist = fntlenb (olist) + rlist = fntopnb (Memc[rejfile], NO) + lrlist = fntlenb (rlist) + + # 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[groupfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and group 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 fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and psf file list lengths") + } + + if ((limlist != lolist) && (strncmp (Memc[nstarfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and nstar file list lengths") + } + + if ((lrlist > 0) && (limlist != lrlist) && (strncmp (Memc[rejfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and rejections file list lengths") + } + + # Open the daophot structure and get some parameters. + call dp_gppars (dao) + + # Set some parameters. + call dp_seti (dao, VERBOSE, verbose) + + # Verify and update the parameters if appropriate. + if (verify == YES) { + call dp_nconfirm (dao) + if (update == YES) + call dp_pppars (dao) + } + + # Get the wcs information. + wcs = clgwrd ("wcsin", Memc[str], SZ_FNAME, WCSINSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the input coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSIN, wcs) + wcs = clgwrd ("wcsout", Memc[str], SZ_FNAME, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSOUT, wcs) + wcs = clgwrd ("wcspsf", Memc[str], SZ_FNAME, WCSPSFSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the psf coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSPSF, wcs) + + # Initialize the photometry structure. + call dp_apsetup (dao) + + # Initialize the PSF structure. + call dp_fitsetup (dao) + + # Initialize the nstar structure. + call dp_nstarsetup (dao) + + # Loop over the images. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the input image. + 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 the input group table. + if (fntgfnb (alist, 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_inname (Memc[image], Memc[outfname], "grp", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[groupfile], Memc[outfname], SZ_FNAME) + ap_text = itob (access (Memc[outfname], 0, TEXT_FILE)) + if (ap_text) + grp = open (Memc[outfname], READ_ONLY, TEXT_FILE) + else + grp = tbtopn (Memc[outfname], READ_ONLY, 0) + call dp_sets (dao, INPHOTFILE, Memc[outfname]) + + # Open and read the PSF image. + 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) + psfim = immap (Memc[outfname], READ_ONLY, 0) + call dp_readpsf (dao, psfim) + call dp_sets (dao, PSFIMAGE, Memc[outfname]) + + # Open the output NSTAR file. If the output is DEF_DEFNAME, + # dir$default or a directory specification then the extension + # "nst" is added to the image name and a suitable version number + # is appended to the output name. + + if (fntgfnb (olist, Memc[nstarfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[nstarfile], SZ_FNAME) + root = fnldir (Memc[nstarfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[nstarfile+root], + DEF_LENDEFNAME) == 0 || (root == strlen (Memc[nstarfile]))) { + call dp_outname (Memc[image], Memc[outfname], "nst", + Memc[outfname], SZ_FNAME) + } else + call strcpy (Memc[nstarfile], Memc[outfname], SZ_FNAME) + if (DP_TEXT(dao) == YES) + nst = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + nst = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTPHOTFILE, Memc[outfname]) + + if (lrlist <= 0) { + rejfd = NULL + Memc[outfname] = EOS + } else { + if (fntgfnb (rlist, Memc[rejfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[rejfile], SZ_FNAME) + root = fnldir (Memc[rejfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[rejfile+root], + DEF_LENDEFNAME) == 0 || (root == strlen (Memc[rejfile]))) + call dp_outname (Memc[image], Memc[outfname], "nrj", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[rejfile], Memc[outfname], SZ_FNAME) + if (DP_TEXT(dao) == YES) + rejfd = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + rejfd = tbtopn (Memc[outfname], NEW_FILE, 0) + } + call dp_sets (dao, OUTREJFILE, Memc[outfname]) + + # Do the PSF fitting. + call dp_nphot (dao, im, grp, nst, rejfd, ap_text) + + # Close the input image. + call imunmap (im) + + # Close the group file. + if (ap_text) + call close (grp) + else + call tbtclo (grp) + + # Close the PSF image. + call imunmap (psfim) + + # Close the output photometry file. + if (DP_TEXT(dao) == YES) + call close (nst) + else + call tbtclo (nst) + + # Close the output rejections file. + if (rejfd != NULL) { + if (DP_TEXT(dao) == YES) + call close (rejfd) + else + call tbtclo (rejfd) + } + + # Uncache memory. + call fixmem (old_size) + + } + + # Close the image/file lists. + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + + # Close the nstar structure. + call dp_nsclose (dao) + + # Close the PSF structure. + call dp_fitclose (dao) + + # Close the photometry structure. + call dp_apclose (dao) + + # Free the daophot structure. + call dp_free (dao) + + call sfree(sp) +end |