diff options
Diffstat (limited to 'noao/digiphot/photcal/mkobsfile/phaigrow.x')
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phaigrow.x | 1635 |
1 files changed, 1635 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/mkobsfile/phaigrow.x b/noao/digiphot/photcal/mkobsfile/phaigrow.x new file mode 100644 index 00000000..7fa32ce0 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phaigrow.x @@ -0,0 +1,1635 @@ +include <time.h> +include <mach.h> +include <gset.h> +include "../lib/apfile.h" + +define HELPFILE "photcal$mkobsfile/apfile.key" + +# PH_AIGROW - Compute the curve of growth using the data from the input +# photometry files, optional airmasses from the obsparams file, and the +# initial values of the parameters using interactive graphics. + +procedure ph_aigrow (gd, apcat, magfd, logfd, mgd, imtable, id, x, y, nap, + rap, mag, merr, naperts, params, lterms, smallap, largeap) + +pointer gd # pointer to the graphics descriptor +int apcat # the aperture correction file descriptor +int magfd # the best magnitudes file descriptor +int logfd # the output log file descriptor +pointer mgd # pointer to the metacode graphics stream +pointer imtable # pointer to the image symbol table +int id[ARB] # the star ids +real x[ARB] # the star x coordinates +real y[ARB] # the star y coordinates +int nap[ARB] # the array of aperture numbers +real rap[naperts,ARB] # input array of aperture radii +real mag[naperts,ARB] # input array of magnitudes / differences +real merr[naperts,ARB] # input array of magnitude errors +int naperts # the number of input apertures +double params[ARB] # the initial values for the parameters +int lterms # the number of terms to be fit +int smallap # the small aperture number +int largeap # the large aperture number + +int k, nimages, newfit, fitok, wcs, key, isymbol, newgraph, newimage +int graphtype, ndata +pointer sp, sym, symbol, agr, cmd, inap +real wx, wy +int stnsymbols(), ph_amfit(), clgcur(), ph_audelete() +pointer sthead(), stnext() + +begin + if (logfd != NULL) + call ph_logtitle (logfd, apcat) + + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) { + call printf ("Error: There is no input data to fit\n") + if (logfd != NULL) + call fprintf (logfd, "Error: There is no input data to fit\n") + return + } + + # Allocate temporary working space + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Order the symbol table. + symbol = sthead (imtable) + ndata = 0 + do k = nimages, 1, -1 { + Memi[sym+k-1] = symbol + ndata = ndata + IMT_NENTRIES(symbol) + symbol = stnext (imtable, symbol) + } + + # Allocate memory required for fitting. + call ph_ameminit (agr, lterms, naperts) + + # Do the initial a parameter and seeing fit. + if (ph_amfit (imtable, agr, nap, rap, mag, merr, naperts, params, + lterms) == ERR) + fitok = NO + else + fitok = YES + newfit = NO + + # Define the current image to be the first image and evaluate the fit. + symbol = Memi[sym] + isymbol = 1 + call ph_a1eval (agr, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_NXAIRMASS(symbol), nap[IMT_OFFSET(symbol)], + rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + merr[1,IMT_OFFSET(symbol)], naperts, IMT_NENTRIES(symbol)) + newimage = NO + + # Plot the initial graph. + call ph_gimfit (gd, agr, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_XAIRMASS(symbol), nap[IMT_OFFSET(symbol)], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, IMT_NENTRIES(symbol), fitok) + graphtype = AGR_FIT + newgraph = NO + + # Print the status report. + call ph_aeprint (agr, IMT_IMNAME(symbol), IMT_RO(symbol), + rap[1,IMT_OFFSET(symbol)], naperts, smallap, largeap, fitok, + newfit) + + # Allocate the rejection array + call malloc (inap, ndata, TY_INT) + call amovi (nap, Memi[inap], ndata) + + # Examine the fit. + while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + switch (key) { + + # Quit the program. + case 'q': + break + + # Print help + case '?': + call gpagefile (gd, HELPFILE, "") + + # Execute colon commands. + case ':': + call ph_acolon (gd, imtable, agr, symbol, isymbol, params, + lterms, naperts, smallap, largeap, Memc[cmd], fitok, newfit, + newimage, newgraph) + + # Compute a new fit. + case 'f': + if (newfit == YES) { + call ph_amemfree (agr) + call ph_ameminit (agr, lterms, naperts) + } + if (ph_amfit (imtable, agr, Memi[inap], rap, mag, merr, naperts, + params, lterms) == ERR) + fitok = NO + else + fitok = YES + newfit = NO + newimage = YES + newgraph = YES + isymbol = 1 + symbol = Memi[sym+isymbol-1] + + + # Print out the id of a particular object + case 'c': + call ph_pnearest (gd, graphtype, wx, wy, agr, + id[IMT_OFFSET(symbol)], x[IMT_OFFSET(symbol)], + y[IMT_OFFSET(symbol)], Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, IMT_NENTRIES(symbol), + smallap, largeap) + + # Write the answer for the current image. + case 'w': + call ph_aeprint (agr, IMT_IMNAME(symbol), IMT_RO(symbol), + rap[1,IMT_OFFSET(symbol)], naperts, smallap, largeap, + fitok, newfit) + + # Delete points. + case 'd': + if (ph_audelete (gd, graphtype, wx, wy, agr, + x[IMT_OFFSET(symbol)], y[IMT_OFFSET(symbol)], + Memi[inap+IMT_OFFSET(symbol)-1], nap[IMT_OFFSET(symbol)], + rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + naperts, IMT_NENTRIES(symbol), smallap, largeap, YES) == + YES) + newfit = YES + + # Undelete points. + case 'u': + if (ph_audelete (gd, graphtype, wx, wy, agr, + x[IMT_OFFSET(symbol)], y[IMT_OFFSET(symbol)], + Memi[inap+IMT_OFFSET(symbol)-1], nap[IMT_OFFSET(symbol)], + rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + naperts, IMT_NENTRIES(symbol), smallap, largeap, NO) == + YES) + newfit = YES + + # Redraw the graph. + case 'g': + newgraph = YES + + # Graph the model fit. + case 'm': + if (graphtype != AGR_FIT) { + graphtype = AGR_FIT + newgraph = YES + } + + # Graph the resdiduals from the adopted model as a function of + # aperture. + case 'r': + if (fitok == YES && ! IS_INDEFR(IMT_RO(symbol)) && + graphtype != AGR_ARESIDUALS) { + graphtype = AGR_ARESIDUALS + newgraph = YES + } + + # Graph the residuals from the adopted model as a function of + # magnitude. + case 'b': + if (fitok == YES && ! IS_INDEFR(IMT_RO(symbol)) && + graphtype != AGR_BRESIDUALS) { + graphtype = AGR_BRESIDUALS + newgraph = YES + } + + # Graph the residuals as a fucntion of x. + case 'x': + if (fitok == YES && ! IS_INDEFR(IMT_RO(symbol)) && + graphtype != AGR_XRESIDUALS) { + graphtype = AGR_XRESIDUALS + newgraph = YES + } + + # Graph the residuals as a function of y. + case 'y': + if (fitok == YES && ! IS_INDEFR(IMT_RO(symbol)) && + graphtype != AGR_YRESIDUALS) { + graphtype = AGR_YRESIDUALS + newgraph = YES + } + + # Graph the cumulative aperture correction. + case 'a': + if (fitok == YES && ! IS_INDEFR(IMT_RO(symbol)) && + graphtype != AGR_CUMULATIVE) { + graphtype = AGR_CUMULATIVE + newgraph = YES + } + + # Move to the previous image. + case 'p': + if (isymbol == 1) + call printf ("Already at beginning of image list\n") + else { + isymbol = isymbol - 1 + symbol = Memi[sym+isymbol-1] + newimage = YES + newgraph = YES + } + + # Move to the next image. + case 'n': + if (isymbol == nimages) + call printf ("Already at end of image list\n") + else { + isymbol = isymbol + 1 + symbol = Memi[sym+isymbol-1] + newimage = YES + newgraph = YES + } + } + + if (newimage == YES) { + call ph_a1eval (agr, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_NXAIRMASS(symbol), Memi[inap+IMT_OFFSET(symbol)-1], + rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + merr[1,IMT_OFFSET(symbol)], naperts, IMT_NENTRIES(symbol)) + newimage = NO + } + + if (newgraph == YES) { + switch (graphtype) { + case AGR_FIT: + call ph_gimfit (gd, agr, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_XAIRMASS(symbol), Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), fitok) + case AGR_ARESIDUALS: + call ph_gaimres (gd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), fitok) + case AGR_BRESIDUALS: + call ph_gbimres (gd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + naperts, IMT_NENTRIES(symbol), fitok) + case AGR_XRESIDUALS: + call ph_gaximres (gd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], x[IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), fitok) + case AGR_YRESIDUALS: + call ph_gayimres (gd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], y[IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), fitok) + case AGR_CUMULATIVE: + call ph_gacum (gd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), smallap, largeap, fitok) + default: + call ph_gimfit (gd, agr, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_XAIRMASS(symbol), Memi[inap+IMT_OFFSET(symbol)-1], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), fitok) + } + call ph_aeprint (agr, IMT_IMNAME(symbol), IMT_RO(symbol), + rap[1,IMT_OFFSET(symbol)], naperts, smallap, largeap, + fitok, newfit) + newgraph = NO + } + } + + # Output the results. + + # Fit the seeing radii and the cog model parameters. + if (fitok == NO) { + call printf ("Error: The cog model fit did not converge\n") + if (logfd != NULL) { + call ph_ishow (logfd, params, lterms) + call ph_pshow (logfd, agr, lterms) + call ph_rshow (logfd, imtable) + call fprintf (logfd, + "Error: The cog model fit did not converge\n") + } + } else { + if (newfit == YES) { + call printf ("Warning: The cog model fit is out-of-date\n") + if (logfd != NULL) + call fprintf (logfd, + "Warning: The cog model fit is out-of-date\n") + } + if (logfd != NULL) { + call ph_ishow (logfd, params, lterms) + call ph_pshow (logfd, agr, lterms) + call ph_rshow (logfd, imtable) + } + call ph_aeval (imtable, agr, apcat, magfd, logfd, mgd, id, x, y, + Memi[inap], rap, mag, merr, naperts, smallap, largeap) + if (logfd != NULL) + call ph_tfshow (logfd, agr, naperts) + } + + # Free memory required for fitting. + call ph_amemfree (agr) + + call mfree (inap, TY_INT) + call sfree (sp) +end + + +# PH_GIMFIT -- Graph the data and fit for a particular image. + +procedure ph_gimfit (gd, agr, image, r0, xairmass, inap, nap, rap, mag, naperts, + npts, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good aperture numbers +int nap[ARB] # the array of aperture numbers +real rap[naperts,ARB] # input array of aperture radii +real mag[naperts,ARB] # input array of magnitudes / differences +int naperts # the number of input apertures +int npts # the number of points +int fitok # is the fit ok + +int i, j +pointer sp, title +real rmin, rmax, dr, mmin, mmax, dm +int strlen() +real ph_achi() + +begin + if (gd == NULL) + return + + # Allocate working space. + call smark (sp) + call salloc (title, 2 * SZ_LINE, TY_CHAR) + + # Compute the data window in x. + call alimr (rap, naperts * npts, rmin, rmax) + dr = rmax - rmin + rmin = rmin - 0.1 * dr + rmax = rmax + 0.1 * dr + + # Compute the data window in y. + mmin = MAX_REAL + mmax = -MAX_REAL + do i = 1, npts { + do j = 2, nap[i] { + if (mag[j,i] < mmin) + mmin = mag[j,i] + if (mag[j,i] > mmax) + mmax = mag[j,i] + } + } + if (mmin > mmax) { + mmin = -0.1 + mmax = 0.1 + } else { + dm = mmax - mmin + mmin = mmin - 0.1 * dm + mmax = mmax + 0.1 * dm + } + + # Clear the screen and set the data window. + call gclear (gd) + call gswind (gd, rmin, rmax, mmin, mmax) + + # Set up the title and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], + 2 * SZ_LINE - strlen(Memc[title]), + "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + if (fitok == YES && ! IS_INDEFR(r0)) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], + Memr[AGR_RESID(agr)], Memr[AGR_RESSQ(agr)], naperts))) + else + call pargr (INDEFR) + if (fitok == YES && ! IS_INDEFR(r0)) { + call sprintf (Memc[title+strlen(Memc[title])], + 2 * SZ_LINE - strlen (Memc[title]), + "\nDashed line: theoretical cog ") + call sprintf (Memc[title+strlen(Memc[title])], + 2 * SZ_LINE - strlen (Memc[title]), "Solid line: adopted cog") + } + call glabax (gd, Memc[title], "Aperture Radius", "Curve of Growth") + + # Draw the fitted and rejected data. + do i = 1, npts { + call gpmark (gd, Memr[AGR_RBAR(agr)+1], mag[2,i], inap[i] - 1, + GM_PLUS, 2.0, 2.0) + if (nap[i] > inap[i]) + call gpmark (gd, Memr[AGR_RBAR(agr)+inap[i]], + mag[inap[i]+1,i], nap[i] - inap[i], GM_CROSS, 2.0, 2.0) + } + + # Plot the theoretical and adopt models. + if (fitok == YES && ! IS_INDEFR(r0)) { + call gseti (gd, G_PLTYPE, GL_DASHED) + call gpline (gd, Memr[AGR_RBAR(agr)+1], Memr[AGR_THEO(agr)+1], + naperts - 1) + call gseti (gd, G_PLTYPE, GL_SOLID) + call gpline (gd, Memr[AGR_RBAR(agr)+1], Memr[AGR_ADOPT(agr)+1], + naperts - 1) + } + + call gflush (gd) + call sfree (sp) +end + + +# PH_GAIMRES -- Graph the residuals from the adopted model for a particular +# image. + +procedure ph_gaimres (gd, agr, image, r0, xairmass, inap, nap, rap, mag, + naperts, npts, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good aperture numbers +int nap[ARB] # the array of aperture numbers +real rap[naperts,ARB] # input array of aperture radii +real mag[naperts,ARB] # input array of magnitudes / differences +int naperts # the number of input apertures +int npts # the number of points +int fitok # is the fit ok + +int i +pointer sp, title, diff +real rmin, rmax, dr, dmin, dmax, ddmin, ddmax, dd +int strlen() +real ph_achi() + +begin + call smark (sp) + call salloc (title, 2 * SZ_LINE, TY_CHAR) + call salloc (diff, naperts, TY_REAL) + + # Set the data window. + call alimr (rap, naperts * npts, rmin, rmax) + dr = rmax - rmin + rmin = rmin - 0.1 * dr + rmax = rmax + 0.1 * dr + + dmin = MAX_REAL + dmax = -MAX_REAL + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + call alimr (Memr[diff+1], nap[i] - 1, ddmin, ddmax) + if (ddmin < dmin) + dmin = ddmin + if (ddmax > dmax) + dmax = ddmax + } + if (dmin > dmax) { + dmin = -0.1 + dmax = 0.1 + } else { + dd = dmax - dmin + dmin = dmin - 0.1 * dd + dmax = dmax + 0.1 * dd + } + + call gclear (gd) + call gswind (gd, rmin, rmax, dmin, dmax) + + # Set up the axes and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen(Memc[title]), "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g\n") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], Memr[AGR_RESID(agr)], + Memr[AGR_RESSQ(agr)], naperts))) + call glabax (gd, Memc[title], "Aperture Radius", "Residuals") + + # Draw the data. + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + call gpmark (gd, Memr[AGR_RBAR(agr)+1], Memr[diff+1], inap[i] - 1, + GM_PLUS, 2.0, 2.0) + if (nap[i] > inap[i]) + call gpmark (gd, Memr[AGR_RBAR(agr)+inap[i]], + Memr[diff+inap[i]], nap[i] - inap[i], GM_CROSS, 2.0, 2.0) + } + + call gseti (gd, G_PLTYPE, GL_SOLID) + call gamove (gd, rmin, 0.0) + call gadraw (gd, rmax, 0.0) + call gpline (gd, Memr[AGR_RBAR(agr)+1], Memr[AGR_AVE(agr)+1], + naperts - 1) + + call gflush (gd) + call sfree (sp) +end + + + +# PH_GBIMRES -- Graph the residuals from the adopted model for a particular +# image as a function of magnitude in the first aperture. + +procedure ph_gbimres (gd, agr, image, r0, xairmass, inap, nap, mag, + naperts, npts, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good aperture numbers +int nap[ARB] # the array of aperture numbers +real mag[naperts,ARB] # input array of magnitudes / differences +int naperts # the number of input apertures +int npts # the number of points +int fitok # is the fit ok + +int i, j +pointer sp, title, diff +real rmin, rmax, dr, dmin, dmax, ddmin, ddmax, dd +int strlen() +real ph_achi() + +begin + call smark (sp) + call salloc (title, 2 * SZ_LINE, TY_CHAR) + call salloc (diff, naperts, TY_REAL) + + # Set the data window. + rmin = MAX_REAL + rmax = -MAX_REAL + do i = 1, npts { + if (mag[1,i] < rmin) + rmin = mag[1,i] + if (mag[1,i] > rmax) + rmax = mag[1,i] + } + if (rmin > rmax) { + rmin = -0.1 + rmax = 0.1 + } else { + dr = rmax - rmin + rmin = rmin - 0.1 * dr + rmax = rmax + 0.1 * dr + } + + dmin = MAX_REAL + dmax = -MAX_REAL + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + call alimr (Memr[diff+1], nap[i] - 1, ddmin, ddmax) + if (ddmin < dmin) + dmin = ddmin + if (ddmax > dmax) + dmax = ddmax + } + if (dmin > dmax) { + dmin = -0.1 + dmax = 0.1 + } else { + dd = dmax - dmin + dmin = dmin - 0.1 * dd + dmax = dmax + 0.1 * dd + } + + call gclear (gd) + call gswind (gd, rmin, rmax, dmin, dmax) + + # Set up the axes and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen(Memc[title]), "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g\n") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], Memr[AGR_RESID(agr)], + Memr[AGR_RESSQ(agr)], naperts))) + call glabax (gd, Memc[title], "Magnitude", "Residuals") + + # Draw the data. + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + do j = 2, inap[i] + call gmark (gd, mag[1,i], Memr[diff+j-1], GM_PLUS, 2.0, 2.0) + do j = inap[i] + 1, nap[i] + call gmark (gd, mag[1,i], Memr[diff+j-1], GM_CROSS, 2.0, 2.0) + } + + call gseti (gd, G_PLTYPE, GL_SOLID) + call gamove (gd, rmin, 0.0) + call gadraw (gd, rmax, 0.0) + + call gflush (gd) + call sfree (sp) +end + + +# PH_GAXIMRES -- Graph the residuals from the adopted model for a particular +# image. + +procedure ph_gaximres (gd, agr, image, r0, xairmass, inap, nap, x, mag, + naperts, npts, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good aperture numbers +int nap[ARB] # the array of aperture numbers +real x[ARB] # the array of x values +real mag[naperts,ARB] # input array of magnitudes / differences +int naperts # the number of input apertures +int npts # the number of points +int fitok # is the fit ok + +int i, j +pointer sp, title, diff +real xmin, xmax, dx, dmin, dmax, ddmin, ddmax, dd +int strlen() +real ph_achi() + +begin + call smark (sp) + call salloc (title, 2 * SZ_LINE, TY_CHAR) + call salloc (diff, naperts, TY_REAL) + + # Set the x data window. + call alimr (x, npts, xmin, xmax) + dx = xmax - xmin + xmin = xmin - 0.1 * dx + xmax = xmax + 0.1 * dx + + # Set the y data window. + dmin = MAX_REAL + dmax = -MAX_REAL + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + call alimr (Memr[diff+1], nap[i] - 1, ddmin, ddmax) + if (ddmin < dmin) + dmin = ddmin + if (ddmax > dmax) + dmax = ddmax + } + if (dmin > dmax) { + dmin = -0.1 + dmax = 0.1 + } else { + dd = dmax - dmin + dmin = dmin - 0.1 * dd + dmax = dmax + 0.1 * dd + } + + # Initialize the plot. + call gclear (gd) + call gswind (gd, xmin, xmax, dmin, dmax) + + # Set up the axes and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen(Memc[title]), "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g\n") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], Memr[AGR_RESID(agr)], + Memr[AGR_RESSQ(agr)], naperts))) + call glabax (gd, Memc[title], "X Coordinate", "Residuals") + + # Draw the data. + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + do j = 2, inap[i] + call gmark (gd, x[i], Memr[diff+j-1], GM_PLUS, 2.0, 2.0) + do j = inap[i] + 1, nap[i] + call gmark (gd, x[i], Memr[diff+j-1], GM_CROSS, 2.0, 2.0) + } + + call gseti (gd, G_PLTYPE, GL_SOLID) + call gamove (gd, xmin, 0.0) + call gadraw (gd, xmax, 0.0) + + call gflush (gd) + call sfree (sp) +end + + +# PH_GAYIMRES -- Graph the residuals from the adopted model for a particular +# image. + +procedure ph_gayimres (gd, agr, image, r0, xairmass, inap, nap, y, mag, + naperts, npts, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good aperture numbers +int nap[ARB] # the array of aperture numbers +real y[ARB] # the array of y values +real mag[naperts,ARB] # input array of magnitudes / differences +int naperts # the number of input apertures +int npts # the number of points +int fitok # is the fit ok + +int i, j +pointer sp, title, diff +real ymin, ymax, dy, dmin, dmax, ddmin, ddmax, dd +int strlen() +real ph_achi() + +begin + call smark (sp) + call salloc (title, 2 * SZ_LINE, TY_CHAR) + call salloc (diff, naperts, TY_REAL) + + # Set the y data window. + call alimr (y, npts, ymin, ymax) + dy = ymax - ymin + ymin = ymin - 0.1 * dy + ymax = ymax + 0.1 * dy + + # Set the y data window. + dmin = MAX_REAL + dmax = -MAX_REAL + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + call alimr (Memr[diff+1], nap[i] - 1, ddmin, ddmax) + if (ddmin < dmin) + dmin = ddmin + if (ddmax > dmax) + dmax = ddmax + } + if (dmin > dmax) { + dmin = -0.1 + dmax = 0.1 + } else { + dd = dmax - dmin + dmin = dmin - 0.1 * dd + dmax = dmax + 0.1 * dd + } + + # Initialize the plot. + call gclear (gd) + call gswind (gd, ymin, ymax, dmin, dmax) + + # Set up the axes and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen(Memc[title]), "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g\n") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], Memr[AGR_RESID(agr)], + Memr[AGR_RESSQ(agr)], naperts))) + call glabax (gd, Memc[title], "Y Coordinate", "Residuals") + + # Draw the data. + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[diff+1], + nap[i] - 1) + do j = 2, inap[i] + call gmark (gd, y[i], Memr[diff+j-1], GM_PLUS, 2.0, 2.0) + do j = inap[i] + 1, nap[i] + call gmark (gd, y[i], Memr[diff+j-1], GM_CROSS, 2.0, 2.0) + } + + call gseti (gd, G_PLTYPE, GL_SOLID) + call gamove (gd, ymin, 0.0) + call gadraw (gd, ymax, 0.0) + + call gflush (gd) + call sfree (sp) +end + + +# PH_GACUM -- Plot the cumulative profile. + +procedure ph_gacum (gd, agr, image, r0, xairmass, inap, nap, rap, mag, naperts, + npts, smallap, largeap, fitok) + +pointer gd # pointer to the graphics descriptor +pointer agr # pointer to the fit structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass of the image +int inap[ARB] # the array of good number apertures +int nap[ARB] # the array of number of apertures +real rap[ARB] # input array of aperture radii +real mag[naperts,ARB] # the array of magnitudes +int naperts # the number of input apertures +int npts # the number of points +int smallap # small aperture number +int largeap # large aperture number +int fitok # is the fit ok + +int i, j +pointer sp, title, cmags, ctheo, cadopt, caderr, cptr +real rmin, rmax, dr, dmin, dmax, ddmin, ddmax, dd +int strlen() +real ph_achi() + +begin + # Allocate working space. + call smark (sp) + call salloc (title, SZ_LINE, TY_CHAR) + call salloc (cmags, (naperts + 1) * npts, TY_REAL) + call salloc (ctheo, naperts + 1, TY_REAL) + call salloc (cadopt, naperts + 1, TY_REAL) + call salloc (caderr, naperts + 1, TY_REAL) + + # Compute the data. + call aclrr (Memr[cmags], (naperts + 1) * npts) + cptr = cmags + do j = 1, npts { + do i = largeap, 2, -1 { + if (i > nap[j]) + Memr[cptr+i-1] = Memr[AGR_ADOPT(agr)+i-1] + Memr[cptr+i] + else + Memr[cptr+i-1] = mag[i,j] + Memr[cptr+i] + } + cptr = cptr + naperts + 1 + } + + # Compute the model. + call aclrr (Memr[ctheo], naperts + 1) + call aclrr (Memr[cadopt], naperts + 1) + call aclrr (Memr[caderr], naperts + 1) + do i = largeap, 2, -1 { + Memr[ctheo+i-1] = Memr[AGR_THEO(agr)+i-1] + Memr[ctheo+i] + Memr[cadopt+i-1] = Memr[AGR_ADOPT(agr)+i-1] + Memr[cadopt+i] + Memr[caderr+i-1] = Memr[AGR_WADO(agr)+i-1] ** 2 + Memr[caderr+i] + } + + # Set the x data window. + call alimr (rap, naperts, rmin, rmax) + dr = rmax - rmin + rmin = rmin - 0.1 * dr + rmax = rmax + 0.1 * dr + + # Set the y data window. + dmin = MAX_REAL + dmax = -MAX_REAL + call alimr (Memr[cmags], (naperts + 1) * npts, ddmin, ddmax) + if (ddmin < dmin) + dmin = ddmin + if (ddmax > dmax) + dmax = ddmax + dd = dmax - dmin + dmin = dmin - 0.1 * dd + dmax = dmax + 0.1 * dd + + call gclear (gd) + call gswind (gd, rmin, rmax, dmax, dmin) + + # Set up the title and the axis labels. + call sysid (Memc[title], 2 * SZ_LINE) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen(Memc[title]), "\nImage: %s Ro: %.3f X: %.3f Rms: %.3g") + call pargstr (image) + call pargr (r0) + call pargr (xairmass) + call pargr (sqrt (ph_achi (Memr[AGR_WR(agr)], Memr[AGR_RESID(agr)], + Memr[AGR_RESSQ(agr)], naperts))) + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen (Memc[title]), "\nDashed line: theoretical apcor ") + call sprintf (Memc[title+strlen(Memc[title])], 2 * SZ_LINE - + strlen (Memc[title]), "Solid line: adopted apcor") + call glabax (gd, Memc[title], "Aperture Radius", "Aperture Correction") + + # Draw the data. + cptr = cmags + do j = 1, npts { + if (largeap < inap[j]) + call gpmark (gd, rap, Memr[cptr+1], largeap - 1, GM_PLUS, + 2.0, 2.0) + else + call gpmark (gd, rap, Memr[cptr+1], inap[j] - 1, GM_PLUS, + 2.0, 2.0) + if (largeap > inap[j]) + call gpmark (gd, rap[inap[j]], Memr[cptr+inap[j]], + largeap - inap[j], GM_CROSS, 2.0, 2.0) + cptr = cptr + naperts + 1 + } + + # Draw the zero correction line. + call gamove (gd, rmin, 0.0) + call gadraw (gd, rmax, 0.0) + + # Draw the aperture limits. + if (smallap > rap[1]) { + call gseti (gd, G_PLTYPE, GL_DASHED) + call gamove (gd, rap[1], dmin) + call gadraw (gd, rap[1], dmax) + call gseti (gd, G_PLTYPE, GL_SOLID) + } + call gamove (gd, rap[smallap], dmin) + call gadraw (gd, rap[smallap], dmax) + call gamove (gd, rap[largeap], dmin) + call gadraw (gd, rap[largeap], dmax) + if (rap[naperts] > largeap) { + call gseti (gd, G_PLTYPE, GL_DASHED) + call gamove (gd, rap[naperts], dmin) + call gadraw (gd, rap[naperts], dmax) + call gseti (gd, G_PLTYPE, GL_SOLID) + } + + # Draw the model. + call gseti (gd, G_PLTYPE, GL_DASHED) + call gpline (gd, rap, Memr[ctheo+1], naperts) + call gseti (gd, G_PLTYPE, GL_SOLID) + call gpline (gd, rap, Memr[cadopt+1], naperts) + + call gflush (gd) + call sfree (sp) +end + + +# PH_PNEAREST -- Print basic information for the nearest object + +procedure ph_pnearest (gd, graphtype, wx, wy, agr, id, x, y, inap, nap, rap, + mag, naperts, npts, smallap, largeap) + +pointer gd # pointer to the graphics stream +int graphtype # the current graphtype +real wx, wy # the coordinates of the point to be un/deleted +pointer agr # pointer to the fitting structure +int id[ARB] # array of star ids +real x[ARB] # the x coordinates +real y[ARB] # the y coordinates +int inap[ARB] # the number of good apertures +int nap[ARB] # the number of measured apertures +real rap[naperts,ARB] # the list of aperture radii +real mag[naperts,ARB] # the list of magnitude differences +int naperts # the number of apertures +int npts # the number of points +int smallap # the small aperture number +int largeap # the large aperture number + +int i, j +pointer sp, xin, yin, iap, inpts, index +real r2min, r2, xc, yc +real ph_nearest() + +begin + # Allocate working space. + call smark (sp) + call salloc (xin, naperts + 1, TY_REAL) + call salloc (yin, naperts + 1, TY_REAL) + + # Initialize. + r2min = MAX_REAL + iap = 0 + inpts = 0 + + # Find the nearest point. + switch (graphtype) { + case AGR_FIT: + + do i = 1, npts { + r2 = ph_nearest (gd, wx, wy, Memr[AGR_RBAR(agr)+1], mag[2,i], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = Memr[AGR_RBAR(agr)+iap-1] + yc = mag[iap,inpts] + } + } + + case AGR_ARESIDUALS: + + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[AGR_RBAR(agr)+1], Memr[yin+1], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = Memr[AGR_RBAR(agr)+iap-1] + yc = mag[iap,inpts] - Memr[AGR_ADOPT(agr)+iap-1] + } + } + + case AGR_BRESIDUALS: + do i = 1, npts { + call amovkr (mag[1,i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = mag[1,inpts] + yc = mag[iap,inpts] - Memr[AGR_ADOPT(agr)+iap-1] + } + } + + case AGR_XRESIDUALS: + + do i = 1, npts { + call amovkr (x[i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = x[inpts] + yc = mag[iap,inpts] - Memr[AGR_ADOPT(agr)+iap-1] + } + } + + case AGR_YRESIDUALS: + + do i = 1, npts { + call amovkr (y[i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = y[inpts] + yc = mag[iap,inpts] - Memr[AGR_ADOPT(agr)+iap-1] + } + } + + case AGR_CUMULATIVE: + + do i = 1, npts { + call aclrr (Memr[yin], naperts + 1) + do j = largeap, 2, -1 { + if (j > nap[i]) + Memr[yin+j-1] = Memr[AGR_ADOPT(agr)+j-1] + Memr[yin+j] + else + Memr[yin+j-1] = mag[j,i] + Memr[yin+j] + } + r2 = ph_nearest (gd, wx, wy, rap[1,i], Memr[yin+1], + nap[i] - 1, nap[i] - 1, index, YES) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + xc = rap[iap,inpts] + yc = Memr[yin+iap-1] + } + } + } + + if (iap != 0 && inpts != 0) { + switch (graphtype) { + case AGR_FIT: + call printf ( + "Star: %d x: %.3f y: %.3f radius: %.3f delta mag: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + case AGR_ARESIDUALS: + call printf ( + "Star: %d x: %.3f y: %.3f radius: %.3f residual: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + case AGR_BRESIDUALS: + call printf ( + "Star: %d x: %.3f y: %.3f mag[1]: %.3f residual: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + case AGR_XRESIDUALS: + call printf ( + "Star: %d x: %.3f y: %.3f x: %.3f residual: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + case AGR_YRESIDUALS: + call printf ( + "Star: %d x: %.3f y: %.3f y: %.3f residual: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + case AGR_CUMULATIVE: + call printf ( + "Star: %d x: %.3f y: %.3f y: %.3f apercor: %.3f\n") + call pargi (id[inpts]) + call pargr (x[inpts]) + call pargr (y[inpts]) + call pargr (xc) + call pargr (yc) + } + } + + call sfree (sp) +end + + +# PH_AUDELETE -- Delete or undelete points from the current graph. + +int procedure ph_audelete (gd, graphtype, wx, wy, agr, x, y, inap, nap, rap, + mag, naperts, npts, smallap, largeap, delete) + +pointer gd # pointer to the graphics stream +int graphtype # the current graphtype +real wx, wy # the coordinates of the point to be un/deleted +pointer agr # pointer to the fitting structure +real x[ARB] # the x coordinates +real y[ARB] # the y coordinates +int inap[ARB] # the number of good apertures +int nap[ARB] # the number of measured apertures +real rap[naperts,ARB] # the list of aperture radii +real mag[naperts,ARB] # the list of magnitude differences +int naperts # the number of apertures +int npts # the number of points +int smallap # the small aperture number +int largeap # the large aperture number +int delete # delete points ?, otherwise undelete + +int i, j, iap, inpts, index, stat +pointer sp, xin, yin +real r2min, r2 +real ph_nearest() + +begin + # Allocate working space. + call smark (sp) + call salloc (xin, naperts + 1, TY_REAL) + call salloc (yin, naperts + 1, TY_REAL) + + # Initialize. + r2min = MAX_REAL + iap = 0 + inpts = 0 + + # Find the nearest point. + switch (graphtype) { + case AGR_FIT: + + do i = 1, npts { + r2 = ph_nearest (gd, wx, wy, Memr[AGR_RBAR(agr)+1], mag[2,i], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + + case AGR_ARESIDUALS: + + do i = 1, npts { + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[AGR_RBAR(agr)+1], Memr[yin+1], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + + case AGR_BRESIDUALS: + do i = 1, npts { + call amovkr (mag[1,i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + + case AGR_XRESIDUALS: + + do i = 1, npts { + call amovkr (x[i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + + case AGR_YRESIDUALS: + + do i = 1, npts { + call amovkr (y[i], Memr[xin+1], nap[i] - 1) + call asubr (mag[2,i], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[i] - 1) + r2 = ph_nearest (gd, wx, wy, Memr[xin+1], Memr[yin+1], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + + case AGR_CUMULATIVE: + + do i = 1, npts { + call aclrr (Memr[yin], naperts + 1) + do j = largeap, 2, -1 { + if (j > nap[i]) + Memr[yin+j-1] = Memr[AGR_ADOPT(agr)+j-1] + Memr[yin+j] + else + Memr[yin+j-1] = mag[j,i] + Memr[yin+j] + } + r2 = ph_nearest (gd, wx, wy, rap[1,i], Memr[yin+1], + inap[i] - 1, nap[i] - 1, index, delete) + if (r2 < r2min) { + r2min = r2 + iap = index + 1 + inpts = i + } + } + } + + # Delete or undelete the point and mark it. + if (iap != 0 && inpts != 0) { + + switch (graphtype) { + case AGR_FIT: + call amovr (Memr[AGR_RBAR(agr)+1], Memr[xin+1], nap[inpts] - 1) + call amovr (mag[2,inpts], Memr[yin+1], nap[inpts] - 1) + case AGR_ARESIDUALS: + call amovr (Memr[AGR_RBAR(agr)+1], Memr[xin+1], nap[inpts] - 1) + call asubr (mag[2,inpts], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[inpts] - 1) + case AGR_BRESIDUALS: + call amovkr (mag[1,inpts], Memr[xin+1], nap[inpts] - 1) + call asubr (mag[2,inpts], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[inpts] - 1) + case AGR_XRESIDUALS: + call amovkr (x[inpts], Memr[xin+1], nap[inpts] - 1) + call asubr (mag[2,inpts], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[inpts] - 1) + case AGR_YRESIDUALS: + call amovkr (y[inpts], Memr[xin+1], nap[inpts] - 1) + call asubr (mag[2,inpts], Memr[AGR_ADOPT(agr)+1], Memr[yin+1], + nap[inpts] - 1) + case AGR_CUMULATIVE: + call amovr (rap[1,inpts], Memr[xin+1], nap[inpts] - 1) + do j = largeap, 2, -1 { + if (j > nap[inpts]) + Memr[yin+j-1] = Memr[AGR_ADOPT(agr)+j-1] + Memr[yin+j] + else + Memr[yin+j-1] = mag[j,inpts] + Memr[yin+j] + } + } + + if (delete == YES) { + inap[inpts] = iap - 1 + do i = inap[inpts] + 1, nap[inpts] { + call gscur (gd, Memr[xin+i-1], Memr[yin+i-1]) + call gmark (gd, Memr[xin+i-1], Memr[yin+i-1], GM_CROSS, + 2.0, 2.0) + } + } else { + do i = inap[inpts]+1, iap { + call gscur (gd, Memr[xin+i-1], Memr[yin+i-1]) + call gseti (gd, G_PMLTYPE, GL_CLEAR) + call gmark (gd, Memr[xin+i-1], Memr[yin+i-1], GM_CROSS, + 2.0, 2.0) + call gseti (gd, G_PMLTYPE, GL_SOLID) + call gmark (gd, Memr[xin+i-1], Memr[yin+i-1], GM_PLUS, + 2.0, 2.0) + } + inap[inpts] = iap + } + + stat = YES + + } else + stat = NO + + call sfree (sp) + return (stat) +end + + +# PH_AEPRINT -- Print the appropriate error message under the plot. + +procedure ph_aeprint (agr, image, r0, rap, naperts, smallap, largeap, fitok, + newfit) + +pointer agr # pointer to the fitting structure +char image[ARB] # the image name +real r0 # the seing rdius +real rap[ARB] # the lists of aperture radii +int naperts # the number of aperture radii +int smallap # the index of the small aperture radius +int largeap # the index of the large aperture radius +int fitok # did the fit converge ? +int newfit # is the fit out-of-date ? + +begin + if (fitok == NO) { + call printf ("Error: The cog model fit did not converge\n") + } else if (newfit == YES) { + call printf ("Warning: The cog model fit is out-of-date\n") + } else if (IS_INDEFR(r0)) { + call printf ("Warning: Unable to fit RO for image %s\n") + call pargstr (image) + } else { + call ph_papcor (STDOUT, image, r0, rap, Memr[AGR_ADOPT(agr)], + Memr[AGR_WADO(agr)], naperts, smallap, largeap) + } +end + + +# PH_NEAREST -- Find the nearest data point to the coordinates. + +real procedure ph_nearest (gd, wx, wy, x, y, inpts, npts, index, delete) + +pointer gd # the graphics descriptor +real wx, wy # the cursor coordinates +real x[ARB] # the x values +real y[ARB] # the y values +int inpts # the number of good points +int npts # the number of points +int index # the index of the nearest point +int delete # delete or undelete points + +int i, ib, ie +real r2min, r2, nwx, nwy, nx, ny + +begin + r2min = MAX_REAL + index = 0 + if (delete == YES) { + ib = 1 + ie = inpts + } else { + ib = inpts + 1 + ie = npts + } + call gctran (gd, wx, wy, nwx, nwy, 1, 0) + + do i = ib, ie { + call gctran (gd, x[i], y[i], nx, ny, 1, 0) + r2 = (nx - nwx) ** 2 + (ny - nwy) ** 2 + if (r2 < r2min) { + r2min = r2 + index = i + } + } + + return (r2min) +end + + +# PH_ACOLON -- Execute the apfile colon commands. + +procedure ph_acolon (gd, imtable, agr, symbol, isymbol, params, lterms, + naperts, smallap, largeap, cmd, fitok, newfit, newimage, newgraph) + +pointer gd # pointer to the graphics stream +pointer imtable # pointer to the image name symbol table +pointer agr # pointer to the fitting structure +pointer symbol # the current image symbol +int isymbol # the index of the current image symbol +double params[ARB] # the initial parameters +int lterms # the number of parameters to fit +int naperts # the number of apertures +int smallap # the small aperture number +int largeap # the large aperture number +char cmd[ARB] # the colon command +int fitok # is the fit ok +int newfit # compute a new fit +int newimage # evaluate the fit for a new image +int newgraph # plot a new graph + +int ncmd, ival +pointer sp, incmd, outcmd, newsymbol +real rval +int strdic(), nscan(), stnsymbols() +pointer stfind(), sthead(), stnext() + +begin + # Allocate working space. + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmd) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + ncmd = strdic (Memc[incmd], Memc[incmd], SZ_LINE, AGR_CMDS) + + switch (ncmd) { + case AGR_CMD_SHOW: + + call gargwrd (Memc[outcmd], SZ_LINE) + ncmd = strdic (Memc[outcmd], Memc[outcmd], SZ_LINE, AGR_SHOWCMDS) + switch (ncmd) { + case AGR_CMD_MODEL: + if (fitok == NO) { + call printf ("Error: The cog model fit did not converge\n") + } else { + if (newfit == YES) + call printf ( + "Warning: The cog model fit is out-of-date\n") + call gdeactivate (gd, 0) + call ph_pshow (STDOUT, agr, lterms) + call greactivate (gd, 0) + } + case AGR_CMD_SEEING: + if (fitok == NO) { + call printf ("Error: The cog model fit did not converge\n") + } else { + if (newfit == YES) + call printf ( + "Warning: The cog model fit is out-of-date\n") + call gdeactivate (gd, 0) + call ph_rshow (STDOUT, imtable) + call greactivate (gd, 0) + } + case AGR_CMD_PARAMETERS: + call gdeactivate (gd, 0) + call ph_ishow (STDOUT, params, lterms) + call greactivate (gd, 0) + default: + if (fitok == NO) { + call printf ("Error: The cog model fit did not converge\n") + } else { + if (newfit == YES) + call printf ( + "Warning: The cog model fit is out-of-date\n") + call gdeactivate (gd, 0) + call ph_pshow (STDOUT, agr, lterms) + call greactivate (gd, 0) + } + } + + case AGR_CMD_IMAGE: + call gargwrd (Memc[outcmd], SZ_LINE) + if (nscan() == 1) { + call printf ("Image: %s R0: %7.3f X: %5.3f\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (IMT_RO(symbol)) + call pargr (IMT_XAIRMASS(symbol)) + } else { + newsymbol = stfind (imtable, Memc[outcmd]) + if (newsymbol == NULL) { + call printf ("Image: %s not found in data\n") + call pargstr (Memc[outcmd]) + } else { + symbol = newsymbol + isymbol = 1 + newsymbol = sthead (imtable) + while (newsymbol != NULL) { + if (newsymbol == symbol) + break + isymbol = isymbol + 1 + newsymbol = stnext (imtable, newsymbol) + } + isymbol = stnsymbols(imtable, 0) - isymbol + 1 + newimage = YES + newgraph = YES + } + } + + case AGR_CMD_MTERMS: + call gargi (ival) + if (nscan() == 1) { + call printf ("nparams = %d\n") + call pargi (lterms) + } else { + lterms = max (1, min (ival, MAX_MTERMS)) + newfit = YES + } + + case AGR_CMD_SWINGS: + call gargr (rval) + if (nscan() == 1) { + call printf ("swings = %g\n") + call pargr (real(params[1])) + } else { + params[1] = max (1.0, rval) + } + + case AGR_CMD_PWINGS: + call gargr (rval) + if (nscan() == 1) { + call printf ("pwings = %g\n") + call pargr (real(params[2])) + } else { + params[2] = max (0.0, min (1.0, rval)) + } + + case AGR_CMD_PGAUSS: + call gargr (rval) + if (nscan() == 1) { + call printf ("pgauss = %g\n") + call pargr (real(params[3])) + } else { + params[3] = max (0.0, min (1.0, rval)) + } + + case AGR_CMD_RGESCALE: + call gargr (rval) + if (nscan() == 1) { + call printf ("rgescale = %g\n") + call pargr (real(params[4])) + } else { + params[4] = max (0.0, rval) + } + + case AGR_CMD_XWINGS: + call gargr (rval) + if (nscan() == 1) { + call printf ("xwings = %g\n") + call pargr (real(params[5])) + } else { + params[5] = rval + } + + case AGR_CMD_SMALLAP: + call gargi (ival) + if (nscan() == 1) { + call printf ("smallap = %d\n") + call pargi (smallap) + } else { + smallap = max (1, min (ival, naperts)) + } + + case AGR_CMD_LARGEAP: + call gargi (ival) + if (nscan() == 1) { + call printf ("largeap = %d\n") + call pargi (largeap) + } else { + largeap = max (1, min (ival, naperts)) + } + + default: + call printf ("Ambiguous or undefined command\7\n") + } + + + call sfree (sp) +end |