diff options
Diffstat (limited to 'noao/digiphot/daophot/daolib')
36 files changed, 5010 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/daolib/bicubic.x b/noao/digiphot/daophot/daolib/bicubic.x new file mode 100644 index 00000000..50bec983 --- /dev/null +++ b/noao/digiphot/daophot/daolib/bicubic.x @@ -0,0 +1,47 @@ +# BICUBIC -- Perform bicubic interpolation on an array. The point at which +# the function is to be estimated lies between the second and third columns +# and second and third rows of f, at a distance of (dx, dy), from the (2,2) +# element of f. + +real procedure bicubic (f, nbox, dx, dy, dfdx, dfdy) + +real f[nbox,nbox] # input real array to be interpolated +int nbox # size of the input array +real dx, dy # point at which array is to be interpolated +real dfdx, dfdy # output derivative of the interpolant + +int jy +real c1, c2, c3, c4, temp[4], dfdxt[4], interp + +begin + # Interpolate first in x. + do jy = 1, 4 { + c1 = 0.5 * (f[3,jy] - f[1,jy]) + c4 = f[3,jy] - f[2,jy] - c1 + c2 = 3.0 * c4 - 0.5 * (f[4,jy] - f[2,jy]) + c1 + c3 = c4 - c2 + c4 = dx * c3 + temp[jy] = dx * (dx * (c4 + c2) + c1) + f[2,jy] + dfdxt[jy] = dx * (c4 * 3.0 + 2.0 * c2) + c1 + } + + # Interpolate next in y. + c1 = 0.5 * (temp[3] - temp[1]) + c4 = temp[3] - temp[2] - c1 + c2 = 3.0 * c4 - 0.5 * (temp[4] - temp[2]) + c1 + c3 = c4 - c2 + c4 = dy * c3 + + # Get the result. + interp = dy * (dy * (c4 + c2) + c1) + temp[2] + + # Compute the derivatives. + dfdy = dy * (c4 * 3.0 + 2.0 * c2) + c1 + c1 = 0.5 * (dfdxt[3] - dfdxt[1]) + c4 = dfdxt[3] - dfdxt[2] - c1 + c2 = 3.0 * c4 - 0.5 * (dfdxt[4] - dfdxt[2]) + c1 + c3 = c4 - c2 + dfdx = dy * (dy * (dy * c3 + c2) + c1) + dfdxt[2] + + return (interp) +end diff --git a/noao/digiphot/daophot/daolib/daoran.x b/noao/digiphot/daophot/daolib/daoran.x new file mode 100644 index 00000000..a63b8b39 --- /dev/null +++ b/noao/digiphot/daophot/daolib/daoran.x @@ -0,0 +1,43 @@ +define LEN_IR 97 +define IC 150889 +define M 714025 +define IA 1366 +define RM 1.400511e-6 + +# DAORAN -- The random number generator RAN2 from Numerical Recipes. + +real procedure daoran (idum) + +int idum # seed for the random number generator + +int j, iff, iy, ir[LEN_IR] +real rnum2 +data iff /0/ + +begin + repeat { + + # Initialize the random number generator. + if ((idum < 0) || (iff == 0)) { + iff = 1 + idum = mod (abs (IC - idum), M) + do j = 1, LEN_IR { + idum = mod (IA * idum + IC, M) + ir[j] = idum + } + idum = mod (IA * idum + IC, M) + iy = idum + } + + # Get the random number + j = 1 + (LEN_IR * iy) / M + j = max (1, min (LEN_IR, j)) + iy = ir[j] + rnum2 = iy * RM + idum = mod (IA * idum + IC, M) + ir[j] = idum + + } until (rnum2 > 0.0) + + return (rnum2) +end diff --git a/noao/digiphot/daophot/daolib/dpairmass.x b/noao/digiphot/daophot/daolib/dpairmass.x new file mode 100644 index 00000000..bad13607 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpairmass.x @@ -0,0 +1,42 @@ +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_AIRMASS -- Set the image airmass. + +procedure dp_airmass (im, dao) + +pointer im # pointer to IRAF image +pointer dao # pointer to the daophot structure + +pointer sp, key +real xair +real imgetr(), dp_statr() + +begin + # Get the airmass keyword. + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call dp_stats (dao, AIRMASS, Memc[key], SZ_FNAME) + + # Get the value. + if (Memc[key] == EOS) + xair = dp_statr (dao, XAIRMASS) + else { + iferr { + xair = imgetr (im, Memc[key]) + } then { + xair = dp_statr (dao, XAIRMASS) + call eprintf ("Warning: Image %s Keyword: %s not found\n") + call pargstr (IM_HDRFILE(im)) + call pargstr (Memc[key]) + } + } + + # Store the value. + if (IS_INDEFR(xair) || xair <= 0.0) + call dp_setr (dao, XAIRMASS, INDEFR) + else + call dp_setr (dao, XAIRMASS, xair) + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpapheader.x b/noao/digiphot/daophot/daolib/dpapheader.x new file mode 100644 index 00000000..898f9dab --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpapheader.x @@ -0,0 +1,56 @@ +# DP_APHEADER -- Copy the text database column headers to another file. +# Consider placing this simple routine in the pttables library at some point. + +procedure dp_apheader (in, out) + +int in # input file descriptor +int out # output file descriptor + +pointer sp, line +int getline() + +begin + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + + while (getline (in, Memc[line]) != EOF) { + if (Memc[line] != '#') + break + if (Memc[line+1] == 'N') + break + call putline (out, Memc[line]) + } + + call seek (in, BOF) + + call sfree (sp) +end + + +# DP_APBANNER -- Copy the text database keyword definitions to another file. +# Consider placing this simple routine in the pttables library at some point. + + +procedure dp_apbanner (in, out) + +int in # input file descriptor +int out # output file descriptor + +pointer sp, line +int getline() + +begin + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + + while (getline (in, Memc[line]) != EOF) { + if (Memc[line] != '#') + break + if (Memc[line+1] == 'K') + next + call putline (out, Memc[line]) + } + call seek (in, BOF) + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpdate.x b/noao/digiphot/daophot/daolib/dpdate.x new file mode 100644 index 00000000..70bf8f41 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpdate.x @@ -0,0 +1,28 @@ +include <time.h> + +# DP_DATE -- Create the date and time strings for the daophot output files. + +procedure dp_date (date, time, maxch) + +char date[ARB] # the date string +char time[ARB] # the time string +int maxch # the maximum number of character in the string + +int tm[LEN_TMSTRUCT] +long ctime +long clktime() +#long lsttogmt() + +begin + ctime = clktime (long(0)) + #ctime = lsttogmt (ctime) + call brktime (ctime, tm) + call sprintf (date, maxch, "%04d-%02d-%02d") + call pargi (TM_YEAR(tm)) + call pargi (TM_MONTH(tm)) + call pargi (TM_MDAY(tm)) + call sprintf (time, maxch, "%02d:%02d:%02d") + call pargi (TM_HOUR(tm)) + call pargi (TM_MIN(tm)) + call pargi (TM_SEC(tm)) +end diff --git a/noao/digiphot/daophot/daolib/dpfilter.x b/noao/digiphot/daophot/daolib/dpfilter.x new file mode 100644 index 00000000..b6932aa1 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpfilter.x @@ -0,0 +1,41 @@ +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_FILTER -- Procedure to set the image airmass. + +procedure dp_filter (im, dao) + +pointer im # pointer to IRAF image +pointer dao # pointer to the daophot structure + +pointer sp, key, filt + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (filt, SZ_FNAME, TY_CHAR) + + call dp_stats (dao, FILTER, Memc[key], SZ_FNAME) + Memc[filt] = EOS + if (Memc[key] == EOS) + call dp_stats (dao, IFILTER, Memc[filt], SZ_FNAME) + else { + iferr { + call imgstr (im, Memc[key], Memc[filt], SZ_FNAME) + } then { + call dp_stats (dao, IFILTER, Memc[filt], SZ_FNAME) + call eprintf ("Warning: Image %s Keyword: %s not found\n") + call pargstr (IM_HDRFILE(im)) + call pargstr (Memc[key]) + } + } + + if (Memc[filt] == EOS) { + call dp_sets (dao, IFILTER, "INDEF") + } else { + call dp_rmwhite (Memc[filt], Memc[filt], SZ_FNAME) + call dp_sets (dao, IFILTER, Memc[filt]) + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpfree.x b/noao/digiphot/daophot/daolib/dpfree.x new file mode 100644 index 00000000..fa7c954d --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpfree.x @@ -0,0 +1,71 @@ +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/allstardef.h" + +# DP_FREE - Procedure to free the the daophot structure. + +procedure dp_free (dp) + +pointer dp # pointer to the daophot structure + +begin + if (DP_MW(dp) != NULL) + call mw_close (DP_MW(dp)) + call mfree (dp, TY_STRUCT) +end + + +# DP_FITCLOSE -- Procedure to close up the psf fitting structure. + +procedure dp_fitclose (dp) + +pointer dp # pointer to the daophot structure + +pointer psffit + +begin + psffit = DP_PSFFIT(dp) + if (DP_PSFLUT(psffit) != NULL) + call mfree (DP_PSFLUT(psffit), TY_REAL) + if (DP_PSFPARS(psffit) != NULL) + call mfree (DP_PSFPARS(psffit), TY_REAL) + call mfree (psffit, TY_STRUCT) +end + + +# DP_APCLOSE -- Procedure to close up the APSEL parameters. + +procedure dp_apclose (dp) + +pointer dp # pointer to daophot structure + +pointer apsel + +begin + apsel = DP_APSEL(dp) + + if (DP_APRESULT(apsel) != NULL) + call mfree (DP_APRESULT(apsel), TY_INT) + if (DP_APID(apsel) != NULL) + call mfree (DP_APID(apsel), TY_INT) + if (DP_APXCEN(apsel) != NULL) + call mfree (DP_APXCEN(apsel), TY_REAL) + if (DP_APYCEN(apsel) != NULL) + call mfree (DP_APYCEN(apsel), TY_REAL) + if (DP_APMAG(apsel) != NULL) + call mfree (DP_APMAG(apsel), TY_REAL) + if (DP_APERR(apsel) != NULL) + call mfree (DP_APERR(apsel), TY_REAL) + if (DP_APMSKY(apsel) != NULL) + call mfree (DP_APMSKY(apsel), TY_REAL) + if (DP_APGROUP(apsel) != NULL) + call mfree (DP_APGROUP(apsel), TY_INT) + if (DP_APNITER(apsel) != NULL) + call mfree (DP_APNITER(apsel), TY_INT) + if (DP_APSHARP(apsel) != NULL) + call mfree (DP_APSHARP(apsel), TY_REAL) + if (DP_APCHI(apsel) != NULL) + call mfree (DP_APCHI(apsel), TY_REAL) + + call mfree (apsel, TY_STRUCT) +end diff --git a/noao/digiphot/daophot/daolib/dpgetapert.x b/noao/digiphot/daophot/daolib/dpgetapert.x new file mode 100644 index 00000000..8268373f --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpgetapert.x @@ -0,0 +1,530 @@ +include <tbset.h> +include "../../lib/ptkeysdef.h" +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +# DP_WGETAPERT -- Read the aperture photometry results and transform the +# input coordinates to the appropriate coordinate system. Works with +# either the "old" APPHOT files or the new ST Tables. + +procedure dp_wgetapert (dao, im, apd, max_nstars, old_ap) + +pointer dao # pointer to the DAOPHOT structure +pointer im # the input image descriptor +int apd # input photometry file descriptor +int max_nstars # maximum number of stars +bool old_ap # YES indicates old APPHOT file + +pointer apsel +int dp_stati() + +begin + # Get the stars. + call dp_getapert (dao, apd, max_nstars, old_ap) + + # Transform the coordinates if necessary. + apsel = DP_APSEL (dao) + if (dp_stati (dao, WCSIN) != WCS_LOGICAL) + call dp_win (dao, im, Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], DP_APNUM(apsel)) +end + + +# DP_GETAPERT -- Read the aperture photometry results. Works with +# either the "old" APPHOT files or the new ST Tables. + +procedure dp_getapert (dao, apd, max_nstars, old_ap) + +pointer dao # pointer to the DAOPHOT structure +int apd # input Photometry file descriptor +int max_nstars # maximum number of stars +bool old_ap # YES indicates old APPHOT file + +int nstars +pointer apsel +int tbpsta(), dp_goldap(), dp_gtabphot() + +begin + # Get APSEL pointer. + apsel = DP_APSEL (dao) + + # Get the required memory. + Memi[DP_APRESULT(apsel)] = DP_PAPID + Memi[DP_APRESULT(apsel)+1] = DP_PAPXCEN + Memi[DP_APRESULT(apsel)+2] = DP_PAPYCEN + Memi[DP_APRESULT(apsel)+3] = DP_PAPMAG1 + Memi[DP_APRESULT(apsel)+4] = DP_PAPSKY + + # Get the results. + if (old_ap) { + call dp_memapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT, + max_nstars) + nstars = dp_goldap (apd, dao, max_nstars) + } else { + call dp_memapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT, + tbpsta (apd, TBL_NROWS)) + nstars = dp_gtabphot (apd, dao, tbpsta (apd, TBL_NROWS)) + } + + # Reallocate to save space if appropopriate. + if (nstars < max_nstars) + call dp_rmemapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT, + nstars) +end + + +# DP_MEMAPSEL -- Procedure to allocate memory for the apselect strucuture. + +procedure dp_memapsel (dao, fields, nfields, max_nstars) + +pointer dao # pointer to the daophot strucuture +int fields[ARB] # array of fields +int nfields # number of fields to allocate space for +int max_nstars # maximum number of stars + +int i +pointer apsel + +begin + apsel = DP_APSEL(dao) + + # Allocate space for results. + do i = 1, nfields { + switch (fields[i]) { + case DP_PAPID: + if (DP_APID(apsel) != NULL) + call mfree (DP_APID(apsel), TY_INT) + call malloc (DP_APID(apsel), max_nstars, TY_INT) + + case DP_PAPXCEN: + if (DP_APXCEN(apsel) != NULL) + call mfree (DP_APXCEN(apsel), TY_REAL) + call malloc (DP_APXCEN(apsel), max_nstars, TY_REAL) + + case DP_PAPYCEN: + if (DP_APYCEN(apsel) != NULL) + call mfree (DP_APYCEN(apsel), TY_REAL) + call malloc (DP_APYCEN(apsel), max_nstars, TY_REAL) + + case DP_PAPSKY: + if (DP_APMSKY(apsel) != NULL) + call mfree (DP_APMSKY(apsel), TY_REAL) + call malloc (DP_APMSKY(apsel), max_nstars, TY_REAL) + + case DP_PAPMAG1: + if (DP_APMAG(apsel) != NULL) + call mfree (DP_APMAG(apsel), TY_REAL) + call malloc (DP_APMAG(apsel), max_nstars, TY_REAL) + + case DP_PAPGROUP: + if (DP_APGROUP(apsel) != NULL) + call mfree (DP_APGROUP(apsel), TY_INT) + #call malloc (DP_APGROUP(apsel), max_nstars, TY_INT) + + case DP_PAPMERR1: + if (DP_APERR(apsel) != NULL) + call mfree (DP_APERR(apsel), TY_REAL) + call malloc (DP_APERR(apsel), max_nstars, TY_REAL) + + case DP_PAPNITER: + if (DP_APNITER(apsel) != NULL) + call mfree (DP_APNITER(apsel), TY_INT) + #call malloc (DP_APNITER(apsel), max_nstars, TY_INT) + + case DP_PAPCHI: + if (DP_APCHI(apsel) != NULL) + call mfree (DP_APCHI(apsel), TY_REAL) + call malloc (DP_APCHI(apsel), max_nstars, TY_REAL) + + case DP_PAPSHARP: + if (DP_APSHARP(apsel) != NULL) + call mfree (DP_APSHARP(apsel), TY_REAL) + call malloc (DP_APSHARP(apsel), max_nstars, TY_REAL) + } + } +end + + +# DP_GOLDAP -- Read in the photometry from an old style APPHOT file + +int procedure dp_goldap (apd, dao, max_nstars) + +int apd # the input file descriptor +pointer dao # pointer to the daophot structure +int max_nstars # maximum number of stars + +int nstars, bufsize, stat +pointer apsel, apkey, sp, fields +int dp_apsel() + +begin + # Define the point to the apselect structure. + apsel = DP_APSEL (dao) + + # Allocate some temporary space. + call smark (sp) + call salloc (fields, SZ_LINE, TY_CHAR) + + # Initialize the keyword structure. + call pt_kyinit (apkey) + + # Set up the fields to be retrieved. + call dp_gappsf (Memi[DP_APRESULT(apsel)], Memc[fields], NAPRESULT) + + # Now read in the results. + nstars = 0 + bufsize = max_nstars + repeat { + + # Read in a group of stars. + while (nstars < bufsize) { + stat = dp_apsel (apkey, apd, Memc[fields], + Memi[DP_APRESULT(apsel)], Memi[DP_APID(apsel)+nstars], + Memr[DP_APXCEN(apsel)+nstars], + Memr[DP_APYCEN(apsel)+nstars], + Memr[DP_APMSKY(apsel)+nstars], + Memr[DP_APMAG(apsel)+nstars]) + if (stat == EOF) + break + nstars = nstars + 1 + } + + # Check the buffer size. + if (stat == EOF) + break + bufsize = bufsize + max_nstars + call dp_rmemapsel (dao, Memi[DP_APRESULT(apsel)], NAPRESULT, + bufsize) + + } + DP_APNUM(apsel) = nstars + + # Free the keyword structure. + call pt_kyfree (apkey) + call sfree (sp) + + return (nstars) +end + + +# DP_GTABPHOT -- Read in the complete photometry from an ST table. + +int procedure dp_gtabphot (tp, dao, max_nstars) + +pointer tp # table descriptor +pointer dao # pointer to daophot structure +int max_nstars # maximum number of stars + +bool nullflag +int record, index, nrow +pointer apsel, idpt, xcenpt, ycenpt, magpt, skypt +int tbpsta() + +begin + # Define the point to the apselect structure. + apsel = DP_APSEL (dao) + + # Find the column pointers + call tbcfnd (tp, ID, idpt, 1) + if (idpt == NULL) + call tbcfnd (tp, "ID", idpt, 1) + if (idpt == NULL) + call printf ("Error reading ID.\n") + + call tbcfnd (tp, XCENTER, xcenpt, 1) + if (xcenpt == NULL) + call tbcfnd (tp, "XCENTER", xcenpt, 1) + if (xcenpt == NULL) + call printf ("Error reading XCENTER.\n") + + call tbcfnd (tp, YCENTER, ycenpt, 1) + if (ycenpt == NULL) + call tbcfnd (tp, "YCENTER", ycenpt, 1) + if (ycenpt == NULL) + call printf ("Error reading YCENTER.\n") + + call tbcfnd (tp, MAG, magpt, 1) + if (magpt == NULL) + call tbcfnd (tp, APMAG, magpt, 1) + if (magpt == NULL) + call printf ("Error reading MAG.\n") + + call tbcfnd (tp, SKY, skypt, 1) + if (skypt == NULL) + call tbcfnd (tp, SKY, skypt, 1) + if (skypt == NULL) + call printf ("Error reading SKY.\n") + + + # Get the results ignoring any record with ID = NULL. + nrow = min (tbpsta (tp, TBL_NROWS), max_nstars) + index = 0 + do record = 1, nrow { + + # Check the ID record. + call tbrgti (tp, idpt, Memi[DP_APID(apsel)+index], nullflag, 1, + record) + if (nullflag) + next + + # Read the remaining records. + call tbrgtr (tp, xcenpt, Memr[DP_APXCEN(apsel)+index], nullflag, 1, + record) + call tbrgtr (tp, ycenpt, Memr[DP_APYCEN(apsel)+index], nullflag, 1, + record) + call tbrgtr (tp, magpt, Memr[DP_APMAG(apsel)+index], nullflag, 1, + record) + call tbrgtr (tp, skypt, Memr[DP_APMSKY(apsel)+index], nullflag, 1, + record) + index = index + 1 + } + + DP_APNUM(apsel) = index + + return (index) +end + + +# DP_RMEMAPSEL -- Procedure to reallocate memory for the apselect strucuture. + +procedure dp_rmemapsel (dao, fields, nfields, max_nstars) + +pointer dao # pointer to the daophot strucuture +int fields[ARB] # integer fields +int nfields # number of fields +int max_nstars # maximum number of stars + +int i +pointer apsel + +begin + # Reallocate space for results. + apsel = DP_APSEL(dao) + + do i = 1, nfields { + switch (fields[i]) { + case DP_PAPID: + call realloc (DP_APID(apsel), max_nstars, TY_INT) + case DP_PAPXCEN: + call realloc (DP_APXCEN(apsel), max_nstars, TY_REAL) + case DP_PAPYCEN: + call realloc (DP_APYCEN(apsel), max_nstars, TY_REAL) + case DP_PAPSKY: + call realloc (DP_APMSKY(apsel), max_nstars, TY_REAL) + case DP_PAPGROUP: + #call realloc (DP_APGROUP(apsel), max_nstars, TY_INT) + case DP_PAPMAG1: + call realloc (DP_APMAG(apsel), max_nstars, TY_REAL) + case DP_PAPMERR1: + call realloc (DP_APERR(apsel), max_nstars, TY_REAL) + case DP_PAPNITER: + #call realloc (DP_APNITER(apsel), max_nstars, TY_INT) + case DP_PAPSHARP: + call realloc (DP_APSHARP(apsel), max_nstars, TY_REAL) + case DP_PAPCHI: + call realloc (DP_APCHI(apsel), max_nstars, TY_REAL) + } + } +end + + +# DP_GAPPSF -- Set up the structures necessary for retrieving the +# aperture phtometery results needed for the PSF task. + +procedure dp_gappsf (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 + sel_fields[1] = EOS + + 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 (APMAG) + } + } + + if (sel_fields[1] != EOS) + sel_fields[strlen(sel_fields)] = EOS +end + + +# DP_APSEL -- Select records from an apphot/daophot text file. + +int procedure dp_apsel (key, fd, fields, indices, id, 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 of fields +int id # star id 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) + } + + 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 table 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 <= 0) + break + } + + # Construct the output record by moving selected fields + # into data structures. + call dp_getap (key, indices, id, x, y, sky, mag) + record = record + 1 + first_rec = NO + + # Record is complete so exist the loop. + break + } + } + + } + + if (nchars == EOF || nselect <= 0) { + first_rec = YES + nunique = 0 + uunique = 0 + funique = 0 + nselect = 0 + call mfree (line, TY_CHAR) + return (EOF) + } else + return (record) +end + + +# DP_GETAP -- Decode the selected daophot/apphot values. + +procedure dp_getap (key, indices, id, x, y, sky, mag) + +pointer key # pointer to keys strucuture +int indices[ARB] # index array +int id # star 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 + + # Extract the appropriate field. + call amovc (Memc[kip], buffer, maxch) + buffer[maxch+1] = EOS + + # Decode the output value. + ip = 1 + switch (indices[i]) { + case DP_PAPID: + if (ctoi (buffer, ip, id) <= 0) + id = 0 + case DP_PAPXCEN: + if (ctor (buffer, ip, x) <= 0) + x = INDEFR + case DP_PAPYCEN: + if (ctor (buffer, ip, y) <= 0) + y = INDEFR + case DP_PAPSKY: + if (ctor (buffer, ip, sky) <= 0) + sky = INDEFR + case DP_PAPMAG1: + if (ctor (buffer, ip, mag) <= 0) + mag = INDEFR + default: + call printf ("Error reading the APPHOT results.\n") + } + + } +end diff --git a/noao/digiphot/daophot/daolib/dpgppars.x b/noao/digiphot/daophot/daolib/dpgppars.x new file mode 100644 index 00000000..f16fc59e --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpgppars.x @@ -0,0 +1,227 @@ +include <ctotok.h> +include "../lib/daophotdef.h" + +# DP_GPPARS -- Procedure to fetch the daophot task parameters. + +procedure dp_gppars (dao) + +pointer dao # pointer to daophot structure + +int dp, dap +pointer mp, str, tstr +real scale, fwhmpsf, psfrad, matchrad, fitrad, annulus, dannulus, mergerad + +bool clgetb(), clgpsetb() +int clgpseti(), btoi(), dp_fctdecode(), dp_strwrd() +pointer clopset() +real clgpsetr() + +begin + # Allocate working space. + call smark (mp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (tstr, SZ_FNAME, TY_CHAR) + + # Open the daophot structure. + call dp_init (dao) + + # Set the package parameter text and initialize the verbose switch. + call dp_seti (dao, TEXT, btoi (clgetb ("text"))) + call dp_seti (dao, VERBOSE, btoi (false)) + + # Open the datapars parameter set. + dp = clopset ("datapars") + + # Set the datapars parameters. + scale = clgpsetr (dp, "scale") + call dp_setr (dao, SCALE, scale) + fwhmpsf = clgpsetr (dp, "fwhmpsf") + call dp_setr (dao, SFWHMPSF, fwhmpsf) + call dp_setr (dao, MAXGDATA, clgpsetr (dp, "datamax")) + call dp_setr (dao, MINGDATA, clgpsetr (dp, "datamin")) + + # Initialize the noise parameters. + call clgpset (dp, "ccdread", Memc[str], SZ_FNAME) + call dp_sets (dao, CCDREAD, Memc[str]) + call dp_setr (dao, READNOISE, clgpsetr (dp, "readnoise")) + call clgpset (dp, "gain", Memc[str], SZ_FNAME) + call dp_sets (dao, CCDGAIN, Memc[str]) + call dp_setr (dao, PHOTADU, clgpsetr (dp, "epadu")) + + # Initialize the observing parameters. Note that whitespace + # is removed from the filter id. + call clgpset (dp, "exposure", Memc[str], SZ_FNAME) + call dp_sets (dao, EXPTIME, Memc[str]) + call dp_setr (dao, ITIME, 1.0) + call clgpset (dp, "airmass", Memc[str], SZ_FNAME) + call dp_sets (dao, AIRMASS, Memc[str]) + call dp_setr (dao, XAIRMASS, clgpsetr (dp, "xairmass")) + call clgpset (dp, "filter", Memc[str], SZ_FNAME) + call dp_sets (dao, FILTER, Memc[str]) + call clgpset (dp, "ifilter", Memc[str], SZ_FNAME) + call dp_rmwhite (Memc[str], Memc[str], SZ_FNAME) + call dp_sets (dao, IFILTER, Memc[str]) + call clgpset (dp, "obstime", Memc[str], SZ_FNAME) + call dp_sets (dao, OBSTIME, Memc[str]) + call clgpset (dp, "otime", Memc[str], SZ_FNAME) + call dp_sets (dao, OTIME, Memc[str]) + + # Close the datapars parameter set. + call clcpset (dp) + + # Open the daopars parameter set. + dap = clopset ("daopars") + + # Set the psf fitting parameters. + call clgpset (dap, "function", Memc[tstr], SZ_FNAME) + if (dp_fctdecode (Memc[tstr], Memc[str], SZ_FNAME) <= 0) + call strcpy (",gauss,", Memc[str], SZ_FNAME) + call dp_sets (dao, FUNCLIST, Memc[str]) + if (dp_strwrd (1, Memc[tstr], SZ_FNAME, Memc[str]) <= 0) + call strcpy ("gauss", Memc[tstr], SZ_FNAME) + call dp_sets (dao, FUNCTION, Memc[tstr]) + call dp_seti (dao, VARORDER, clgpseti (dap, "varorder")) + #call dp_seti (dao, FEXPAND, btoi (clgpsetb (dap, "fexpand"))) + call dp_seti (dao, FEXPAND, NO) + call dp_seti (dao, NCLEAN, clgpseti (dap, "nclean")) + call dp_seti (dao, SATURATED, btoi (clgpsetb (dap, "saturated"))) + psfrad = clgpsetr (dap, "psfrad") + call dp_setr (dao, RPSFRAD, psfrad) + call dp_setr (dao, SPSFRAD, psfrad) + matchrad = clgpsetr (dap, "matchrad") + call dp_setr (dao, SMATCHRAD, matchrad) + + # Set the fitting parameters. + fitrad = clgpsetr (dap, "fitrad") + call dp_setr (dao, SFITRAD, fitrad) + annulus = clgpsetr (dap, "sannulus") + call dp_setr (dao, SANNULUS, annulus) + dannulus = clgpsetr (dap, "wsannulus") + call dp_setr (dao, SDANNULUS, dannulus) + call dp_setr (dao, CRITSNRATIO, clgpsetr (dap, "critsnratio")) + call dp_seti (dao, MAXITER, clgpseti (dap, "maxiter")) + call dp_seti (dao, MAXGROUP, clgpseti (dap, "maxgroup")) + call dp_seti (dao, MAXNSTAR, clgpseti (dap, "maxnstar")) + call dp_seti (dao, RECENTER, btoi (clgpsetb (dap, "recenter"))) + call dp_seti (dao, FITSKY, btoi (clgpsetb (dap, "fitsky"))) + call dp_seti (dao, GROUPSKY, btoi (clgpsetb (dap, "groupsky"))) + call dp_setr (dao, FLATERR, clgpsetr (dap, "flaterr")) + call dp_setr (dao, PROFERR, clgpsetr (dap, "proferr")) + call dp_setr (dao, CLIPRANGE, clgpsetr (dap, "cliprange")) + call dp_seti (dao, CLIPEXP, clgpseti (dap, "clipexp")) + mergerad = clgpsetr (dap, "mergerad") + call dp_setr (dao, SMERGERAD, mergerad) + + # Close the daopars pset file. + call clcpset (dap) + + # Compute the fwhmpsf, psf radius, fitting radius and matching radius + # in pixels and store. + + call dp_setr (dao, FWHMPSF, fwhmpsf / scale) + call dp_setr (dao, PSFRAD, psfrad / scale) + call dp_setr (dao, MATCHRAD, matchrad / scale) + call dp_setr (dao, FITRAD, fitrad / scale) + call dp_setr (dao, ANNULUS, annulus / scale) + call dp_setr (dao, DANNULUS, dannulus / scale) + if (IS_INDEFR(mergerad)) + call dp_setr (dao, MERGERAD, INDEFR) + else + call dp_setr (dao, MERGERAD, mergerad / scale) + + call sfree (mp) +end + + +# DP_FCTDECODE -- Decode and re-encode the list of analytic functions to be +# fit in a from suitable for use by strdic. If no valid psf types are included +# in the list set the dictionary to the gaussian function. + +int procedure dp_fctdecode (instr, outstr, maxch) + +char instr[ARB] # the input list of functions +char outstr[ARB] # the output list of functions +int maxch # maximum size of the output string + +int ip, op, ntok, tok +pointer sp, token +int ctotok(), strdic(), gstrcpy, gstrcat() + +begin + call smark (sp) + call salloc (token, maxch, TY_CHAR) + + outstr[1] = ',' + outstr[2] = EOS + op = 2 + + ntok = 0 + ip = 1 + while (instr[ip] != EOS) { + tok = ctotok (instr, ip, Memc[token], maxch) + if (tok != TOK_IDENTIFIER) + next + if (strdic (Memc[token], Memc[token], maxch, FCTN_FTYPES) <= 0) + next + ntok = ntok + 1 + op = op + gstrcpy (Memc[token], outstr[op], maxch - op + 1) + op = op + gstrcat (",", outstr[op], maxch - op + 1) + } + + call sfree (sp) + + return (ntok) +end + + +# DP_STRWRD -- Search a dictionary string for a given string index number. +# This is the opposite function of strdic(), that returns the index for +# given string. The entries in the dictionary string are separated by +# a delimiter character which is the first character of the dictionary +# string. The index of the string found is returned as the function value. +# Otherwise, if there is no string for that index, a zero is returned. + +int procedure dp_strwrd (index, outstr, maxch, dict) + +int index # String index +char outstr[ARB] # Output string as found in dictionary +int maxch # Maximum length of output string +char dict[ARB] # Dictionary string + +int i, len, start, count + +int strlen() + +begin + # Clear output string + outstr[1] = EOS + + # Return if the dictionary is not long enough + if (dict[1] == EOS) + return (0) + + # Initialize counters + count = 1 + len = strlen (dict) + + # Search the dictionary string. This loop only terminates + # successfully if the index is found. Otherwise the procedure + # returns with and error condition. + for (start = 2; count < index; start = start + 1) { + if (dict[start] == dict[1]) + count = count + 1 + if (start == len) + return (0) + } + + # Extract the output string from the dictionary + for (i = start; dict[i] != EOS && dict[i] != dict[1]; i = i + 1) { + if (i - start + 1 > maxch) + break + outstr[i - start + 1] = dict[i] + } + outstr[i - start + 1] = EOS + + # Return index for output string + return (count) +end diff --git a/noao/digiphot/daophot/daolib/dpgsubrast.x b/noao/digiphot/daophot/daolib/dpgsubrast.x new file mode 100644 index 00000000..6794824b --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpgsubrast.x @@ -0,0 +1,32 @@ +include <imhdr.h> + +# DP_GSUBRAST -- Extract a sub-raster around a specified position such that +# all pixels within the specified radius are included. + +pointer procedure dp_gsubrast (im, x, y, radius, lowx, lowy, nxpix, nypix) + +pointer im # image descriptor +real x,y # position +real radius # radius +int lowx, lowy # lower boundaries of subraster +int nxpix, nypix # number of pixels in subraster + +pointer imgs2r() + +begin + # Determine the boundaries of the area. + lowx = int (x - radius) + 1 + lowy = int (y - radius) + 1 + lowx = max (1, lowx) + lowy = max (1, lowy) + + # Compute the size of the area. + nxpix = min (int (x + radius), IM_LEN(im, 1)) - lowx + 1 + nypix = min (int (y + radius), IM_LEN(im, 2)) - lowy + 1 + + # Return a pointer to the pixel buffer. + if ((nxpix < 1) || (nypix < 1)) + return (NULL) + else + return (imgs2r (im, lowx, lowx + nxpix - 1, lowy, lowy + nypix - 1)) +end diff --git a/noao/digiphot/daophot/daolib/dpgsvw.x b/noao/digiphot/daophot/daolib/dpgsvw.x new file mode 100644 index 00000000..289b7d15 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpgsvw.x @@ -0,0 +1,162 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imio.h> +include <imhdr.h> +include <math.h> + +# DP_GSWV -- Set the data window and viewport for the image display. + +procedure dp_gswv (id, image, im, max_nframes) + +pointer id # pointer to the image display graphics stream +char image # the input image name +pointer im # pointer to the input image +int max_nframes # the maximum number of display frames + +real vx1, vx2, vy1, vy2 + +begin + if (id == NULL) + return + call dp_gsview (image, im, max_nframes, vx1, vx2, vy1, vy2) + call gsview (id, vx1, vx2, vy1, vy2) + call gswind (id, 1.0, real (IM_LEN(im,1)), 1.0, real (IM_LEN(im,2))) +end + + +# DP_GSVIEW -- Map the viewport and window of the image display. + +procedure dp_gsview (image, im, max_nframes, vx1, vx2, vy1, vy2) + +char image # the input image name +pointer im # pointer to the input image +int max_nframes # the maximum number of display frames +real vx1, vx2, vy1, vy2 # the output viewport + +int i, frame, wcs_status, dim1, dim2, step1, step2 +pointer sp, rimname, frimage, frimname, frim, iw +real x1, x2, y1, y2, fx1, fx2, fy1, fy2, junkx, junky +bool streq() +pointer imd_mapframe(), iw_open() + +begin + # Allocate some memory. + call smark (sp) + call salloc (rimname, SZ_FNAME, TY_CHAR) + call salloc (frimage, SZ_FNAME, TY_CHAR) + call salloc (frimname, SZ_FNAME, TY_CHAR) + + # Get the root image name. + call imgimage (image, Memc[rimname], SZ_FNAME) + + # Loop through the defined image frames searching for the one + # which has the image loaded. + + frame = 0 + do i = 1, max_nframes { + frim = imd_mapframe (i, READ_ONLY, NO) + iw = iw_open (frim, i, Memc[frimage], SZ_FNAME, wcs_status) + call imgimage (Memc[frimage], Memc[frimname], SZ_FNAME) + if (streq (Memc[rimname], Memc[frimname])) { + frame = i + break + } else { + call iw_close (iw) + call imunmap (frim) + } + } + + # Default to current frame if the image has not been displayes? + if (frame == 0) { + call eprintf ("Warning: image %s is not loaded in the display\n") + call pargstr (Memc[rimname]) + vx1 = 0.0 + vx2 = 1.0 + vy1 = 0.0 + vy2 = 1.0 + call sfree (sp) + return + } + + # Find the beginning and end points of the requested image section. + # We already know at this point that the input logical image is + # 2-dimensional. However this 2-dimensional section may be part of + # n-dimensional image. + + # X dimension. + dim1 = IM_VMAP(im,1) + step1 = IM_VSTEP(im,dim1) + if (step1 >= 0) { + x1 = IM_VOFF(im,dim1) + 1 + x2 = x1 + IM_LEN(im,1) - 1 + } else { + x1 = IM_VOFF(im,dim1) - 1 + x2 = x1 - IM_LEN(im,1) + 1 + } + + # Y dimension. + dim2 = IM_VMAP(im,2) + step2 = IM_VSTEP(im,dim2) + if (step2 >= 0) { + y1 = IM_VOFF(im,dim2) + 1 + y2 = y1 + IM_LEN(im,2) - 1 + } else { + y1 = IM_VOFF(im,dim2) - 1 + y2 = y1 - IM_LEN(im,2) + 1 + } + + # Get the frame buffer coordinates corresponding to the lower left + # and upper right corners of the image section. + + call iw_im2fb (iw, x1, y1, fx1, fy1) + call iw_im2fb (iw, x2, y2, fx2, fy2) + if (fx1 > fx2) { + junkx = fx1 + fx1 = fx2 + fx2 = junkx + } + if (fy1 > fy2) { + junky = fy1 + fy1 = fy2 + fy2 = junky + } + + # Check that some portion of the input image is in the display. + # If not select the default viewport and window coordinates. + if (fx1 > IM_LEN(frim,1) || fx2 < 1.0 || fy1 > IM_LEN(frim,2) || + fy2 < 1.0) { + vx1 = 0.0 + vx2 = 1.0 + vy1 = 0.0 + vy2 = 1.0 + call iw_close (iw) + call imunmap (frim) + call sfree (sp) + return + } + + # Compute a new viewport and window for X. + if (fx1 >= 1.0) + vx1 = max (0.0, min (1.0, (fx1 - 0.5) / IM_LEN(frim,1))) + else + vx1 = 0.0 + if (fx2 <= IM_LEN(frim,1)) + vx2 = max (0.0, min (1.0, (fx2 + 0.5) / IM_LEN(frim,1))) + else + vx2 = 1.0 + + # Compute a new viewport and window for Y. + if (fy1 >= 1.0) + vy1 = max (0.0, min (1.0, (fy1 - 0.5) / IM_LEN(frim,2))) + else + vy1 = 0.0 + if (fy2 <= IM_LEN(frim,2)) + vy2 = max (0.0, min (1.0, (fy2 + 0.5) / IM_LEN(frim,2))) + else + vy2 = 1.0 + + # Clean up. + call iw_close (iw) + call imunmap (frim) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpimkeys.x b/noao/digiphot/daophot/daolib/dpimkeys.x new file mode 100644 index 00000000..9cfd1943 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpimkeys.x @@ -0,0 +1,71 @@ +include "../lib/daophotdef.h" + +# DP_IMKEYS - Set the image name and keyword parameters after an image +# is mapped. + +procedure dp_imkeys (dp, im) + +pointer dp # pointer to the daophot structure +pointer im # the image descriptor + +pointer mw, ct +int dp_stati() +pointer mw_openim(), mw_sctran() +errchk mw_openim(), mw_sctran() + +begin + # Set the wcs descriptors. + mw = dp_stati (dp, MW) + if (mw != NULL) + call mw_close (mw) + iferr { + mw = mw_openim (im) + } then { + call dp_seti (dp, MW, NULL) + call dp_seti (dp, CTIN, NULL) + call dp_seti (dp, CTOUT, NULL) + call dp_seti (dp, CTPSF, NULL) + } else { + call dp_seti (dp, MW, mw) + switch (dp_stati (dp, WCSIN)) { + case WCS_WORLD: + iferr (ct = mw_sctran (mw, "world", "logical", 03B)) + ct = NULL + case WCS_PHYSICAL: + iferr (ct = mw_sctran (mw, "physical", "logical", 03B)) + ct = NULL + case WCS_TV, WCS_LOGICAL: + ct = NULL + default: + ct = NULL + } + call dp_seti (dp, CTIN, ct) + switch (dp_stati (dp, WCSOUT)) { + case WCS_PHYSICAL: + iferr (ct = mw_sctran (mw, "logical", "physical", 03B)) + ct = NULL + case WCS_TV, WCS_LOGICAL: + ct = NULL + default: + ct = NULL + } + call dp_seti (dp, CTOUT, ct) + switch (dp_stati (dp, WCSPSF)) { + case WCS_PHYSICAL: + iferr (ct = mw_sctran (mw, "logical", "physical", 03B)) + ct = NULL + case WCS_TV, WCS_LOGICAL: + ct = NULL + default: + ct = NULL + } + call dp_seti (dp, CTPSF, ct) + } + + # Get the proper values from the image header if an image is defined. + call dp_padu (im, dp) + call dp_rdnoise (im, dp) + call dp_filter (im, dp) + call dp_airmass (im, dp) + call dp_otime (im, dp) +end diff --git a/noao/digiphot/daophot/daolib/dpinit.x b/noao/digiphot/daophot/daolib/dpinit.x new file mode 100644 index 00000000..ceac884b --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpinit.x @@ -0,0 +1,225 @@ +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/allstardef.h" + +# DP_INIT - Procedure to initialize the daophot structure. + +procedure dp_init (dp) + +pointer dp # pointer to the daophot structure + +begin + # Set the daophot structure. + call calloc (dp, LEN_DPSTRUCT, TY_STRUCT) + + # Initalize the output type parameters. + DP_TEXT(dp) = YES + DP_VERBOSE(dp) = YES + + # Initialize the wcs parameters. + DP_MW(dp) = NULL + DP_WCSIN(dp) = WCS_LOGICAL + DP_WCSOUT(dp) = WCS_LOGICAL + DP_WCSPSF(dp) = WCS_LOGICAL + DP_CTIN(dp) = NULL + DP_CTOUT(dp) = NULL + DP_CTPSF(dp) = NULL + + # Initialize the data depedent parameters. + DP_SCALE(dp) = DEF_SCALE + DP_SFWHMPSF(dp) = DEF_FWHMPSF + DP_FWHMPSF(dp) = DEF_FWHMPSF / DEF_SCALE + DP_MINGDATA(dp) = INDEFR + DP_MAXGDATA(dp) = INDEFR + + # Initialize the noise parameters + DP_CCDGAIN(dp) = EOS + DP_PHOTADU(dp) = INDEFR + DP_CCDREAD(dp) = EOS + DP_READNOISE(dp) = INDEFR + + # Initialize the observing parameters. + DP_EXPTIME(dp) = EOS + DP_ITIME(dp) = INDEFR + DP_AIRMASS(dp) = EOS + DP_XAIRMASS(dp) = INDEFR + DP_FILTER(dp) = EOS + DP_IFILTER(dp) = EOS + DP_OBSTIME(dp) = EOS + DP_OTIME(dp) = EOS + + # Initialize the psf fitting parameters. + call strcpy (",gauss,", DP_FUNCLIST(dp), SZ_FNAME) + call strcpy ("gauss", DP_FUNCTION(dp), SZ_FNAME) + DP_VARORDER(dp) = -1 + DP_FEXPAND(dp) = NO + DP_SATURATED(dp) = NO + DP_RPSFRAD(dp) = DEF_PSFRAD / DEF_SCALE + DP_SPSFRAD (dp) = DEF_PSFRAD / DEF_SCALE + DP_PSFRAD (dp) = DEF_PSFRAD + + # Initialize the fitting parameters. + DP_SFITRAD(dp) = DEF_FITRAD / DEF_SCALE + DP_FITRAD(dp) = DEF_FITRAD + DP_SANNULUS(dp) = DEF_ANNULUS / DEF_SCALE + DP_ANNULUS(dp) = DEF_ANNULUS + DP_SDANNULUS(dp) = DEF_DANNULUS / DEF_SCALE + DP_DANNULUS(dp) = DEF_DANNULUS + DP_MAXITER(dp) = DEF_MAXITER + DP_MAXGROUP(dp) = DEF_MAXGROUP + DP_MAXNSTAR(dp) = DEF_MAXNSTAR + DP_RECENTER(dp) = YES + DP_FITSKY(dp) = NO + DP_GROUPSKY(dp) = YES + + # Initialize the file names. + DP_INIMAGE(dp) = EOS + DP_COORDS(dp) = EOS + DP_INPHOTFILE(dp) = EOS + DP_PSFIMAGE(dp) = EOS + DP_OUTPHOTFILE(dp) = EOS + DP_OUTIMAGE(dp) = EOS + DP_OUTREJFILE(dp) = EOS + + # Initialize the substructure pointers. + DP_PSF(dp) = NULL + DP_PSFFIT(dp) = NULL + DP_GROUP(dp) = NULL + DP_PEAK(dp) = NULL + DP_NSTAR(dp) = NULL + DP_SUBSTAR(dp) = NULL + DP_ADDSTAR(dp) = NULL + DP_ALLSTAR(dp) = NULL + DP_APSEL(dp) = NULL +end + + +# DP_FITSETUP -- Setup the current PSF parameters given that the fitting +# parameters have been correctly read in. + +procedure dp_fitsetup (dp) + +pointer dp # pointer to daophot structure + +pointer psffit +bool streq() + +begin + # Set up the psf fit structure. + call malloc (DP_PSFFIT(dp), LEN_PSFFIT, TY_STRUCT) + psffit = DP_PSFFIT(dp) + call calloc (DP_PSFPARS(psffit), MAX_NFCTNPARS, TY_REAL) + + # Define the psf function. The if user entered string is "auto" + # or a list then intialize the psf function to "gauss" or the + # first function in the list respectively. + + if (streq (DP_FUNCTION(dp), "gauss")) { + DP_PSFUNCTION(psffit) = FCTN_GAUSS + } else if (streq (DP_FUNCTION(dp), "moffat25")) { + DP_PSFUNCTION(psffit) = FCTN_MOFFAT25 + } else if (streq (DP_FUNCTION(dp), "penny1")) { + DP_PSFUNCTION(psffit) = FCTN_PENNY1 + } else if (streq (DP_FUNCTION(dp), "moffat15")) { + DP_PSFUNCTION(psffit) = FCTN_MOFFAT15 + } else if (streq (DP_FUNCTION(dp), "penny2")) { + DP_PSFUNCTION(psffit) = FCTN_PENNY2 + } else if (streq (DP_FUNCTION(dp), "lorentz")) { + DP_PSFUNCTION(psffit) = FCTN_LORENTZ + } else if (streq (DP_FUNCTION(dp), "auto")) { + DP_PSFUNCTION(psffit) = FCTN_GAUSS + } else { + call error (0, "Unknown analytic PSF function") + } + + switch (DP_VARORDER(dp)) { + case -1: + DP_NVLTABLE(psffit) = 0 + case 0: + DP_NVLTABLE(psffit) = 1 + case 1: + DP_NVLTABLE(psffit) = 3 + case 2: + DP_NVLTABLE(psffit) = 6 + } + if (DP_FEXPAND(dp) == NO) + DP_NFEXTABLE(psffit) = 0 + else + DP_NFEXTABLE(psffit) = 5 + + # Set the initial values of the function parameters. + switch (DP_PSFUNCTION(psffit)) { + case FCTN_GAUSS: + DP_PSFNPARS(psffit) = 2 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + case FCTN_MOFFAT25: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + Memr[DP_PSFPARS(psffit)+3] = 2.5 + case FCTN_PENNY1: + DP_PSFNPARS(psffit) = 4 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.75 + Memr[DP_PSFPARS(psffit)+3] = 0.0 + case FCTN_MOFFAT15: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + Memr[DP_PSFPARS(psffit)+3] = 1.5 + case FCTN_PENNY2: + DP_PSFNPARS(psffit) = 5 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.75 + Memr[DP_PSFPARS(psffit)+3] = 0.0 + Memr[DP_PSFPARS(psffit)+4] = 0.0 + case FCTN_LORENTZ: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dp) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + default: + call error (0, "Unknown analytic PSF function") + } + + DP_PSFHEIGHT(psffit) = INDEFR + DP_PSFMAG(psffit) = INDEFR + DP_PSFX(psffit) = INDEFR + DP_PSFY(psffit) = INDEFR + + DP_PSFSIZE(psffit) = 0 + DP_PSFLUT(psffit) = NULL +end + + +# DP_APSELSETUP -- Procedure to set up the APSEL parameters. + +procedure dp_apselsetup (dp) + +pointer dp # pointer to daophot structure + +pointer apsel + +begin + # APSEL structure + call malloc (DP_APSEL(dp), LEN_DPAPSTRUCT, TY_STRUCT) + apsel = DP_APSEL(dp) + call malloc (DP_APRESULT(apsel), NAPPAR, TY_INT) + + # Set the default values for the apsel parameters. + DP_APID(apsel) = NULL + DP_APXCEN(apsel)= NULL + DP_APYCEN(apsel)= NULL + DP_APMAG(apsel) = NULL + DP_APERR(apsel) = NULL + DP_APMSKY(apsel)= NULL + DP_APGROUP(apsel) = NULL + DP_APNITER(apsel) = NULL + DP_APCHI(apsel) = NULL + DP_APSHARP(apsel) = NULL +end diff --git a/noao/digiphot/daophot/daolib/dpnames.x b/noao/digiphot/daophot/daolib/dpnames.x new file mode 100644 index 00000000..ca1fa858 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpnames.x @@ -0,0 +1,415 @@ + +# DP_IIMNAME -- Procedure to construct an daophot input image name. +# If input is null or a directory a name is constructed from the root +# of the image name and the extension. The disk is searched to avoid +# name collisions. + +procedure dp_iimname (image, input, ext, name, maxch) + +char image[ARB] # image name +char input[ARB] # input directory or name +char ext[ARB] # extension +char name[ARB] # input name +int maxch # maximum size of name + +int ndir, nimdir, clindex, clsize +pointer sp, root, str +int fnldir(), strlen() + +begin + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + ndir = fnldir (input, name, maxch) + if (strlen (input) == ndir) { + call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME, + Memc[str], SZ_FNAME, clindex, clsize) + nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME) + if (clindex >= 0) { + call sprintf (name[ndir+1], maxch, "%s%d.%s.*") + call pargstr (Memc[root+nimdir]) + call pargi (clindex) + call pargstr (ext) + } else { + call sprintf (name[ndir+1], maxch, "%s.%s.*") + call pargstr (Memc[root+nimdir]) + call pargstr (ext) + } + + call dp_iimversion (name, name, maxch) + } else + call strcpy (input, name, maxch) + + call sfree (sp) +end + + +# DP_IMROOT -- Fetch the root image name minus the directory specification +# and the section notation. + +procedure dp_imroot (image, root, maxch) + +char image[ARB] # image specification +char root[ARB] # output root name +int maxch # maximum number of characters + +pointer sp, imroot, kernel, section, str +int clindex, clsize, nchars +int fnldir() + +begin + call smark (sp) + call salloc (imroot, SZ_PATHNAME, TY_CHAR) + call salloc (kernel, SZ_FNAME, TY_CHAR) + call salloc (section, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_PATHNAME, TY_CHAR) + + call imparse (image, Memc[imroot], SZ_PATHNAME, Memc[kernel], SZ_FNAME, + Memc[section], SZ_FNAME, clindex, clsize) + nchars = fnldir (Memc[imroot], Memc[str], SZ_PATHNAME) + if (clindex >= 0) { + call sprintf (root, maxch, "%s[%d]%s%s") + call pargstr (Memc[imroot+nchars]) + call pargi (clindex) + call pargstr (Memc[kernel]) + call pargstr (Memc[section]) + } else { + call sprintf (root, maxch, "%s%s%s") + call pargstr (Memc[imroot+nchars]) + call pargstr (Memc[kernel]) + call pargstr (Memc[section]) + } + + call sfree (sp) +end + + +# DP_OIMVERSION -- Routine to compute the next available version number of +# a given file name template and output the new files name. + +procedure dp_oimversion (template, filename, maxch) + +char template[ARB] # name template +char filename[ARB] # output name +int maxch # maximum number of characters + +char period +int newversion, version, len +pointer sp, list, name +int imtopen(), imtgetim(), strldx(), ctoi() + +begin + # Allocate temporary space + call smark (sp) + call salloc (name, maxch, TY_CHAR) + period = '.' + list = imtopen (template) + + # Loop over the names in the list searchng for the highest version. + newversion = 0 + while (imtgetim (list, Memc[name], maxch) != EOF) { + len = strldx (period, Memc[name]) + Memc[name+len-1] = EOS + len = strldx (period, Memc[name]) + len = len + 1 + if (ctoi (Memc[name], len, version) <= 0) + next + newversion = max (newversion, version) + } + + # Make new output file name. + len = strldx (period, template) + call strcpy (template, filename, len) + call sprintf (filename[len+1], maxch, "%d") + call pargi (newversion + 1) + + call imtclose (list) + call sfree (sp) +end + + +# DP_IIMVERSION -- Routine to compute the next available version number of +# a given file name template and output the new files name. + +procedure dp_iimversion (template, filename, maxch) + +char template[ARB] # name template +char filename[ARB] # output name +int maxch # maximum number of characters + +char period +int newversion, version, len +pointer sp, list, name +int imtopen(), imtgetim() strldx(), ctoi() + +begin + # Allocate temporary space + call smark (sp) + call salloc (name, maxch, TY_CHAR) + period = '.' + list = imtopen (template) + + # Loop over the names in the list searchng for the highest version. + newversion = 0 + while (imtgetim (list, Memc[name], maxch) != EOF) { + len = strldx (period, Memc[name]) + Memc[name+len-1] = EOS + len = strldx (period, Memc[name]) + len = len + 1 + if (ctoi (Memc[name], len, version) <= 0) + next + newversion = max (newversion, version) + } + + # Make new output file name. + len = strldx (period, template) + call strcpy (template, filename, len) + call sprintf (filename[len+1], maxch, "%d") + call pargi (newversion) + + call imtclose (list) + call sfree (sp) +end + + +# DP_OIMNAME -- Procedure to construct an daophot output image name. +# If output is null or a directory a name is constructed from the root +# of the image name and the extension. The disk is searched to avoid +# name collisions. + +procedure dp_oimname (image, output, ext, name, maxch) + +char image[ARB] # image name +char output[ARB] # output directory or name +char ext[ARB] # extension +char name[ARB] # output name +int maxch # maximum size of name + +int ndir, nimdir, clindex, clsize +pointer sp, root, str +int fnldir(), strlen() + +begin + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + ndir = fnldir (output, name, maxch) + if (strlen (output) == ndir) { + call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME, + Memc[str], SZ_FNAME, clindex, clsize) + nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME) + if (clindex >= 0) { + call sprintf (name[ndir+1], maxch, "%s%d.%s.*") + call pargstr (Memc[root+nimdir]) + call pargi (clindex) + call pargstr (ext) + } else { + call sprintf (name[ndir+1], maxch, "%s.%s.*") + call pargstr (Memc[root+nimdir]) + call pargstr (ext) + } + call dp_oimversion (name, name, maxch) + } else + call strcpy (output, name, maxch) + + call sfree (sp) +end + + +# DP_INNAME -- Procedure to construct an daophot input file name. +# If input is null or a directory a name is constructed from the root +# of the image name and the extension. The disk is searched to avoid +# name collisions. + +procedure dp_inname (image, input, ext, name, maxch) + +char image[ARB] # image name +char input[ARB] # input directory or name +char ext[ARB] # extension +char name[ARB] # input name +int maxch # maximum size of name + +int ndir, nimdir, clindex, clsize +pointer sp, root, str +int fnldir(), strlen() + +begin + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + ndir = fnldir (input, name, maxch) + if (strlen (input) == ndir) { + call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME, + Memc[str], SZ_FNAME, clindex, clsize) + nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME) + if (clindex >= 0) { + call sprintf (name[ndir+1], maxch, "%s%d.%s.*") + call pargstr (Memc[root+nimdir]) + call pargi (clindex) + call pargstr (ext) + } else { + call sprintf (name[ndir+1], maxch, "%s.%s.*") + call pargstr (Memc[root+nimdir]) + call pargstr (ext) + } + call dp_iversion (name, name, maxch) + } else + call strcpy (input, name, maxch) + + call sfree (sp) +end + + +# DP_OUTNAME -- Procedure to construct an daophot output file name. +# If output is null or a directory a name is constructed from the root +# of the image name and the extension. The disk is searched to avoid +# name collisions. + +procedure dp_outname (image, output, ext, name, maxch) + +char image[ARB] # image name +char output[ARB] # output directory or name +char ext[ARB] # extension +char name[ARB] # output name +int maxch # maximum size of name + +int ndir, nimdir, clindex, clsize +pointer sp, root, str +int fnldir(), strlen() + +begin + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + ndir = fnldir (output, name, maxch) + if (strlen (output) == ndir) { + call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME, + Memc[str], SZ_FNAME, clindex, clsize) + nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME) + if (clindex >= 0) { + call sprintf (name[ndir+1], maxch, "%s%d.%s.*") + call pargstr (Memc[root+nimdir]) + call pargi (clindex) + call pargstr (ext) + } else { + call sprintf (name[ndir+1], maxch, "%s.%s.*") + call pargstr (Memc[root+nimdir]) + call pargstr (ext) + } + call dp_oversion (name, name, maxch) + } else + call strcpy (output, name, maxch) + + call sfree (sp) +end + + +# DP_OVERSION -- Routine to compute the next available version number of a given +# file name template and output the new files name. + +procedure dp_oversion (template, filename, maxch) + +char template[ARB] # name template +char filename[ARB] # output name +int maxch # maximum number of characters + +char period +int newversion, version, len +pointer sp, list, name +int fntgfnb() strldx(), ctoi(), fntopnb() + +begin + # Allocate temporary space + call smark (sp) + call salloc (name, maxch, TY_CHAR) + period = '.' + list = fntopnb (template, NO) + + # Loop over the names in the list searchng for the highest version. + newversion = 0 + while (fntgfnb (list, Memc[name], maxch) != EOF) { + len = strldx (period, Memc[name]) + len = len + 1 + if (ctoi (Memc[name], len, version) <= 0) + next + newversion = max (newversion, version) + } + + # Make new output file name. + len = strldx (period, template) + call strcpy (template, filename, len) + call sprintf (filename[len+1], maxch, "%d") + call pargi (newversion + 1) + + call fntclsb (list) + call sfree (sp) +end + + +# DP_IVERSION -- Routine to compute the last version number of a given +# file name template and output the new files name. + +procedure dp_iversion (template, filename, maxch) + +char template[ARB] # name template +char filename[ARB] # output name +int maxch # maximum number of characters + +char period +int newversion, version, len +pointer sp, list, name +int fntgfnb() strldx(), ctoi(), fntopnb() + +begin + # Allocate temporary space + call smark (sp) + call salloc (name, maxch, TY_CHAR) + period = '.' + list = fntopnb (template, NO) + + # Loop over the names in the list searchng for the highest version. + newversion = 0 + while (fntgfnb (list, Memc[name], maxch) != EOF) { + len = strldx (period, Memc[name]) + len = len + 1 + if (ctoi (Memc[name], len, version) <= 0) + next + newversion = max (newversion, version) + } + + # Make new output file name. + len = strldx (period, template) + call strcpy (template, filename, len) + call sprintf (filename[len+1], maxch, "%d") + call pargi (newversion) + + call fntclsb (list) + call sfree (sp) +end + + +# DP_FROOT -- Fetch the file name minus the directory specification, + +procedure dp_froot (filename, root, maxch) + +char filename[ARB] # input file name +char root[ARB] # output root file name +int maxch # maximum number of characters + +pointer sp, str +int nchars +int fnldir() + +begin + call smark (sp) + call salloc (str, SZ_PATHNAME, TY_CHAR) + + nchars = fnldir (filename, Memc[str], SZ_PATHNAME) + call strcpy (filename[nchars+1], root, maxch) + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpotime.x b/noao/digiphot/daophot/daolib/dpotime.x new file mode 100644 index 00000000..0c8a7925 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpotime.x @@ -0,0 +1,51 @@ +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_OTIME -- Read the epoch of the observation from the image header. + +procedure dp_otime (im, dao) + +pointer im # pointer to IRAF image +pointer dao # pointer to the daophot structure + +char timechar +int index +pointer sp, key, otime +bool streq() +int strldx() + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (otime, SZ_FNAME, TY_CHAR) + + call dp_stats (dao, OBSTIME, Memc[key], SZ_FNAME) + Memc[otime] = EOS + if (Memc[key] == EOS) + call dp_stats (dao, OTIME, Memc[otime], SZ_FNAME) + else { + iferr { + call imgstr (im, Memc[key], Memc[otime], SZ_FNAME) + } then { + call dp_stats (dao, OTIME, Memc[otime], SZ_FNAME) + call eprintf ("Warning: Image %s Keyword: %s not found\n") + call pargstr (IM_HDRFILE(im)) + call pargstr (Memc[key]) + } + } + if (Memc[otime] == EOS) { + call dp_sets (dao, OTIME, "INDEF") + } else if (streq ("DATE-OBS", Memc[key]) || streq ("date-obs", + Memc[key])) { + timechar = 'T' + index = strldx (timechar, Memc[otime]) + if (index > 0) + call dp_sets (dao, OTIME, Memc[otime+index]) + else + call dp_sets (dao, OTIME, "INDEF") + } else { + call dp_sets (dao, OTIME, Memc[otime]) + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dppadu.x b/noao/digiphot/daophot/daolib/dppadu.x new file mode 100644 index 00000000..b1ca7c25 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dppadu.x @@ -0,0 +1,36 @@ +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_PADU -- Read the gain value from the image header. + +procedure dp_padu (im, dao) + +pointer im # pointer to IRAF image +pointer dao # pointer to the daophot structure + +pointer sp, key +real padu +real imgetr(), dp_statr() + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call dp_stats (dao, CCDGAIN, Memc[key], SZ_FNAME) + if (Memc[key] == EOS) + padu = dp_statr (dao, PHOTADU) + else { + iferr { + padu = imgetr (im, Memc[key]) + } then { + padu = dp_statr (dao, PHOTADU) + call eprintf ("Warning: Image %s Keyword %s not found.\n") + call pargstr (IM_HDRFILE(im)) + call pargstr (Memc[key]) + } + } + if (IS_INDEFR(padu) || padu <= 0.0) + call dp_setr (dao, PHOTADU, 1.0) + else + call dp_setr (dao, PHOTADU, padu) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dppcache.x b/noao/digiphot/daophot/daolib/dppcache.x new file mode 100644 index 00000000..81ca7a3d --- /dev/null +++ b/noao/digiphot/daophot/daolib/dppcache.x @@ -0,0 +1,83 @@ +include <imhdr.h> +include <imset.h> + +# DP_MEMSTAT -- Figure out if there is enough memory to cache the image +# pixels. If it is necessary to request more memory and the memory is +# avalilable return YES otherwise return NO. + +int procedure dp_memstat (cache, req_size, old_size) + +int cache #I cache memory ? +int req_size #I the requested working set size in chars +int old_size #O the original working set size in chars + +int cur_size, max_size +int begmem() + +begin + # Find the default working set size. + cur_size = begmem (0, old_size, max_size) + + # If cacheing is disabled return NO regardless of the working set size. + if (cache == NO) + return (NO) + + # If the requested working set size is less than the current working + # set size return YES. + if (req_size <= cur_size) + return (YES) + + # Reset the current working set size. + cur_size = begmem (req_size, old_size, max_size) + if (req_size <= cur_size) { + return (YES) + } else { + return (NO) + } +end + + +# DP_PCACHE -- Cache the image pixels im memory by resetting the default image +# buffer size. If req_size is INDEF the size of the image is used to determine +# the size of the image i/o buffers. + +procedure dp_pcache (im, req_size, buf_size) + +pointer im #I the input image point +int req_size #I the requested working set size in chars +int buf_size #O the new image buffer size + +int def_size, new_imbufsize +int sizeof(), imstati() + +begin + # Find the default buffer size. + def_size = imstati (im, IM_BUFSIZE) + + # Return if the image is not 2-dimensional. + if (IM_NDIM(im) != 2) { + buf_size = def_size + return + } + + # Compute the new required image i/o buffer size in chars. + if (IS_INDEFI(req_size)) { + new_imbufsize = IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + } else { + new_imbufsize = req_size + } + + # If the default image i/o buffer size is already bigger than + # the requested size do nothing. + if (def_size >= new_imbufsize) { + buf_size = def_size + return + } + + # Reset the image i/o buffer. + call imseti (im, IM_BUFSIZE, new_imbufsize) + call imseti (im, IM_BUFFRAC, 0) + buf_size = new_imbufsize + return +end diff --git a/noao/digiphot/daophot/daolib/dpppars.x b/noao/digiphot/daophot/daolib/dpppars.x new file mode 100644 index 00000000..e1c86f58 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpppars.x @@ -0,0 +1,94 @@ +include "../lib/daophotdef.h" + +# DP_PPPARS -- Store the daophot package parameters in the pset files. + +procedure dp_pppars (dao) + +pointer dao # pointer to daophot structure + +pointer sp, str, fstr, dap +bool itob() +int dp_stati(), strlen() +pointer clopset() +real dp_statr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (fstr, SZ_FNAME, TY_CHAR) + + # Set the package parameter text. + call clputb ("text", itob (dp_stati (dao, TEXT))) + + # Get and set the datapars parameter sets. + dap = clopset ("datapars") + + # Store the data dependent parameters. + call clppsetr (dap, "scale", dp_statr (dao, SCALE)) + call clppsetr (dap, "fwhmpsf", dp_statr (dao, SFWHMPSF)) + call clppsetr (dap, "datamin", dp_statr (dao, MINGDATA)) + call clppsetr (dap, "datamax", dp_statr (dao, MAXGDATA)) + + # Store the noise parameters. + call dp_stats (dao, CCDGAIN, Memc[str], SZ_FNAME) + call clppset (dap, "gain", Memc[str]) + call clppsetr (dap, "epadu", dp_statr (dao, PHOTADU)) + call dp_stats (dao, CCDREAD, Memc[str], SZ_FNAME) + call clppset (dap, "ccdread", Memc[str]) + call clppsetr (dap, "readnoise", dp_statr (dao, READNOISE)) + + # Store the observing parameters. + call dp_stats (dao, EXPTIME, Memc[str], SZ_FNAME) + call clppset (dap, "exposure", Memc[str]) + call clppsetr (dap, "itime", dp_statr (dao, ITIME)) + call dp_stats (dao, AIRMASS, Memc[str], SZ_FNAME) + call clppset (dap, "airmass", Memc[str]) + call clppsetr (dap, "xairmass", dp_statr (dao, XAIRMASS)) + call dp_stats (dao, FILTER, Memc[str], SZ_FNAME) + call clppset (dap, "filter", Memc[str]) + call dp_stats (dao, IFILTER, Memc[str], SZ_FNAME) + call clppset (dap, "ifilter", Memc[str]) + call dp_stats (dao, OBSTIME, Memc[str], SZ_FNAME) + call clppset (dap, "obstime", Memc[str]) + call dp_stats (dao, OTIME, Memc[str], SZ_FNAME) + call clppset (dap, "otime", Memc[str]) + + # Close the datapars parameter set. + call clcpset (dap) + + # Open the daopars parameter set. + dap = clopset ("daopars") + + # Store the psf function parameters. + call dp_stats (dao, FUNCLIST, Memc[fstr], SZ_FNAME) + call strcpy (Memc[fstr+1], Memc[str], strlen(Memc[fstr+1]) - 1) + call clppset (dap, "function", Memc[str]) + call clppseti (dap, "varorder", dp_stati (dao, VARORDER)) + #call clppsetb (dap, "fexpand", itob (dp_stati (dao, FEXPAND))) + call clppseti (dap, "nclean", dp_stati (dao, NCLEAN)) + call clppsetb (dap, "saturated", itob (dp_stati (dao, SATURATED))) + call clppsetr (dap, "psfrad", dp_statr (dao, SPSFRAD)) + call clppsetr (dap, "matchrad", dp_statr (dao, SMATCHRAD)) + + # Store the fitting algorithm parameters. + call clppsetr (dap, "fitrad", dp_statr (dao, SFITRAD)) + call clppsetr (dap, "sannulus", dp_statr (dao, SANNULUS)) + call clppsetr (dap, "wsannulus", dp_statr (dao, SDANNULUS)) + call clppsetr (dap, "critsnratio", dp_statr (dao, CRITSNRATIO)) + call clppseti (dap, "maxiter", dp_stati (dao, MAXITER)) + call clppseti (dap, "maxgroup", dp_stati (dao, MAXGROUP)) + call clppseti (dap, "maxnstar", dp_stati (dao, MAXNSTAR)) + call clppsetb (dap, "recenter", itob (dp_stati (dao, RECENTER))) + call clppsetb (dap, "fitsky", itob (dp_stati (dao, FITSKY))) + call clppsetb (dap, "groupsky", itob (dp_stati (dao, GROUPSKY))) + call clppsetr (dap, "flaterr", dp_statr (dao, FLATERR)) + call clppsetr (dap, "proferr", dp_statr (dao, PROFERR)) + call clppsetr (dap, "cliprange", dp_statr (dao, CLIPRANGE)) + call clppseti (dap, "clipexp", dp_stati (dao, CLIPEXP)) + call clppsetr (dap, "mergerad", dp_statr (dao, SMERGERAD)) + + # Close the daopars parameter set. + call clcpset (dap) + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dprdnoise.x b/noao/digiphot/daophot/daolib/dprdnoise.x new file mode 100644 index 00000000..9e4baa3c --- /dev/null +++ b/noao/digiphot/daophot/daolib/dprdnoise.x @@ -0,0 +1,36 @@ +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_RDNOISE - Read the readout noise value from the image header. + +procedure dp_rdnoise (im, dao) + +pointer im # pointer to IRAF image +pointer dao # pointer to the daophot structure + +pointer sp, key +real rdnoise +real imgetr(), dp_statr() + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call dp_stats (dao, CCDREAD, Memc[key], SZ_FNAME) + if (Memc[key] == EOS) + rdnoise = dp_statr (dao, READNOISE) + else { + iferr { + rdnoise = imgetr (im, Memc[key]) + } then { + rdnoise = dp_statr (dao, READNOISE) + call eprintf ("Warning: Image %s Keyword %s not found.\n") + call pargstr (IM_HDRFILE(im)) + call pargstr (Memc[key]) + } + } + if (IS_INDEFR(rdnoise) || rdnoise <= 0.0) + call dp_setr (dao, READNOISE, 0.0) + else + call dp_setr (dao, READNOISE, rdnoise) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dpreadpsf.x b/noao/digiphot/daophot/daolib/dpreadpsf.x new file mode 100644 index 00000000..77aa0b74 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpreadpsf.x @@ -0,0 +1,138 @@ +include <error.h> +include <imhdr.h> +include "../lib/daophotdef.h" + +# DP_READPSF -- Read in the PSF from the specified image. + +procedure dp_readpsf (dao, im) + +pointer dao # pointer to the DAOPHOT Structure +pointer im # image descriptor + +int i, ival, npsfstars +pointer sp, str, v, psffit, psflut, buf +real scale, rval +bool imgetb(), streq() +int imgeti(), imgnlr(), btoi() +real imgetr() +errchk imgetr(), imgeti() + +begin + # Allocate working memory. + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (v, IM_MAXDIM, TY_LONG) + + # Set up needed daophot pointers. + psffit = DP_PSFFIT(dao) + + # Read in the function parameters. + call imgstr (im, "FUNCTION", Memc[str], SZ_FNAME) + if (streq (Memc[str], "gauss")) { + DP_PSFUNCTION(psffit) = FCTN_GAUSS + call strcpy ("gauss", DP_FUNCTION(dao), SZ_FNAME) + } else if (streq (Memc[str], "moffat25")) { + DP_PSFUNCTION(psffit) = FCTN_MOFFAT25 + call strcpy ("moffat25", DP_FUNCTION(dao), SZ_FNAME) + } else if (streq (Memc[str], "penny1")) { + DP_PSFUNCTION(psffit) = FCTN_PENNY1 + call strcpy ("penny1", DP_FUNCTION(dao), SZ_FNAME) + } else if (streq (Memc[str], "moffat15")) { + DP_PSFUNCTION(psffit) = FCTN_MOFFAT15 + call strcpy ("moffat15", DP_FUNCTION(dao), SZ_FNAME) + } else if (streq (Memc[str], "penny2")) { + DP_PSFUNCTION(psffit) = FCTN_PENNY2 + call strcpy ("penny2", DP_FUNCTION(dao), SZ_FNAME) + } else if (streq (Memc[str], "lorentz")) { + DP_PSFUNCTION(psffit) = FCTN_LORENTZ + call strcpy ("lorentz", DP_FUNCTION(dao), SZ_FNAME) + } else + call error (0, "Unknown PSF function in PSF image\n") + + # Read in the position and brightness parameters. + DP_PSFX (psffit) = imgetr (im, "PSFX") + DP_PSFY (psffit) = imgetr (im, "PSFY") + DP_PSFHEIGHT(psffit) = imgetr (im, "PSFHEIGHT") + DP_PSFMAG (psffit) = imgetr (im, "PSFMAG") + + DP_PSFNPARS(psffit) = imgeti (im, "NPARS") + do i = 1, DP_PSFNPARS(psffit) { + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (i) + Memr[DP_PSFPARS(psffit)+i-1] = imgetr (im, Memc[str]) + } + + # Get the psfradius with which the psf was made. Make sure the + # psf radius requested by the user is less than or equal to the + # stored psf radius. + + iferr { + scale = imgetr (im, "SCALE") + } then { + DP_PSFRAD(dao) = min (DP_RPSFRAD(dao) / DP_SCALE(dao), + imgetr (im, "PSFRAD")) + DP_SPSFRAD(dao) = DP_SCALE(dao) * DP_PSFRAD(dao) + } else { + DP_PSFRAD(dao) = min (DP_RPSFRAD(dao) / DP_SCALE(dao), + imgetr (im, "PSFRAD") / scale) + DP_SPSFRAD(dao) = DP_SCALE(dao) * DP_PSFRAD(dao) + } + + # Get the lookup table(s) parameters. + DP_VARORDER(dao) = imgeti (im, "VARORDER") + switch (DP_VARORDER(dao)) { + case -1: + DP_NVLTABLE(psffit) = 0 + case 0: + DP_NVLTABLE(psffit) = 1 + case 1: + DP_NVLTABLE(psffit) = 3 + case 2: + DP_NVLTABLE(psffit) = 6 + } + DP_FEXPAND(dao) = btoi (imgetb (im, "FEXPAND")) + if (DP_FEXPAND(dao) == NO) + DP_NFEXTABLE(psffit) = 0 + else + DP_NFEXTABLE(psffit) = 5 + + # Read in the lookup table(s). + if ((DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit)) <= 0) { + DP_PSFSIZE(psffit) = 0 + } else { + DP_PSFSIZE(psffit) = IM_LEN(im,1) + call realloc (DP_PSFLUT(psffit), DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit) * IM_LEN(im,3), TY_REAL) + psflut = DP_PSFLUT (psffit) + call amovkl (long(1), Meml[v], IM_MAXDIM) + while (imgnlr (im, buf, Meml[v]) != EOF) { + call amovr (Memr[buf], Memr[psflut], DP_PSFSIZE(psffit)) + psflut = psflut + DP_PSFSIZE(psffit) + } + } + + # Check that the complete header can be read. + iferr { + npsfstars = imgeti (im, "NPSFSTAR") + call sprintf (Memc[str], SZ_FNAME, "ID%d") + call pargi (npsfstars) + ival = imgeti (im, Memc[str]) + call sprintf (Memc[str], SZ_FNAME, "X%d") + call pargi (npsfstars) + rval = imgetr (im, Memc[str]) + call sprintf (Memc[str], SZ_FNAME, "Y%d") + call pargi (npsfstars) + rval = imgetr (im, Memc[str]) + call sprintf (Memc[str], SZ_FNAME, "MAG%d") + call pargi (npsfstars) + rval = imgetr (im, Memc[str]) + } then { + call eprintf ("PSF image header is too long to be read in.\n") + call eprintf ( + "Reset min_lenusearea environment variable and try again.") + call erract (EA_ERROR) + } + + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daolib/dprmwhite.x b/noao/digiphot/daophot/daolib/dprmwhite.x new file mode 100644 index 00000000..9c61969c --- /dev/null +++ b/noao/digiphot/daophot/daolib/dprmwhite.x @@ -0,0 +1,22 @@ +include <ctype.h> + +# DP_RMWHITE -- Remove whitespace from a string. + +procedure dp_rmwhite (instr, outstr, maxch) + +char instr[ARB] # the input string +char outstr[ARB] # the output string, may be the same as instr +int maxch # maximum number of characters in outstr + +int ip, op + +begin + op = 1 + for (ip = 1; (instr[ip] != EOS) && (op <= maxch); ip = ip + 1) { + if (IS_WHITE(instr[ip])) + next + outstr[op] = instr[ip] + op = op + 1 + } + outstr[op] = EOS +end diff --git a/noao/digiphot/daophot/daolib/dpset.x b/noao/digiphot/daophot/daolib/dpset.x new file mode 100644 index 00000000..ca5b1ad4 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpset.x @@ -0,0 +1,181 @@ +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +# DP_SETS -- Set a daophot string parameter. + +procedure dp_sets (dao, param, str) + +pointer dao # pointer to daophot structure +int param # parameter +char str[ARB] # string value + +begin + switch (param) { + case INIMAGE: + call strcpy (str, DP_INIMAGE(dao), SZ_FNAME) + case INPHOTFILE: + call strcpy (str, DP_INPHOTFILE(dao), SZ_FNAME) + case COORDS: + call strcpy (str, DP_COORDS(dao), SZ_FNAME) + case PSFIMAGE: + call strcpy (str, DP_PSFIMAGE(dao), SZ_FNAME) + case OUTPHOTFILE: + call strcpy (str, DP_OUTPHOTFILE(dao), SZ_FNAME) + case OUTREJFILE: + call strcpy (str, DP_OUTREJFILE(dao), SZ_FNAME) + case OUTIMAGE: + call strcpy (str, DP_OUTIMAGE(dao), SZ_FNAME) + case IFILTER: + call strcpy (str, DP_IFILTER(dao), SZ_FNAME) + case OTIME: + call strcpy (str, DP_OTIME(dao), SZ_FNAME) + case CCDGAIN: + call strcpy (str, DP_CCDGAIN(dao), SZ_FNAME) + case CCDREAD: + call strcpy (str, DP_CCDREAD(dao), SZ_FNAME) + case EXPTIME: + call strcpy (str, DP_EXPTIME(dao), SZ_FNAME) + case OBSTIME: + call strcpy (str, DP_OBSTIME(dao), SZ_FNAME) + case AIRMASS: + call strcpy (str, DP_AIRMASS(dao), SZ_FNAME) + case FILTER: + call strcpy (str, DP_FILTER(dao), SZ_FNAME) + case FUNCTION: + call strcpy (str, DP_FUNCTION(dao), SZ_FNAME) + case FUNCLIST: + call strcpy (str, DP_FUNCLIST(dao), SZ_FNAME) + default: + call error (0, "DP_SETS: Unknown daophot string parameter") + } +end + + +# DP_SETI -- Set a daophot integer parameter. + +procedure dp_seti (dao, param, ival) + +pointer dao # pointer to daophot structure +int param # parameter +int ival # integer value + +pointer apsel + +begin + apsel = DP_APSEL(dao) + + switch (param) { + case MW: + DP_MW(dao) = ival + case WCSIN: + DP_WCSIN(dao) = ival + case WCSOUT: + DP_WCSOUT(dao) = ival + case WCSPSF: + DP_WCSPSF(dao) = ival + case CTIN: + DP_CTIN(dao) = ival + case CTOUT: + DP_CTOUT(dao) = ival + case CTPSF: + DP_CTPSF(dao) = ival + case MAXITER: + DP_MAXITER(dao) = ival + case VERBOSE: + DP_VERBOSE(dao) = ival + case TEXT: + DP_TEXT(dao) = ival + case MAXNSTAR: + DP_MAXNSTAR(dao) = ival + case MAXGROUP: + DP_MAXGROUP(dao) = ival + case CLIPEXP: + DP_CLIPEXP(dao) = ival + case RECENTER: + DP_RECENTER(dao) = ival + case FITSKY: + DP_FITSKY(dao) = ival + case GROUPSKY: + DP_GROUPSKY(dao) = ival + case VARORDER: + DP_VARORDER(dao) = ival + case FEXPAND: + DP_FEXPAND(dao) = ival + case SATURATED: + DP_SATURATED(dao) = ival + case NCLEAN: + DP_NCLEAN(dao) = ival + case APNUM: + DP_APNUM(apsel) = ival + default: + call error (0, "DP_SETI: Unknown integer daophot parameter") + } +end + + +# DP_SETR -- Set a real daophot parameter. + +procedure dp_setr (dao, param, rval) + +pointer dao # pointer to daophot structure +int param # parameter +real rval # real value + +begin + switch (param) { + case SCALE: + DP_SCALE(dao) = rval + case SFWHMPSF: + DP_SFWHMPSF(dao) = rval + case FWHMPSF: + DP_FWHMPSF(dao) = rval + case MAXGDATA: + DP_MAXGDATA(dao) = rval + case MINGDATA: + DP_MINGDATA(dao) = rval + case READNOISE: + DP_READNOISE(dao) = rval + case PHOTADU: + DP_PHOTADU(dao) = rval + case RPSFRAD: + DP_RPSFRAD(dao) = rval + case SPSFRAD: + DP_SPSFRAD(dao) = rval + case PSFRAD: + DP_PSFRAD(dao) = rval + case SFITRAD: + DP_SFITRAD(dao) = rval + case FITRAD: + DP_FITRAD(dao) = rval + case SMATCHRAD: + DP_SMATCHRAD(dao) = rval + case MATCHRAD: + DP_MATCHRAD(dao) = rval + case SANNULUS: + DP_SANNULUS(dao) = rval + case ANNULUS: + DP_ANNULUS(dao) = rval + case SDANNULUS: + DP_SDANNULUS(dao) = rval + case DANNULUS: + DP_DANNULUS(dao) = rval + case CRITSNRATIO: + DP_CRITSNRATIO(dao) = rval + case CLIPRANGE: + DP_CLIPRANGE(dao) = rval + case XAIRMASS: + DP_XAIRMASS(dao) = rval + case ITIME: + DP_ITIME(dao) = rval + case FLATERR: + DP_FLATERR(dao) = rval + case PROFERR: + DP_PROFERR(dao) = rval + case SMERGERAD: + DP_SMERGERAD(dao) = rval + case MERGERAD: + DP_MERGERAD(dao) = rval + default: + call error (0, "DP_SETR: Unknown real daophot parameter") + } +end diff --git a/noao/digiphot/daophot/daolib/dpstat.x b/noao/digiphot/daophot/daolib/dpstat.x new file mode 100644 index 00000000..45bcf82e --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpstat.x @@ -0,0 +1,180 @@ +include "../lib/daophotdef.h" +include "../lib/apseldef.h" + +# DP_STATS -- Fetch a daophot string parameter. + +procedure dp_stats (dao, param, str, maxch) + +pointer dao # pointer to daophot structure +int param # parameter +char str[ARB] # string value +int maxch # maximum number of characters + +begin + switch (param) { + case INIMAGE: + call strcpy (DP_INIMAGE(dao), str, maxch) + case INPHOTFILE: + call strcpy (DP_INPHOTFILE(dao), str, maxch) + case COORDS: + call strcpy (DP_COORDS(dao), str, maxch) + case PSFIMAGE: + call strcpy (DP_PSFIMAGE(dao), str, maxch) + case OUTPHOTFILE: + call strcpy (DP_OUTPHOTFILE(dao), str, maxch) + case OUTIMAGE: + call strcpy (DP_OUTIMAGE(dao), str, maxch) + case OUTREJFILE: + call strcpy (DP_OUTREJFILE(dao), str, maxch) + case IFILTER: + call strcpy (DP_IFILTER(dao), str, maxch) + case OTIME: + call strcpy (DP_OTIME(dao), str, maxch) + case CCDGAIN: + call strcpy (DP_CCDGAIN(dao), str, maxch) + case CCDREAD: + call strcpy (DP_CCDREAD(dao), str, maxch) + case EXPTIME: + call strcpy (DP_EXPTIME(dao), str, maxch) + case OBSTIME: + call strcpy (DP_OBSTIME(dao), str, maxch) + case FILTER: + call strcpy (DP_FILTER(dao), str, maxch) + case AIRMASS: + call strcpy (DP_AIRMASS(dao), str, maxch) + case FUNCTION: + call strcpy (DP_FUNCTION(dao), str, maxch) + case FUNCLIST: + call strcpy (DP_FUNCLIST(dao), str, maxch) + default: + call error (0, "DP_STATS: Unknown daophot string parameter") + } +end + + +# DP_STATI -- Fetch a daophot integer parameter. + +int procedure dp_stati (dao, param) + +pointer dao # pointer to daophot structure +int param # parameter + +pointer apsel + +begin + apsel = DP_APSEL(dao) + + switch (param) { + case MW: + return (DP_MW(dao)) + case WCSIN: + return (DP_WCSIN(dao)) + case WCSOUT: + return (DP_WCSOUT(dao)) + case WCSPSF: + return (DP_WCSPSF(dao)) + case CTIN: + return (DP_CTIN(dao)) + case CTOUT: + return (DP_CTOUT(dao)) + case CTPSF: + return (DP_CTPSF(dao)) + case MAXITER: + return (DP_MAXITER(dao)) + case VERBOSE: + return (DP_VERBOSE(dao)) + case TEXT: + return (DP_TEXT(dao)) + case MAXNSTAR: + return (DP_MAXNSTAR(dao)) + case MAXGROUP: + return (DP_MAXGROUP(dao)) + case CLIPEXP: + return (DP_CLIPEXP(dao)) + case RECENTER: + return (DP_RECENTER(dao)) + case FITSKY: + return (DP_FITSKY(dao)) + case GROUPSKY: + return (DP_GROUPSKY(dao)) + case VARORDER: + return (DP_VARORDER(dao)) + case FEXPAND: + return (DP_FEXPAND(dao)) + case SATURATED: + return (DP_SATURATED(dao)) + case NCLEAN: + return (DP_NCLEAN(dao)) + case APNUM: + return (DP_APNUM(apsel)) + default: + call error (0, "DP_STATI: Unknown integer daophot parameter") + } +end + + +# DP_STATR -- Fetch a daophot real parameter. + +real procedure dp_statr (dao, param) + +pointer dao # pointer to daophot structure +int param # parameter + +begin + switch (param) { + case SCALE: + return (DP_SCALE(dao)) + case FWHMPSF: + return (DP_FWHMPSF(dao)) + case SFWHMPSF: + return (DP_SFWHMPSF(dao)) + case MAXGDATA: + return (DP_MAXGDATA(dao)) + case MINGDATA: + return (DP_MINGDATA(dao)) + case READNOISE: + return (DP_READNOISE(dao)) + case PHOTADU: + return (DP_PHOTADU(dao)) + case RPSFRAD: + return (DP_RPSFRAD(dao)) + case SPSFRAD: + return (DP_SPSFRAD(dao)) + case PSFRAD: + return (DP_PSFRAD(dao)) + case SFITRAD: + return (DP_SFITRAD(dao)) + case FITRAD: + return (DP_FITRAD(dao)) + case SMATCHRAD: + return (DP_SMATCHRAD(dao)) + case MATCHRAD: + return (DP_MATCHRAD(dao)) + case SANNULUS: + return (DP_SANNULUS(dao)) + case ANNULUS: + return (DP_ANNULUS(dao)) + case SDANNULUS: + return (DP_SDANNULUS(dao)) + case DANNULUS: + return (DP_DANNULUS(dao)) + case CRITSNRATIO: + return (DP_CRITSNRATIO(dao)) + case CLIPRANGE: + return (DP_CLIPRANGE(dao)) + case XAIRMASS: + return (DP_XAIRMASS(dao)) + case ITIME: + return (DP_ITIME(dao)) + case FLATERR: + return (DP_FLATERR(dao)) + case PROFERR: + return (DP_PROFERR(dao)) + case SMERGERAD: + return (DP_SMERGERAD(dao)) + case MERGERAD: + return (DP_MERGERAD(dao)) + default: + call error (0, "DP_STATR: Unknown real daophot parameter") + } +end diff --git a/noao/digiphot/daophot/daolib/dpverify.x b/noao/digiphot/daophot/daolib/dpverify.x new file mode 100644 index 00000000..dcc721ad --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpverify.x @@ -0,0 +1,563 @@ +include "../lib/daophotdef.h" + +# DP_VFUNCTION -- Verify the analytic psf function. + +procedure dp_vfunction (dao) + +pointer dao # pointer to the daophot structure. + +int len +pointer sp, str, pstr +int scan(), nscan(), strlen(), dp_fctdecode(), dp_strwrd() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (pstr, SZ_FNAME, TY_CHAR) + + # Print the current function list. + len = strlen (DP_FUNCLIST(dao)) + call strcpy (DP_FUNCLIST(dao), Memc[pstr], len - 1) + call printf ("Analytic psf function(s) (%s): ") + call pargstr (Memc[pstr+1]) + call flush (STDOUT) + + # Confirm the PSF function type. + if (scan() == EOF) + ; + else { + call gargstr (Memc[str], SZ_FNAME) + if (nscan () != 1) + ; + else if (dp_fctdecode (Memc[str], Memc[pstr], SZ_FNAME) <= 0) + ; + else + call strcpy (Memc[pstr], DP_FUNCLIST(dao), SZ_FNAME) + } + + # Print the confirmed function list. + len = strlen (DP_FUNCLIST(dao)) + call strcpy (DP_FUNCLIST(dao), Memc[pstr], len - 1) + call printf ( "\tAnalytic psf function(s): %d\n") + call pargstr (Memc[pstr+1]) + + # Set the function type. + if (dp_strwrd (1, Memc[str], SZ_FNAME, DP_FUNCLIST(dao)) <= 0) + call strcpy ("gauss", DP_FUNCTION(dao), SZ_FNAME) + else + call strcpy (Memc[str], DP_FUNCTION(dao), SZ_FNAME) + + call sfree (sp) +end + + +# DP_VVARORDER -- Verify the order of variability of the psf function. + +procedure dp_vvarorder (dao) + +pointer dao # pointer to the daophot structure. + +int varorder +int scan(), nscan() + +begin + # Confirm that the psf is variable. + call printf ("Order of variable psf (%d): ") + call pargi (DP_VARORDER(dao)) + call flush (STDOUT) + if (scan() == EOF) + varorder = DP_VARORDER(dao) + else { + call gargi (varorder) + if (nscan () != 1) + varorder = DP_VARORDER(dao) + else if (varorder < -1 || varorder > 2) + varorder = DP_VARORDER(dao) + } + DP_VARORDER(dao) = varorder + call printf ( "\tOrder of variable psf: %d\n") + call pargi (varorder) +end + + +# DP_VFEXPAND -- Verify whether or not to expand the analytics function. + +procedure dp_vfexpand (dao) + +pointer dao # pointer to the daophot structure. + +bool fexpand +bool itob() +int scan(), nscan(), btoi() + +begin + # Confirm whether of not to expand the analytic function. + call printf ("Expand the analytic function (%b): ") + call pargb (itob (DP_FEXPAND(dao))) + call flush (STDOUT) + if (scan() == EOF) + fexpand = itob (DP_FEXPAND(dao)) + else { + call gargb (fexpand) + if (nscan () != 1) + fexpand = itob (DP_FEXPAND(dao)) + } + DP_FEXPAND(dao) = btoi (fexpand) + call printf ( "\tExpand analytic fucntion: %b\n") + call pargb (fexpand) +end + + +# DP_VSATURATED -- Verify whether or not to use saturated stars in the +# psf computation. + +procedure dp_vsaturated (dao) + +pointer dao # pointer to the daophot structure. + +bool saturated +bool itob() +int scan(), nscan(), btoi() + +begin + # Confirm whether of not to use saturated psf stars. + call printf ("Use saturated psf stars (%b): ") + call pargb (itob (DP_SATURATED(dao))) + call flush (STDOUT) + if (scan() == EOF) + saturated = itob (DP_SATURATED(dao)) + else { + call gargb (saturated) + if (nscan () != 1) + saturated = itob (DP_SATURATED(dao)) + } + DP_SATURATED(dao) = btoi (saturated) + call printf ( "\tUse saturated psf stars: %b\n") + call pargb (saturated) +end + + +# DP_VFWHMPSF -- Confirm the fwhm of the psf. + +procedure dp_vfwhmpsf (dao) + +pointer dao # pointer to the daophot structure + +real sfwhmpsf, fwhmpsf +int scan(), nscan() + +begin + # Confirm the psf radius. + call printf ( "Fwhm of psf in scale units (%g): ") + call pargr (DP_SFWHMPSF(dao)) + call flush (STDOUT) + if (scan() == EOF) + sfwhmpsf = DP_SFWHMPSF(dao) + else { + call gargr (sfwhmpsf) + if (nscan () != 1) + sfwhmpsf = DP_SFWHMPSF(dao) + } + fwhmpsf = sfwhmpsf / DP_SCALE(dao) + + DP_SFWHMPSF(dao) = sfwhmpsf + DP_FWHMPSF(dao) = fwhmpsf + + call printf ( "\tNew fwhm: %g scale units %g pixels\n") + call pargr (sfwhmpsf) + call pargr (fwhmpsf) +end + + +# DP_VPSFRAD -- Confirm the psf radius. + +procedure dp_vpsfrad (dao) + +pointer dao # pointer to the daophot structure + +real rpsfrad, psfrad +int scan(), nscan() + +begin + # Confirm the psf radius. + call printf ( "Psf radius in scale units (%g): ") + call pargr (DP_RPSFRAD(dao)) + call flush (STDOUT) + if (scan() == EOF) + rpsfrad = DP_RPSFRAD(dao) + else { + call gargr (rpsfrad) + if (nscan () != 1) + rpsfrad = DP_RPSFRAD(dao) + } + psfrad = rpsfrad / DP_SCALE(dao) + + DP_RPSFRAD(dao) = rpsfrad + DP_SPSFRAD(dao) = rpsfrad + DP_PSFRAD(dao) = psfrad + + call printf ( "\tNew psf radius: %g scale units %g pixels\n") + call pargr (rpsfrad) + call pargr (psfrad) +end + + +# DP_VFITRAD -- Confirm the fitting radius. + +procedure dp_vfitrad (dao) + +pointer dao # pointer to the daophot structures + +real sfitrad, fitrad +int scan(), nscan() + +begin + # Confirm the fitting radius. + call printf ( "Fitting radius in scale units (%g): ") + call pargr (DP_SFITRAD(dao)) + call flush (STDOUT) + if (scan() == EOF) + sfitrad = DP_SFITRAD(dao) + else { + call gargr (sfitrad) + if (nscan () != 1) + sfitrad = DP_SFITRAD(dao) + } + fitrad = sfitrad / DP_SCALE(dao) + + DP_SFITRAD(dao) = sfitrad + DP_FITRAD(dao) = fitrad + + call printf ( "\tNew fitting radius: %g scale units %g pixels\n") + call pargr (sfitrad) + call pargr (fitrad) +end + + +# DP_VMATCHRAD -- Confirm the matching radius. + +procedure dp_vmatchrad (dao) + +pointer dao # pointer to the daophot structure + +real smatchrad, matchrad +int scan(), nscan() + +begin + # Confirm the matching radius. + call printf ( "Matching radius in scale units (%g): ") + call pargr (DP_SMATCHRAD(dao)) + call flush (STDOUT) + if (scan() == EOF) + smatchrad = DP_SMATCHRAD(dao) + else { + call gargr (smatchrad) + if (nscan () != 1) + smatchrad = DP_SMATCHRAD(dao) + } + matchrad = smatchrad / DP_SCALE(dao) + + DP_SMATCHRAD(dao) = smatchrad + DP_MATCHRAD(dao) = matchrad + + call printf ( "\tNew matching radius: %g scale units %g pixels\n") + call pargr (smatchrad) + call pargr (matchrad) +end + + +# DP_VMERGERAD -- Confirm the merging radius. + +procedure dp_vmergerad (dao) + +pointer dao # pointer to the daophot structure + +real smergerad, mergerad +int scan(), nscan() + +begin + # Confirm the merging radius. + call printf ( "Merging radius in scale units (%g): ") + call pargr (DP_SMERGERAD(dao)) + call flush (STDOUT) + if (scan() == EOF) + smergerad = DP_SMERGERAD(dao) + else { + call gargr (smergerad) + if (nscan () != 1) + smergerad = DP_SMERGERAD(dao) + } + if (IS_INDEFR(smergerad)) + mergerad = INDEFR + else + mergerad = smergerad / DP_SCALE(dao) + + DP_SMERGERAD(dao) = smergerad + DP_MERGERAD(dao) = mergerad + + call printf ( "\tNew merging radius: %g scale units %g pixels\n") + call pargr (smergerad) + call pargr (mergerad) +end + + +# DP_VFITRAD -- Confirm the fitting radius. + + +# DP_VDATAMIN-- Verify the minimum good data value. + +procedure dp_vdatamin (dao) + +pointer dao # pointer to the daophot structure + +real datamin +int scan(), nscan() + +begin + # Confirm the threshold parameter. + call printf ("Minimum good data value (%g) (CR or value): ") + call pargr (DP_MINGDATA(dao)) + call flush (STDOUT) + if (scan() == EOF) + datamin = DP_MINGDATA(dao) + else { + call gargr (datamin) + if (nscan () != 1) + datamin = DP_MINGDATA(dao) + } + DP_MINGDATA(dao) = datamin + + call printf ("\tNew minimum good data value: %g counts\n") + call pargr (datamin) +end + + +# DP_VDATAMAX-- Verify the maximum good data value. + +procedure dp_vdatamax (dao) + +pointer dao # pointer to the daophot structure + +real datamax +int scan(), nscan() + +begin + # Confirm the threshold parameter. + call printf ("Maximum good data value (%g) (CR or value): ") + call pargr (DP_MAXGDATA(dao)) + call flush (STDOUT) + if (scan() == EOF) + datamax = DP_MAXGDATA(dao) + else { + call gargr (datamax) + if (nscan () != 1) + datamax = DP_MAXGDATA(dao) + } + DP_MAXGDATA(dao) = datamax + + call printf ("\tNew maximum good data value: %g counts\n") + call pargr (datamax) +end + + +# DP_VMAXGROUP -- Verify the maximum group size. + +procedure dp_vmaxgroup (dao) + +pointer dao # pointer to the daophot strucuture + +int maxgroup +int scan(), nscan() + +begin + call printf ( "Maximum group size in number of stars (%d): ") + call pargi (DP_MAXGROUP(dao)) + call flush (STDOUT) + if (scan() == EOF) + maxgroup = DP_MAXGROUP(dao) + else { + call gargi (maxgroup) + if (nscan () != 1) + maxgroup = DP_MAXGROUP(dao) + } + DP_MAXGROUP(dao) = maxgroup + + call printf ( "\tNew maximum group size: %d stars\n") + call pargi (maxgroup) +end + + +# DP_VCRITSNRATIO -- Verify the critical signal-to-noise ratio. + +procedure dp_vcritnsratio (dao) + +pointer dao # pointer to the daophot structure + +real critsnratio +int scan(), nscan() + +begin + # Confirm the critical signal-to-noise ratio. + call printf ( "Critical S/N ratio in stdevs per pixel (%g): ") + call pargr (DP_CRITSNRATIO(dao)) + call flush (STDOUT) + if (scan() == EOF) + critsnratio = DP_CRITSNRATIO(dao) + else { + call gargr (critsnratio) + if (nscan () != 1) + critsnratio = DP_CRITSNRATIO(dao) + } + DP_CRITSNRATIO(dao) = critsnratio + + call printf ( "\tNew critical S/N ratio: %g stdevs per pixel\n") + call pargr (critsnratio) +end + + +# DP_VRECENTER -- Verify whether or not to recenter the stars. + +procedure dp_vrecenter (dao) + +pointer dao # pointer to the daophot structure. + +bool recenter +bool itob() +int scan(), nscan(), btoi() + +begin + # Confirm whether of not to recenter the stars. + call printf ("Recenter the stars (%b): ") + call pargb (itob (DP_RECENTER(dao))) + call flush (STDOUT) + if (scan() == EOF) + recenter = itob (DP_RECENTER(dao)) + else { + call gargb (recenter) + if (nscan () != 1) + recenter = itob (DP_RECENTER(dao)) + } + DP_RECENTER(dao) = btoi (recenter) + call printf ( "\tRecenter the stars: %b\n") + call pargb (recenter) +end + + +# DP_VFITSKY -- Verify whether or not to refit the sky value. + +procedure dp_vfitsky (dao) + +pointer dao # pointer to the daophot structure. + +bool fitsky +bool itob() +int scan(), nscan(), btoi() + +begin + # Confirm whether of not to refit the sky. + call printf ("Refit the sky (%b): ") + call pargb (itob (DP_FITSKY(dao))) + call flush (STDOUT) + if (scan() == EOF) + fitsky = itob (DP_FITSKY(dao)) + else { + call gargb (fitsky) + if (nscan () != 1) + fitsky = itob (DP_FITSKY(dao)) + } + DP_FITSKY(dao) = btoi (fitsky) + call printf ( "\tRefit the sky: %b\n") + call pargb (fitsky) +end + + +# DP_VGROUPSKY -- Verify whether or not to fit group sky values. + +procedure dp_vgroupsky (dao) + +pointer dao # pointer to the daophot structure + +bool groupsky +bool itob() +int scan(), nscan(), btoi() + +begin + # Confirm whether of not to use group sky values. + call printf ("Use group sky values (%b): ") + call pargb (itob (DP_GROUPSKY(dao))) + call flush (STDOUT) + if (scan() == EOF) + groupsky = itob (DP_GROUPSKY(dao)) + else { + call gargb (groupsky) + if (nscan () != 1) + groupsky = itob (DP_GROUPSKY(dao)) + } + DP_GROUPSKY(dao) = btoi (groupsky) + call printf ( "\tUse group sky values: %b\n") + call pargb (groupsky) +end + + +# DP_VSANNULUS -- Confirm the inner radius of the sky fitting annulus. + +procedure dp_vsannulus (dao) + +pointer dao # pointer to the daophot structure + +real sannulus, annulus +int scan(), nscan() + +begin + # Confirm the inner radius of the sky annulus. + call printf ( "Inner radius of sky annulus in scale units (%g): ") + call pargr (DP_SANNULUS(dao)) + call flush (STDOUT) + if (scan() == EOF) + sannulus = DP_SANNULUS(dao) + else { + call gargr (sannulus) + if (nscan () != 1) + sannulus = DP_SANNULUS(dao) + } + annulus = sannulus / DP_SCALE(dao) + + DP_SANNULUS(dao) = sannulus + DP_ANNULUS(dao) = annulus + + call printf ( "\tNew inner radius: %g scale units %g pixels\n") + call pargr (sannulus) + call pargr (annulus) +end + + +# DP_VWSANNULUS -- Confirm the width of the sky fitting annulus. + +procedure dp_vwsannulus (dao) + +pointer dao # pointer to the daophot structure + +real sdannulus, dannulus +int scan(), nscan() + +begin + # Confirm the inner radius of the sky annulus. + call printf ( "Width of sky annulus in scale units (%g): ") + call pargr (DP_SDANNULUS(dao)) + call flush (STDOUT) + if (scan() == EOF) + sdannulus = DP_SDANNULUS(dao) + else { + call gargr (sdannulus) + if (nscan () != 1) + sdannulus = DP_SDANNULUS(dao) + } + dannulus = sdannulus / DP_SCALE(dao) + + DP_SDANNULUS(dao) = sdannulus + DP_DANNULUS(dao) = dannulus + + call printf ( "\tNew annulus width: %g scale units %g pixels\n") + call pargr (sdannulus) + call pargr (dannulus) +end diff --git a/noao/digiphot/daophot/daolib/dpwcs.x b/noao/digiphot/daophot/daolib/dpwcs.x new file mode 100644 index 00000000..9961b933 --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpwcs.x @@ -0,0 +1,234 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imio.h> +include "../lib/daophotdef.h" + +# DP_WIN -- Convert the input coordinates to logical coordinates. + +procedure dp_win (dp, im, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +pointer im # the input image descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +int dp_stati() + +begin + # Transform the input coordinates. + switch (dp_stati (dp, WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call dp_itol (dp, xin, yin, xout, yout, npts) + case WCS_TV: + call dp_vtol (im, xin, yin, xout, yout, npts) + default: + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + } +end + + +# DP_WOUT -- Convert the logical coordinates to output coordinates. + +procedure dp_wout (dp, im, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +pointer im # the input image descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +int dp_stati() + +begin + # Transform the input coordinates. + switch (dp_stati (dp, WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call dp_ltoo (dp, xin, yin, xout, yout, npts) + case WCS_TV: + call dp_ltov (im, xin, yin, xout, yout, npts) + default: + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + } +end + + +# DP_WPSF -- Convert the logical coordinates to psf coordinates. + +procedure dp_wpsf (dp, im, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +pointer im # the input image descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +int dp_stati() + +begin + # Transform the input coordinates. + switch (dp_stati (dp, WCSPSF)) { + case WCS_WORLD, WCS_PHYSICAL: + call dp_ltop (dp, xin, yin, xout, yout, npts) + case WCS_TV: + call dp_ltov (im, xin, yin, xout, yout, npts) + default: + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + } +end + + +# DP_ITOL -- Convert coordinates from the input coordinate system to the +# logical coordinate system. + +procedure dp_itol (dp, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +double xt, yt +pointer ct +int i +int dp_stati() + +begin + ct = dp_stati (dp, CTIN) + if (ct == NULL) { + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + return + } + + do i = 1, npts { + call mw_c2trand (ct, double (xin[i]), double (yin[i]), xt, yt) + xout[i] = xt + yout[i] = yt + } +end + + +# DP_LTOO -- Convert coordinates from the logical coordinate system to the +# output coordinate system. + +procedure dp_ltoo (dp, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +double xt, yt +pointer ct +int i +int dp_stati() + +begin + ct = dp_stati (dp, CTOUT) + if (ct == NULL) { + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + return + } + + do i = 1, npts { + call mw_c2trand (ct, double (xin[i]), double (yin[i]), xt, yt) + xout[i] = xt + yout[i] = yt + } +end + + +# DP_LTOP -- Convert coordinates from the logical coordinate system to the +# output coordinate system. + +procedure dp_ltop (dp, xin, yin, xout, yout, npts) + +pointer dp # the apphot package descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +double xt, yt +pointer ct +int i +int dp_stati() + +begin + ct = dp_stati (dp, CTPSF) + if (ct == NULL) { + call amovr (xin, xout, npts) + call amovr (yin, yout, npts) + return + } + + do i = 1, npts { + call mw_c2trand (ct, double (xin[i]), double (yin[i]), xt, yt) + xout[i] = xt + yout[i] = yt + } +end + + +# DP_LTOV -- Convert coordinate from the logical coordinate system to the +# output coordinate system. + +procedure dp_ltov (im, xin, yin, xout, yout, npts) + +pointer im # the input image descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +int i, index1, index2 + +begin + index1 = IM_VMAP(im,1) + index2 = IM_VMAP(im,2) + do i = 1, npts { + xout[i] = xin[i] * IM_VSTEP(im,index1) + IM_VOFF(im,index1) + yout[i] = yin[i] * IM_VSTEP(im,index2) + IM_VOFF(im,index2) + } +end + + +# DP_VTOL -- Convert coordinate from the tv coordinate system to the +# logical coordinate system. + +procedure dp_vtol (im, xin, yin, xout, yout, npts) + +pointer im # the input image descriptor +real xin[ARB] # the input x coordinate +real yin[ARB] # the input y coordinate +real xout[ARB] # the output x coordinate +real yout[ARB] # the output y coordinate +int npts # the number of coordinates to convert + +int i, index1, index2 + +begin + index1 = IM_VMAP(im,1) + index2 = IM_VMAP(im,2) + do i = 1, npts { + xout[i] = (xin[i] - IM_VOFF(im,index1)) / IM_VSTEP(im,index1) + yout[i] = (yin[i] - IM_VOFF(im,index2)) / IM_VSTEP(im,index2) + } +end diff --git a/noao/digiphot/daophot/daolib/dpwparam.x b/noao/digiphot/daophot/daolib/dpwparam.x new file mode 100644 index 00000000..6e122dee --- /dev/null +++ b/noao/digiphot/daophot/daolib/dpwparam.x @@ -0,0 +1,98 @@ +# DP_RPARAM -- Encode a daophot real parameter. + +procedure dp_rparam (out, keyword, value, units, comments) + +int out # output file descriptor +char keyword[ARB] # keyword string +real value # parameter value +char units[ARB] # units string +char comments[ARB] # comment string + +begin + if (out == NULL) + return + + call strupr (keyword) + call fprintf (out, + "#K%4t%-10.10s%14t = %17t%-23.7g%41t%-10.10s%52t%-10s\n") + call pargstr (keyword) + call pargr (value) + call pargstr (units) + call pargstr ("%-23.7g") + call pargstr (comments) +end + + +# DP_IPARAM -- Encode a daophot integer parameter. + +procedure dp_iparam (out, keyword, value, units, comments) + +int out # output file descriptor +char keyword[ARB] # keyword string +int value # parameter value +char units[ARB] # units string +char comments[ARB] # comment string + +begin + if (out == NULL) + return + + call strupr (keyword) + call fprintf (out, + "#K%4t%-10.10s%14t = %17t%-23d%41t%-10.10s%52t%-10s\n") + call pargstr (keyword) + call pargi (value) + call pargstr (units) + call pargstr ("%-23d") + call pargstr (comments) +end + + +# DP_BPARAM -- Encode a daophot boolean parameter. + +procedure dp_bparam (out, keyword, value, units, comments) + +int out # output file descriptor +char keyword[ARB] # keyword string +bool value # parameter value +char units[ARB] # units string +char comments[ARB] # comment string + +begin + if (out == NULL) + return + + call strupr (keyword) + call fprintf (out, + "#K%4t%-10.10s%14t = %17t%-23b%41t%-10.10s%52t%-10s\n") + call pargstr (keyword) + call pargb (value) + call pargstr (units) + call pargstr ("%-23b") + call pargstr (comments) +end + + +# DP_SPARAM -- Encode a daophot string parameter. + +procedure dp_sparam (out, keyword, value, units, comments) + +int out # output file descriptor +char keyword[ARB] # keyword string +char value[ARB] # parameter value +char units[ARB] # units string +char comments[ARB] # comment string + +begin + if (out == NULL) + return + + call strupr (keyword) + call fprintf (out, + "#K%4t%-10.10s%14t = %17t%-23.23s%41t%-10.10s%52t%-10s\n") + call pargstr (keyword) + call pargstr (value) + call pargstr (units) + call pargstr ("%-23s") + call pargstr (comments) +end diff --git a/noao/digiphot/daophot/daolib/erf.x b/noao/digiphot/daophot/daolib/erf.x new file mode 100644 index 00000000..d8a48f50 --- /dev/null +++ b/noao/digiphot/daophot/daolib/erf.x @@ -0,0 +1,81 @@ +define NGL 4 + +# DAOERF -- Numerically integrate a Gaussian function from xin-0.5 to xin+0.5 +# using 4 point Gauss-Legendre integration. Beta is the half-width at +# half-maximum which is equal to 1.17741 * sigma. The Gaussian function is +# shown below. +# +# erf = exp (-0.5 *((x - x0) / beta) ** 2)) +# +# or +# +# erf = exp (-0.6931472 * [(x - xo) / beta] ** 2) +# +# Also provide the first derivative of the integral with respect to xo and beta. + +real procedure daoerf (xin, x0, beta, dfdx0, dfdbet) + +real xin # the input value +real x0 # position of Gaussian peak +real beta # sigma of the Gaussian +real dfdx0 # derivative of Gaussian wrt x0 +real dfdbet # derivative of Gaussian wrt beta + +int i, npt +real betasq, deltax, erfval, xsq, f, x, wf +real dx[NGL,NGL], wt[NGL,NGL] +data dx / 0.0, 0.0, 0.0, 0.0, + -0.28867513, 0.28867513, 0.0, 0.0, + -0.38729833, 0.0, 0.38729833, 0.0, + -0.43056816, -0.16999052, 0.16999052, 0.43056816 / +data wt / 1.0, 0.0, 0.0, 0.0, + 0.5, 0.5, 0.0, 0.0, + 0.27777778, 0.44444444, 0.27777778, 0.0, + 0.17392742, 0.32607258, 0.32607258, 0.17392742 / + +begin + # Compute some constants. + betasq = beta ** 2 + deltax = xin - x0 + + # Initialize. + erfval = 0.0 + dfdx0 = 0.0 + dfdbet = 0.0 + + xsq = deltax ** 2 + f = xsq / betasq + if (f > 34.0) + return (erfval) + + f = exp (-0.6931472 * f) + if (f >= 0.046) { + npt = 4 + } else if (f >= 0.0022) { + npt = 3 + } else if (f >= 0.0001) { + npt = 2 + } else if (f >= 1.0e-10) { + erfval = f + dfdx0 = 1.3862944 * deltax * f / betasq + dfdbet = 1.3862944 * xsq * f / (betasq * beta) + return (erfval) + } else { + return (erfval) + } + + do i = 1, npt { + x = deltax + dx[i,npt] + xsq = x ** 2 + f = exp (-0.6931472 * xsq / betasq) + wf = wt[i,npt] * f + erfval = erfval + wf + dfdx0 = dfdx0 + x * wf + dfdbet = dfdbet + xsq * wf + } + + dfdx0 = 1.3862944 * dfdx0 / betasq + dfdbet = 1.3862944 * dfdbet / (betasq * beta) + + return (erfval) +end diff --git a/noao/digiphot/daophot/daolib/invers.f b/noao/digiphot/daophot/daolib/invers.f new file mode 100644 index 00000000..97d9b582 --- /dev/null +++ b/noao/digiphot/daophot/daolib/invers.f @@ -0,0 +1,112 @@ +c subroutine invers (a, max, n, iflag) +c +c Although it seems counter-intuitive, the tests that I have run +c so far suggest that the 180 x 180 matrices that NSTAR needs can +c be inverted with sufficient accuracy if the elements are REAL*4 +c rather than REAL*8. +c +c Arguments +c +c a (input/output) is a square matrix of dimension N. The inverse +c of the input matrix A is returned in A. +c +c max (input) is the size assigned to the matrix A in the calling +c routine. It is needed for the dimension statement below. +c +c iflag (output) is an error flag. iflag = 1 if the matrix could not +c be inverted; iflag = 0 if it could. +c + subroutine invers (a, max, n, iflag) +c + implicit none + integer max, n, iflag + real a(max,max) + integer i, j, k +c + iflag = 0 + i = 1 + 300 if (a(i,i) .eq. 0.0e0) go to 9100 + a(i,i) = 1.0e0 / a(i,i) + j = 1 + 301 if (j .eq. i) go to 304 + a(j,i) = -a(j,i) * a(i,i) + k = 1 + 302 if (k .eq. i) go to 303 + a(j,k) = a(j,k) + a(j,i) * a(i,k) + 303 if (k .eq. n) go to 304 + k = k + 1 + go to 302 + 304 if (j .eq. n) go to 305 + j = j + 1 + go to 301 + 305 k = 1 + 306 if (k .eq. i) go to 307 + a(i,k) = a(i,k) * a(i,i) + 307 if (k .eq. n) go to 308 + k = k + 1 + go to 306 + 308 if (i .eq. n) return + i = i+1 + go to 300 +c +c Error: zero on the diagonal. +c + 9100 iflag = 1 + return +c + end +c +c +c +c subroutine dinvers (a, max, n, iflag) +c +c Arguments +c +c a (input/output) is a square matrix of dimension N. The inverse +c of the input matrix A is returned in A. +c +c max (input) is the size assigned to the matrix A in the calling +c routine. It's needed for the dimension statement below. +c +c iflag (output) is an error flag. iflag = 1 if the matrix could not +c be inverted; iflag = 0 if it could. +c + subroutine dinvers (a, max, n, iflag) +c + implicit none + integer max, n, iflag + double precision a(max,max) + integer i, j, k +c + iflag = 0 + i = 1 + 300 if (a(i,i) .eq. 0.0e0) go to 9100 + a(i,i) = 1.0e0 / a(i,i) + j = 1 + 301 if (j .eq. i) go to 304 + a(j,i) = -a(j,i) * a(i,i) + k = 1 + 302 if (k .EQ. i) go to 303 + a(j,k) = a(j,k) + a(j,i) * a(i,k) + 303 if (k .eq. n) go to 304 + k = k + 1 + go to 302 + 304 if (j .eq. n) go to 305 + j = j + 1 + go to 301 + 305 k = 1 + 306 if (k .eq. i) go to 307 + a(i,k) = a(i,k) * a(i,i) + 307 if (k .eq. n) go to 308 + k = k + 1 + go to 306 + 308 if (i .eq. n) return + i = i+1 + go to 300 +c +c Error: zero on the diagonal. +c + 9100 iflag = 1 + return +c + end diff --git a/noao/digiphot/daophot/daolib/invers2.x b/noao/digiphot/daophot/daolib/invers2.x new file mode 100644 index 00000000..5393c6ea --- /dev/null +++ b/noao/digiphot/daophot/daolib/invers2.x @@ -0,0 +1,72 @@ +define TINY (1.0e-15) + +# INVERS +# +# Although it seems counter-intuitive, the tests that I have run +# so far suggest that the 180 x 180 matrices that NSTAR needs can +# be inverted with sufficient accuracy if the elements are REAL*4 +# rather than REAL*8. +# +# Arguments +# +# a (input/output) is a square matrix of dimension N. The inverse +# of the input matrix A is returned in A. +# +# nmax (input) is the size assigned to the matrix A in the calling +# routine. It is needed for the dimension statement below. +# +# iflag (output) is an error flag. iflag = 1 if the matrix could not +# be inverted; iflag = 0 if it could. +# +# This is an SPP translation of the original fortran version with the +# addition of a check for tiny numbers which could cause an FPE. + +procedure invers2 (a, nmax, n, iflag) + +real a[nmax,nmax] +int nmax +int n +int iflag + +int i, j, k + +begin + # Check for tiny numbers. + do i = 1, n + do j = 1, n + if (abs (a[i,j]) < TINY) + a[i,j] = 0e0 + + # Original code. + iflag = 0 + i = 1 + 30 if (a[i,i] .eq. 0.0e0) goto 91 + a[i,i] = 1.0e0 / a[i,i] + j = 1 + 31 if (j .eq. i) goto 34 + a[j,i] = -a[j,i] * a[i,i] + k = 1 + 32 if (k .eq. i) goto 33 + a[j,k] = a[j,k] + a[j,i] * a[i,k] + 33 if (k .eq. n) goto 34 + k = k + 1 + goto 32 + 34 if (j .eq. n) goto 35 + j = j + 1 + goto 31 + 35 k = 1 + 36 if (k .eq. i) goto 37 + a[i,k] = a[i,k] * a[i,i] + 37 if (k .eq. n) goto 38 + k = k + 1 + goto 36 + 38 if (i .eq. n) return + i = i+1 + goto 30 + +# Error: zero on the diagonal. + + 91 iflag = 1 + return + +end diff --git a/noao/digiphot/daophot/daolib/mkpkg b/noao/digiphot/daophot/daolib/mkpkg new file mode 100644 index 00000000..44253943 --- /dev/null +++ b/noao/digiphot/daophot/daolib/mkpkg @@ -0,0 +1,48 @@ +# DAOPHOT Library Tools + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + dpairmass.x <imhdr.h> ../lib/daophotdef.h + dpapheader.x + dpfilter.x <imhdr.h> ../lib/daophotdef.h + dpfree.x ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/allstardef.h + dpgetapert.x ../lib/apseldef.h ../lib/daophotdef.h \ + <tbset.h> ../../lib/ptkeysdef.h + dpgppars.x <ctotok.h> ../lib/daophotdef.h + dpgsvw.x <imio.h> <imhdr.h> \ + <math.h> + dpppars.x ../lib/daophotdef.h + dpimkeys.x ../lib/daophotdef.h + dpinit.x ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/allstardef.h + dpreadpsf.x <imhdr.h> ../lib/daophotdef.h \ + <error.h> + dpset.x ../lib/daophotdef.h ../lib/apseldef.h + dpstat.x ../lib/daophotdef.h ../lib/apseldef.h + dpdate.x <time.h> + dpgsubrast.x <imhdr.h> + dpnames.x + dpotime.x <imhdr.h> ../lib/daophotdef.h + dppadu.x <imhdr.h> ../lib/daophotdef.h + dppcache.x <imhdr.h> <imset.h> + dprdnoise.x <imhdr.h> ../lib/daophotdef.h + dprmwhite.x <ctype.h> + dpverify.x ../lib/daophotdef.h + dpwparam.x + dpwcs.x <imio.h> ../lib/daophotdef.h + bicubic.x + daoran.x + erf.x + invers.f + invers2.x + mvmul.x + quick.f + pctile.f + profile.x ../lib/daophotdef.h + usepsf.x ../lib/daophotdef.h + ; diff --git a/noao/digiphot/daophot/daolib/mvmul.x b/noao/digiphot/daophot/daolib/mvmul.x new file mode 100644 index 00000000..352def99 --- /dev/null +++ b/noao/digiphot/daophot/daolib/mvmul.x @@ -0,0 +1,48 @@ +# MVMUL -- Multply a matrix (left-hand side) by a one dimensional vector +# (right-hand side) and return the resultant vector. + +procedure mvmul (matrix, maxdim, dim, vector, result) + +real matrix [maxdim, maxdim] # input matrix +int maxdim # maximum size of input matrix +int dim # dimension of matrix and vectors +real vector[maxdim] # input vector +real result[maxdim] # iutput vector + +double sum +int i, j + +begin + do i = 1, dim { + sum = 0.0 + do j = 1, dim + sum = sum + double (matrix[j,i]) * double(vector[j]) + result[i] = sum + } + +end + + +# DMVMUL -- Multply a matrix (left-hand side) by a one dimensional vector +# (right-hand side) and return the resultant vector. + +procedure dmvmul (matrix, maxdim, dim, vector, result) + +double matrix [maxdim, maxdim] # input matrix +int maxdim # maximum size of input matrix +int dim # dimension of matrix and vectors +double vector[maxdim] # input vector +double result[maxdim] # iutput vector + +double sum +int i, j + +begin + do i = 1, dim { + sum = 0.0d0 + do j = 1, dim + sum = sum + (matrix[j,i] * vector[j]) + result[i] = sum + } + +end diff --git a/noao/digiphot/daophot/daolib/pctile.f b/noao/digiphot/daophot/daolib/pctile.f new file mode 100644 index 00000000..e8a3f2d7 --- /dev/null +++ b/noao/digiphot/daophot/daolib/pctile.f @@ -0,0 +1,91 @@ +c function pctile (datum, n, npct) +c +c This is a modification of a quick-sorting algorithm, which is intended +c to take in a vector of numbers, and return the value of the npct-th +c element in that vector: +c +c dataum input vector +c n number of elements in dataum +c npct npct-th element +c pctile output value of function +c +c +c The array datum contains randomly ordered data +c +c + real function pctile (datum, n, npct) +c + implicit none + integer n, npct + real datum(1) + integer min0, max0 + real dkey + integer lo, hi, limlo, limhi +c +c Initialize the pointers. +c + npct = max0 (1, min0 (n,npct)) + limlo = 1 + limhi = n +c +c Compare all elements in the sub-vector between limlo and limhi with +c the current key datum. +c + 100 dkey = datum (limlo) + lo = limlo + hi = limhi +c +c If lo equals hi, we have tested all the elements in the current search +c interval. +c + 101 continue + if (lo .eq. hi) go to 200 + if (datum(hi) .le. dkey) go to 109 +c +c The pointer hi is to be left pointing at a datum SMALLER than the +c key, which is intended to be overwritten. +c + hi = hi - 1 +c + goto 101 + 109 datum(lo) = datum(hi) + lo = lo + 1 + 110 continue + if (lo .eq. hi) goto 200 + if (datum(lo) .ge. dkey) go to 119 + lo = lo + 1 +c + goto 110 + 119 datum(hi) = datum(lo) +c +c The pointer LO is to be left pointing at a datum LARGER than the +c key, which is intended to be overwritten. +c + hi = hi - 1 +c + go to 101 +c +c lo and hi are equal, and point at a value which is intended to +c be overwritten. Since all values below this point are less than +c the key and all values above this point are greater than the key, +c this is where we stick the key back into the vector. +c + 200 continue +c +c At this point in the subroutine, all data between limlo and lo-1, +c inclusive, are less than datum (lo), and all data between lo+1 and +c limhi are larger than dataum(lo). If lo = npct, then datum(lo) is +c the value we are looking for. If npct < lo, then we want to sort the +c values of datum from limlo to lo-1, inclusive, whereas if npct > lo, +c then we want to sort the values of datum from lo+1 to limhi, +c inclusive. +c + datum(lo) = dkey + if (npct - lo) 300, 900, 400 + 300 limhi = lo - 1 + go to 100 + 400 limlo = lo + 1 + go to 100 + 900 pctile = datum(lo) + return + end diff --git a/noao/digiphot/daophot/daolib/profile.x b/noao/digiphot/daophot/daolib/profile.x new file mode 100644 index 00000000..96e1a531 --- /dev/null +++ b/noao/digiphot/daophot/daolib/profile.x @@ -0,0 +1,506 @@ +include "../lib/daophotdef.h" + +define NGL 4 + +# DP_PROFILE -- Evaluate the analytic part of the psf function and its +# derivatives. + +real procedure dp_profile (ipstyp, dx, dy, par, dhdxc, dhdyc, term, ideriv) + +int ipstyp # the analytic psf function type +real dx, dy # distance of point from center of function +real par[ARB] # the current parameter values +real dhdxc, dhdyc # derivatives of the function integral wrt x,y +real term[ARB] # derivatives of the function wrt parameters +int ideriv # compute the derivatives ? + +int npt, ix, iy +real d[NGL,NGL], w[NGL,NGL], x[NGL], xsq[NGL], p1xsq[NGL] +real p1p2, dhdsx, dhdsy, erfx, erfy, p1sq, p2sq, y, ysq, p2ysq, profile +real xy, wt, denom, alpha, func, p4fod, wp4fod, f, e, rsq, wf +real wfsq, onemp3, dfby, deby, dbyx0, dbyy0 +real daoerf() + +data d / 0.0, 0.0, 0.0, 0.0, + -0.28867513, 0.28867513, 0.0, 0.0, + -0.38729833, 0.0, 0.38729833, 0.0, + -0.43056816, -0.16999052, 0.16999052, 0.43056816 / +data w / 1.0, 0.0, 0.0, 0.0, + 0.5, 0.5, 0.0, 0.0, + 0.27777778, 0.44444444, 0.27777778, 0.0, + 0.17392742, 0.32607258, 0.32607258, 0.17392742 / +begin + # Initialize. + profile = 0.0 + dhdxc = 0.0 + dhdyc = 0.0 + + # Compute the analytic part of the profile for a given x and y. + + switch (ipstyp) { + + # Evaluate the Gaussian function + # f = erfx * erfy / (par[1] * par[2]) + # par[1] is the hwhm in x; sigma(x) = 0.8493218 * hwhm + # par[2] is the hwhm in y; sigma(y) = 0.8493218 * hwhm + + case FCTN_GAUSS: + + p1p2 = par[1] * par[2] + erfx = daoerf (dx, 0.0, par[1], dhdxc, dhdsx) + erfy = daoerf (dy, 0.0, par[2], dhdyc, dhdsy) + profile = erfx * erfy / p1p2 + dhdxc = dhdxc * erfy / p1p2 + dhdyc = dhdyc * erfx / p1p2 + if (ideriv > 0) { + term[1] = (dhdsx - erfx / par[1]) * erfy / p1p2 + term[2] = (dhdsy - erfy / par[2]) * erfx / p1p2 + } + + # Evaluate the Moffat25 function + # + # f = (beta - 1) / (ax * ay * (1 + (x/ax) ** 2 + (y/ay) ** 2 + + # (xy * axy)) ** beta) + # + # par[1] is the hwhm in x at y = 0.0 + # 1/2 = 1 / (1 + (par[1] / ax) ** 2) ** beta + # so + # 2 ** (1/ beta) - 1 = (par[1] / ax) ** 2 + # ax ** 2 = par[1] ** 2 / (2 ** (1 / beta) - 1) + # + # when beta = 2.5 ax ** 2 = 3.129813 * par[1] ** 2 + # + # f = 1 / par[1] * par[2] * (1 + 0.3195079 * ((x / par[1]) ** 2 + + # (y / par[2]) ** 2 + xy * par[3])) ** 2.5 + # + + case FCTN_MOFFAT25: + + alpha = 0.3195079 + #talpha = 0.6390158 + p1sq = par[1] ** 2 + p2sq = par[2] ** 2 + p1p2 = par[1] * par[2] + xy = dx * dy + if (ideriv > 0) + call aclrr (term, 4) + + denom = 1.0 + alpha * (dx ** 2 / p1sq + dy ** 2 / p2sq + xy * + par[3]) + if (denom > 1.0e4) + return (profile) + func = 1.0 / (p1p2 * denom ** par[4]) + if (func >= 0.046) { + npt = 4 + } else if (func >= 0.0022) { + npt = 3 + } else if (func >= 0.0001) { + npt = 2 + } else if (func >= 1.0e-10) { + profile = (par[4] - 1.0) * func + p4fod = par[4] * alpha * profile / denom + dhdxc = p4fod * (2.0 * dx / p1sq + dy * par[3]) + dhdyc = p4fod * (2.0 * dy / p2sq + dx * par[3]) + if (ideriv > 0) { + term[1] = (2.0 * p4fod * dx ** 2 / p1sq - profile) / + par[1] + term[2] = (2.0 * p4fod * dy ** 2 / p2sq - profile) / + par[2] + term[3] = - p4fod * xy + term[4] = profile * (1.0 / (par[4] - 1.0) - log (denom)) + } + return (profile) + } else { + return (profile) + } + + do ix = 1, npt { + x[ix] = dx + d[ix,npt] + xsq[ix] = x[ix] ** 2 + p1xsq[ix] = xsq[ix] / p1sq + } + + do iy = 1, npt { + y = dy + d[iy,npt] + ysq = y ** 2 + p2ysq = ysq / p2sq + do ix = 1, npt { + wt = w[iy,npt] * w[ix,npt] + xy = x[ix] * y + denom = 1.0 + alpha * (p1xsq[ix] + p2ysq + xy * par[3]) + func = (par[4] - 1.0) / (p1p2 * denom ** par[4]) + p4fod = par[4] * alpha * func / denom + wp4fod = wt * p4fod + wf = wt * func + profile = profile + wf + dhdxc = dhdxc + wp4fod * (2.0 * x[ix] / p1sq + + y * par[3]) + dhdyc = dhdyc + wp4fod * (2. * y / p2sq + x[ix] * + par[3]) + if (ideriv > 0) { + term[1] = term[1] + (2.0 * wp4fod * p1xsq[ix] - wf) / + par[1] + term[2] = term[2] + (2.0 * wp4fod * p2ysq - wf) / + par[2] + term[3] = term[3] - wp4fod * xy + term[4] = term[4] + wf * (1.0 / (par[4] - 1.0) - + log (denom)) + } + } + } + + # + # Penny function which has a gaussian core and lorentzian wings. + # The lorentzian is elongated along the x and y axes. + + case FCTN_PENNY1: + + p1sq = par[1] ** 2 + p2sq = par[2] ** 2 + onemp3 = 1.0 - par[3] + xy = dx * dy + if (ideriv > 0) + call aclrr (term, 4) + + rsq = dx ** 2 / p1sq + dy ** 2 / p2sq + if (rsq > 1.0e10) + return (profile) + + f = 1.0 / (1.0 + rsq) + rsq = rsq + xy * par[4] + if (rsq < 34.0) { + e = exp (-0.6931472 * rsq) + func = par[3] * e + onemp3 * f + } else { + e = 0.0 + func = onemp3 * f + } + + if (func >= 0.046) { + npt = 4 + } else if (func >= 0.0022) { + npt = 3 + } else if (func >= 0.0001) { + npt = 2 + } else if (func >= 1.0e-10) { + profile = func + dfby = onemp3 * f ** 2 + deby = 0.6931472 * par[3] * e + dbyx0 = 2.0 * dx / p1sq + dbyy0 = 2.0 * dy / p2sq + dhdxc = deby * (dbyx0 + dy * par[4]) + dfby * dbyx0 + dhdyc = deby * (dbyy0 + dx * par[4]) + dfby * dbyy0 + if (ideriv > 0) { + dbyx0 = dbyx0 * dx / par[1] + dbyy0 = dbyy0 * dy / par[2] + dfby = dfby + deby + term[1] = dfby * dbyx0 + term[2] = dfby * dbyy0 + term[3] = e - f + term[4] = - deby * xy / (0.5 - abs(par[4])) + } + return (profile) + } else { + return (profile) + } + + do ix = 1, npt { + x[ix] = dx + d[ix,npt] + p1xsq[ix] = x[ix] / p1sq + } + + do iy = 1, npt { + y = dy + d[iy,npt] + p2ysq = y / p2sq + do ix = 1, npt { + wt = w[iy,npt] * w[ix,npt] + xy = x[ix] * y + rsq = p1xsq[ix] * x[ix] + p2ysq * y + f = 1.0 / (1.0 + rsq) + rsq = rsq + xy * par[4] + if (rsq < 34.0) { + e = exp (- 0.6931472 * rsq) + func = par[3] * e + onemp3 * f + deby = 0.6931472 * wt * par[3] * e + } else { + e = 0.0 + func = onemp3 * f + deby = 0.0 + } + profile = profile + wt * func + dfby = wt * onemp3 * f ** 2 + dbyx0 = 2.0 * p1xsq[ix] + dbyy0 = 2.0 * p2ysq + dhdxc = dhdxc + deby * (dbyx0 + dy * par[4]) + dfby * dbyx0 + dhdyc = dhdyc + deby * (dbyy0 + dx * par[4]) + dfby * dbyy0 + if (ideriv > 0) { + dbyx0 = dbyx0 * dx / par[1] + dbyy0 = dbyy0 * dy / par[2] + term[1] = term[1] + (dfby + deby) * dbyx0 + term[2] = term[2] + (dfby + deby) * dbyy0 + term[3] = term[3] + wt * (e - f) + term[4] = term[4] - deby * xy + } + } + } + + # Penny function which has a gaussian core and lorentzian wings. + # The gaussian and lorentzian may be tilted in different directions. + + case FCTN_PENNY2: + + p1sq = par[1] ** 2 + p2sq = par[2] ** 2 + onemp3 = 1.0 - par[3] + xy = dx * dy + if (ideriv > 0) + call aclrr (term, 5) + + rsq = dx ** 2 / p1sq + dy ** 2 / p2sq + dfby = rsq + par[5] * xy + if (dfby > 1.0e10) + return (profile) + f = 1.0 / (1.0 + dfby) + + deby = rsq + par[4] * xy + if (deby < 34.0) + e = exp (-0.6931472 * deby) + else + e = 0.0 + + func = par[3] * e + onemp3 * f + if (func >= 0.046) { + npt = 4 + } else if (func >= 0.0022) { + npt = 3 + } else if (func >= 0.0001) { + npt = 2 + } else if (func >= 1.0e-10) { + profile = func + dfby = onemp3 * f ** 2 + deby = 0.6931472 * par[3] * e + dbyx0 = 2.0 * dx / p1sq + dbyy0 = 2.0 * dy / p2sq + dhdxc = deby * (dbyx0 + dy * par[4]) + dfby * (dbyx0 + dy * + par[5]) + dhdyc = deby * (dbyy0 + dx * par[4]) + dfby * (dbyy0 + dx * + par[5]) + if (ideriv > 0) { + dbyx0 = dbyx0 * dx / par[1] + dbyy0 = dbyy0 * dy / par[2] + term[5] = -dfby * xy + dfby = dfby + deby + term[1] = dfby * dbyx0 + term[2] = dfby * dbyy0 + term[3] = e - f + term[4] = - deby * xy + } + return (profile) + } else { + return (profile) + } + + do ix = 1, npt { + x[ix] = dx + d[ix,npt] + p1xsq[ix] = x[ix] / p1sq + } + + do iy = 1, npt { + y = dy + d[iy,npt] + p2ysq = y / p2sq + do ix = 1, npt { + wt = w[iy,npt] * w[ix,npt] + xy = x[ix] * y + rsq = p1xsq[ix] * x[ix] + p2ysq * y + f = 1.0 / (1.0 + rsq + par[5] * xy) + deby = rsq + par[4] * xy + if (deby < 34.0) { + e = exp (- 0.6931472 * deby) + func = par[3] * e + onemp3 * f + deby = 0.6931472 * wt * par[3] * e + } else { + e = 0.0 + func = onemp3 * f + deby = 0.0 + } + profile = profile + wt * func + dfby = wt * onemp3 * f ** 2 + dbyx0 = 2.0 * p1xsq[ix] + dbyy0 = 2.0 * p2ysq + dhdxc = dhdxc + deby * (dbyx0 + dy * par[4]) + dfby * + (dbyx0 + dy * par[5]) + dhdyc = dhdyc + deby * (dbyy0 + dx * par[4]) + dfby * + (dbyy0 + dx * par[5]) + if (ideriv > 0) { + dbyx0 = dbyx0 * dx / par[1] + dbyy0 = dbyy0 * dy / par[2] + term[1] = term[1] + (dfby + deby) * dbyx0 + term[2] = term[2] + (dfby + deby) * dbyy0 + term[3] = term[3] + wt * (e - f) + term[4] = term[4] - deby * xy + term[5] = term[5] - dfby * xy + } + } + } + + # Evaluate the Moffat15 function + # + # f = (beta - 1) / (ax * ay * (1 + (x/ax) ** 2 + (y/ay) ** 2 + + # (xy * axy)) ** beta) + # + # par[1] is the hwhm in x at y = 0.0 + # 1/2 = 1 / (1 + (par[1] / ax) ** 2) ** beta + # so + # 2 ** (1/ beta) - 1 = (par[1] / ax) ** 2 + # ax ** 2 = par[1] ** 2 / (2 ** (1 / beta) - 1) + # + # when beta = 1.5 ax ** 2 = 1.7024144 * par[1] ** 2 + # + # f = 1 / par[1] * par[2] * (1 + 0.5874011 * ((x / par[1]) ** 2 + + # (y / par[2]) ** 2 + xy * par[3])) ** 2.5 + # + + case FCTN_MOFFAT15: + + alpha = 0.5874011 + #talpha = 1.1748021 + p1sq = par[1] ** 2 + p2sq = par[2] ** 2 + p1p2 = par[1] * par[2] + xy = dx * dy + if (ideriv > 0) + call aclrr (term, 4) + + denom = 1.0 + alpha * (dx ** 2 / p1sq + dy ** 2 / p2sq + xy * + par[3]) + if (denom > 5.0e6) + return (profile) + func = 1.0 / (p1p2 * denom ** par[4]) + if (func >= 0.046) { + npt = 4 + } else if (func >= 0.0022) { + npt = 3 + } else if (func >= 0.0001) { + npt = 2 + } else if (func >= 1.0e-10) { + profile = (par[4] - 1.0) * func + p4fod = par[4] * alpha * profile / denom + dhdxc = p4fod * (2.0 * dx / p1sq + dy * par[3]) + dhdyc = p4fod * (2.0 * dy / p2sq + dx * par[3]) + if (ideriv > 0) { + term[1] = (2.0 * p4fod * dx ** 2 / p1sq - profile) / + par[1] + term[2] = (2.0 * p4fod * dy ** 2 / p2sq - profile) / + par[2] + term[3] = - p4fod * xy + term[4] = profile * (1.0 / (par[4] - 1.0) - log (denom)) + } + return (profile) + } else { + return (profile) + } + + do ix = 1, npt { + x[ix] = dx + d[ix,npt] + xsq[ix] = x[ix] ** 2 + p1xsq[ix] = xsq[ix] / p1sq + } + + do iy = 1, npt { + y = dy + d[iy,npt] + ysq = y ** 2 + p2ysq = ysq / p2sq + do ix = 1, npt { + wt = w[iy,npt] * w[ix,npt] + xy = x[ix] * y + denom = 1.0 + alpha * (p1xsq[ix] + p2ysq + xy * + par[3]) + func = (par[4] - 1.0) / (p1p2 * denom ** par[4]) + p4fod = par[4] * alpha * func / denom + wp4fod = wt * p4fod + wf = wt * func + profile = profile + wf + dhdxc = dhdxc + wp4fod * (2.0 * x[ix] / p1sq + y * + par[3]) + dhdyc = dhdyc + wp4fod * (2. * y / p2sq + x[ix] * + par[3]) + if (ideriv > 0) { + term[1] = term[1] + (2.0 * wp4fod * p1xsq[ix] - wf) / + par[1] + term[2] = term[2] + (2.0 * wp4fod * p2ysq - wf) / + par[2] + term[3] = term[3] - wp4fod * xy + term[4] = term[4] + wf * (1.0 / (par[4] - 1.0) - + log (denom)) + } + } + } + + case FCTN_LORENTZ: + + p1sq = par[1] ** 2 + p2sq = par[2] ** 2 + p1p2 = par[1] * par[2] + xy = dx * dy + if (ideriv > 0) + call aclrr (term, 3) + + denom = 1.0 + dx ** 2 / p1sq + dy ** 2 / p2sq + xy * par[3] + if (denom > 1.0e10) + return (profile) + func = 1.0 / denom + if (func >= 0.046) { + npt = 4 + } else if (func >= 0.0022) { + npt = 3 + } else if (func >= 0.0001) { + npt = 2 + } else if (func >= 1.0e-10) { + profile = func + wfsq = func ** 2 + dhdxc = wfsq * (2.0 * dx / p1sq + dy * par[3]) + dhdyc = wfsq * (2.0 * dy / p2sq + dx * par[3]) + if (ideriv > 0) { + term[1] = wfsq * (2.0 * dx ** 2 / p1sq) / par[1] + term[2] = wfsq * (2.0 * dy ** 2 / p2sq) / par[2] + term[3] = - wfsq * xy + } + return (profile) + } else { + return (profile) + } + + do ix = 1, npt { + x[ix] = dx + d[ix,npt] + xsq[ix] = x[ix] ** 2 + p1xsq[ix] = xsq[ix] / p1sq + } + + do iy = 1, npt { + y = dy + d[iy,npt] + ysq = y ** 2 + p2ysq = ysq / p2sq + do ix = 1, npt { + wt = w[iy,npt] * w[ix,npt] + xy = x[ix] * y + denom = 1.0 + p1xsq[ix] + p2ysq + xy * par[3] + func = 1.0 / denom + wf = wt * func + wfsq = wf * func + profile = profile + wf + dhdxc = dhdxc + wfsq * (2.0 * x[ix] / p1sq + y * par[3]) + dhdyc = dhdyc + wfsq * (2.0 * y / p2sq + x[ix] * par[3]) + if (ideriv > 0) { + term[1] = term[1] + wfsq * (2.0 * p1xsq[ix]) / par[1] + term[2] = term[2] + wfsq * (2.0 * p2ysq) / par[2] + term[3] = term[3] - wfsq * xy + } + } + } + + default: + profile = INDEFR + } + + return (profile) +end diff --git a/noao/digiphot/daophot/daolib/quick.f b/noao/digiphot/daophot/daolib/quick.f new file mode 100644 index 00000000..72a58fe5 --- /dev/null +++ b/noao/digiphot/daophot/daolib/quick.f @@ -0,0 +1,202 @@ +c subroutine quick (dataum, n, index) +c +c +c A quick-sorting algorithm suggested by the discussion on pages 114-119 +c of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND SEARCHING, by +c D.E. Knuth, which was referenced in Don Wells' subroutine QUIK. This +c is my own attempt at encoding a quicksort-- PBS. +c +c Arguments +c +c datum (input / output) is a vector of dimension n containing randomly +c ordered real data upon input. Upon output the elements of +c dataum will be in order of increasing value. +c +c +c index (output) is an integer vector of dimension n. Upon return to +c the calling program the i-th element of index will tell where +c the i-th element of the sorted vector datum had been before +c datum was sorted. +c +c +c +c parameter (maxstack = 28) +c +c Parameters +c +c maxstack is the maximum number of entries the stack can contain. +c A limiting stack length of 28 restricts this quicksort + + subroutine quick (datum, n, index, ier) +c + implicit none + integer n, index(n), ier + real datum(n) +c + real dkey + integer stklo(28), stkhi(28), i, lo, hi, nstak, limlo, limhi, ikey +c +c Initialize error code + + ier = 0 +c +c Initialize index. +c + do 50 i = 1, n + 50 index(i) = i +c +c Initialize the pointers. +c + nstak = 0 + limlo = 1 + limhi = n +c + 100 dkey = datum(limlo) + ikey = index(limlo) +c +c Compare all elements in the sub-vector between limlo and limhi with +c the current key datum. +c + lo = limlo + hi = limhi + 101 continue +c + if (lo .eq. hi) go to 200 +c + if (datum(hi) .le. dkey) go to 109 + hi = hi - 1 +c +c The pointer hi is to be left pointing at a datum smaller than the +c key, which is intended to be overwritten. +c + go to 101 +c + 109 datum(lo) = datum(hi) + index(lo) = index(hi) + lo = lo + 1 + 110 continue +c + if (lo .eq. hi) go to 200 +c + if (datum(lo) .ge. dkey) go to 119 +c + lo = lo + 1 + go to 110 +c + 119 datum(hi) = datum(lo) + index(hi) = index(lo) + hi = hi - 1 +c +c The pointer LO is to be left pointing at a datum LARGER than the +c key, which is intended to be overwritten. +c + go to 101 +c + 200 continue +c +c lo and hi are equal, and point at a value which is intended to +c be overwritten. Since all values below this point are less than +c the key and all values above this point are greater than the key, +c this is where we stick the key back into the vector. +c + datum(lo) = dkey + index(lo) = ikey +c +c At this point in the subroutine, all data between limlo and LO-1, +c inclusive, are less than datum(LO), and all data between LO+1 and +c limhi are larger than datum(LO). +c +c If both subarrays contain no more than one element, then take the most +c recent interval from the stack (if the stack is empty, we're done). +c If the larger of the two subarrays contains more than one element, and +c if the shorter subarray contains one or no elements, then forget the +c shorter one and reduce the other subarray. If the shorter subarray +c contains two or more elements, then place the larger subarray on the +c stack and process the subarray. +c + if ((limhi - lo) .gt. (lo - limlo)) go to 300 +c +c Case 1: the lower subarray is longer. If it contains one or no +c elements then take the most recent interval from the stack and go +c back and operate on it. +c + if ((lo - limlo) .le. 1) go to 400 +c +c If the upper (shorter) subinterval contains one or no elements, then +c process the lower (longer) one, but if the upper subinterval contains +c more than one element, then place the lower (longer) subinterval on +c the stack and process the upper one. +c + if ((limhi - lo) .ge. 2) go to 250 +c +c Case 1a: the upper (shorter) subinterval contains no or one elements, +c so we go back and operate on the lower (longer) subinterval. +c + limhi = lo - 1 + go to 100 +c + 250 continue +c +c Case 1b: the upper (shorter) subinterval contains at least two +c elements, so we place the lower (longer) subinterval on the stack and +c then go back and operate on the upper subinterval. +c + nstak = nstak + 1 + if (nstak .gt. 28) then + ier = -1 + return + endif + stklo(nstak) = limlo + stkhi(nstak) = lo - 1 + limlo = lo + 1 + go to 100 +c + 300 continue +c +c Case 2: the upper subarray is longer. If it contains one or no +c elements then take the most recent interval from the stack and +c operate on it. +c + if ((limhi - lo) .le. 1) go to 400 +c +c If the lower (shorter) subinterval contains one or no elements, then +c process the upper (longer) one, but if the lower subinterval contains +c more than one element, then place the upper (longer) subinterval on +c the stack and process the lower one. +c + if ((lo - limlo) .ge. 2) go to 350 +c +c Case 2a: the lower (shorter) subinterval contains no or one elements, +c so we go back and operate on the upper (longer) subinterval. +c + limlo = lo + 1 + go to 100 +c + 350 continue +c +c Case 2b: the lower (shorter) subinterval contains at least two +c elements, so we place the upper (longer) subinterval on the stack and +c then go back and operate on the lower subinterval. +c + nstak = nstak + 1 + if (nstak .gt. 28) then + ier = -1 + return + endif + stklo(nstak) = lo + 1 + stkhi(nstak) = limhi + limhi = lo - 1 + go to 100 +c + 400 continue +c +c Take the most recent interval from the stack. If the stack happens +c to be empty, we are done. +c + if (nstak .le. 0) return + limlo = stklo(nstak) + limhi = stkhi(nstak) + nstak = nstak - 1 + go to 100 +c + end diff --git a/noao/digiphot/daophot/daolib/ran3.x b/noao/digiphot/daophot/daolib/ran3.x new file mode 100644 index 00000000..c8915e70 --- /dev/null +++ b/noao/digiphot/daophot/daolib/ran3.x @@ -0,0 +1,63 @@ +define MBIG 1000000000 +define MSEED 161803398 +define MZ 0.0 +define FAC 1.0 / MBIG + +# RAN3 -- Returns a uniform random deviate between 0.0 and 1.0. Set +# 'idum' to any negative value to initialize or reinitialize the sequence. +# From Numerical Recipes (originally attributed to Donald Knuth, 1981; +# Seminumerical Algorithms, 2nd edition, volume 2 of 'The Art of Computer +# Programming' - Section 3.2-3.3. + +real procedure ran3 (idum) + +int idum + +int ma[55] +int mj, mk, i, k, ii +int iff, inext, inextp +data iff /0/ + +begin + if(idum < 0 || iff == 0) { + iff = 1 + mj = MSEED - iabs(idum) + mj = mod(mj, MBIG) + ma[55] = mj + mk = 1 + + do i = 1, 54 { + ii = mod(21 * i , 55) + ma[ii] = mk + mk = mj - mk + if (mk < MZ) + mk = mk + MBIG + mj = ma[ii] + } + + do k = 1, 4 { + do i = 1, 55 { + ma[i] = ma[i] - ma[1+mod(i+30, 55)] + if (ma[i] < MZ) + ma[i] = ma[i] + MBIG + } + } + + inext = 0 + inextp = 31 + idum = 1 + } + + inext = inext + 1 + if (inext == 56) + inext = 1 + inextp = inextp + 1 + if (inextp == 56) + inextp = 1 + mj = ma[inext] - ma[inextp] + if (mj < MZ) + mj = mj + MBIG + ma[inext]= mj + return (mj * FAC) + +end diff --git a/noao/digiphot/daophot/daolib/usepsf.x b/noao/digiphot/daophot/daolib/usepsf.x new file mode 100644 index 00000000..f7a7db2f --- /dev/null +++ b/noao/digiphot/daophot/daolib/usepsf.x @@ -0,0 +1,81 @@ +include "../lib/daophotdef.h" + +# DP_USEPSF -- Evaluate the psf at a given point. + +real procedure dp_usepsf (ipstyp, dx, dy, bright, par, psf, npsf, nexp, + nfrac, deltax, deltay, dvdxc, dvdyc) + +int ipstyp # analytic psf function type +real dx, dy # distance for center of function +real bright # the relative brightness of the object +real par[ARB] # current values of the parameters +real psf[npsf,npsf,ARB] # the psf lookup tables +int npsf # size of the psf look-up table +int nexp # number pf look-up tables +int nfrac # fractional pixel expansion +real deltax, deltay # distance from center of look-up tables +real dvdxc, dvdyc # derivatives with respect to position + +int nterm, j, k, lx, ly +real psfval, middle, junk[MAX_NEXPTERMS], xx, yy, corr, dfdx, dfdy +real dp_profile(), bicubic() + +begin + nterm = nexp + nfrac + psfval = bright * dp_profile (ipstyp, dx, dy, par, dvdxc, dvdyc, + junk, 0) + dvdxc = bright * dvdxc + dvdyc = bright * dvdyc + if (nterm <= 0) + return (psfval) + + # The PSF look-up tables are centered at (MIDDLE, MIDDLE). + + switch (nexp) { + case 1: + junk[1] = 1. + case 3: + junk[1] = 1. + junk[2] = deltax + junk[3] = deltay + case 6: + junk[1] = 1. + junk[2] = deltax + junk[3] = deltay + junk[4] = 1.5 * deltax ** 2 - 0.5 + junk[5] = deltax * deltay + junk[6] = 1.5 * deltay ** 2 - 0.5 + } + + if (nfrac > 0) { + j = nexp + 1 + junk[j] = - 2. * (dx - real(nint(dx))) + j = j + 1 + junk[j] = - 2. * (dy - real(nint(dy))) + j = j + 1 + junk[j] = 1.5 * junk[j-2] ** 2 - 0.5 + j = j + 1 + junk[j] = junk[j-3] * junk[j-2] + j = j + 1 + junk[j] = 1.5 * junk[j-3] ** 2 - 0.5 + } + + # This point in the stellar profile lies between columns LX and LX+1, + # and between rows LY and LY+1 in the look-up tables. + + middle = (npsf + 1) / 2 + xx = (2. * dx) + middle + lx = int (xx) + yy = (2. * dy) + middle + ly = int (yy) + + do k = 1, nterm { + corr = bicubic (psf[lx-1,ly-1,k], npsf, xx - real(lx), + yy - real(ly), dfdx, dfdy) + psfval = psfval + junk[k] * corr + dvdxc = dvdxc - junk[k] * dfdx + dvdyc = dvdyc - junk[k] * dfdy + } + + return (psfval) +end |