diff options
Diffstat (limited to 'noao/digiphot/photcal/mkobsfile')
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/apfile.key | 36 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/dinvers.f | 53 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/mkpkg | 19 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phagrow.x | 1982 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phaigrow.x | 1635 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phaimtable.x | 409 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phimtable.x | 971 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phmatch.x | 487 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/phsort.x | 528 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/t_apfile.x | 727 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/t_mkphotcors.x | 273 | ||||
| -rw-r--r-- | noao/digiphot/photcal/mkobsfile/t_obsfile.x | 1131 |
12 files changed, 8251 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/mkobsfile/apfile.key b/noao/digiphot/photcal/mkobsfile/apfile.key new file mode 100644 index 00000000..b6cdd1f6 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/apfile.key @@ -0,0 +1,36 @@ + Keystroke Commands + +? Print help +w Print computed aperture correction +c Print coordinates of star nearest cursor +f Compute a new fit +d Delete point(s) nearest the cursor +u Undelete point(s) nearest the cursor +m Plot the observed and model cog versus radius +r Plot the cog residuals versus radius +b Plot the cog residuals versus magnitude +x Plot the cog residuals versus the x coordinate +y Plot the cog residuals versus the y coordinate +a Plot the aperture correction versus radius +g Redraw the current plot +n Move to the next image +p Move to the previous image +q Quit task + + Colon commands + +:show parameters Show the initial cog model parameter values +:show model Show the fitted cog model parameters +:show seeing Show the computed seeing radii for all images +:image [value] Show/set the image to be analyzed + + Colon Parameter Editing Commands + +:nparams [value] Show/set the number of parameters in model to fit +:swings [value] Show/set initial power law slope of stellar wings +:pwings [value] Show/set fraction of total power in stellar wings +:pgauss [value] Show/set fraction of total core power in gaussian +:rgescale [value] Show/set ratio of exponential to gaussian radial scales +:xwings [values] Show/set the extinction coefficient +:smallap [value] The index of the smallest aperture +:largeap [value] The index of the largest aperture diff --git a/noao/digiphot/photcal/mkobsfile/dinvers.f b/noao/digiphot/photcal/mkobsfile/dinvers.f new file mode 100644 index 00000000..5b1d8a5f --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/dinvers.f @@ -0,0 +1,53 @@ +c +c subroutine dinver (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 dinver (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.0d0) 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/photcal/mkobsfile/mkpkg b/noao/digiphot/photcal/mkobsfile/mkpkg new file mode 100644 index 00000000..6ba9935f --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/mkpkg @@ -0,0 +1,19 @@ +# The MKPKG file for the mkobsfile subdirectory. + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + t_obsfile.x "../lib/obsfile.h" <fset.h> <ctotok.h> <ctype.h> + t_mkphotcors.x "../lib/obsfile.h" + phimtable.x "../lib/obsfile.h" <fset.h> + phsort.x "../lib/obsfile.h" + phmatch.x <mach.h> "../lib/obsfile.h" + t_apfile.x "../lib/apfile.h" <fset.h> <ctotok.h> <ctype.h> + phaimtable.x "../lib/apfile.h" <fset.h> + phagrow.x "../lib/apfile.h" <time.h> <mach.h> <fset.h> + phaigrow.x "../lib/apfile.h" <time.h> <mach.h> <gset.h> + dinvers.f + ; diff --git a/noao/digiphot/photcal/mkobsfile/phagrow.x b/noao/digiphot/photcal/mkobsfile/phagrow.x new file mode 100644 index 00000000..a8c6dfce --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phagrow.x @@ -0,0 +1,1982 @@ +include <time.h> +include <mach.h> +include <fset.h> +include "../lib/apfile.h" + + +# PH_AGROW - Compute the curve of growth using the data from the input +# photometry files, optional airmasses from the airmass file, and the +# initial values of the parameters. + +procedure ph_agrow (apcat, magfd, logfd, mgd, imtable, id, x, y, nap, rap, + mag, merr, naperts, params, lterms, smallap, largeap) + +int apcat # the aperture correction file descriptor +int magfd # the output magnitude file descriptor +int logfd # the log file descriptor +pointer mgd # pointer to the plot metacode file +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 of the parameters +int lterms # the number of terms to be fit +int smallap # the small aperture number +int largeap # the large aperture number + +int nimages +pointer agr +int stnsymbols() +int ph_amfit() + +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, "There is no input data to fit\n") + return + } + + # Allocate memory required for fitting. + call ph_ameminit (agr, lterms, naperts) + + # Fit the seeing radii and the cog model parameters. + if (ph_amfit (imtable, agr, nap, rap, mag, merr, naperts, params, + lterms) == ERR) { + 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 (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, + nap, rap, mag, merr, naperts, smallap, largeap) + if (logfd != NULL) + call ph_tfshow (logfd, agr, naperts) + } + + # Free memory required for fitting. + call ph_amemfree (agr) +end + + +# PH_LOGTITLE -- Print the title for the new log file entry + +procedure ph_logtitle (logfd, apcat) + +int logfd # the log file descriptor +int apcat # the aperture correction file descriptor + +pointer sp, afname, date +long clktime() + +begin + # Allocate working space. + call smark (sp) + call salloc (afname, SZ_FNAME, TY_CHAR) + call salloc (date, SZ_TIME, TY_CHAR) + + # Get the file name and the times. + call fstats (apcat, F_FILENAME, Memc[afname], SZ_FNAME) + call cnvtime (clktime (long(0)), Memc[date], SZ_TIME) + + call fprintf (logfd, "NEW LOGFILE ENTRY AT: %s\n") + call pargstr (Memc[date]) + call fprintf (logfd, "APFILE: %s\n\n") + call pargstr (Memc[afname]) + + call sfree (sp) +end + + +# PH_AMEMINIT -- Allocate memory for doing the fit. + +procedure ph_ameminit (agr, lterms, naperts) + +pointer agr # pointer to the fitting structure +int lterms # the number of terms to fit +int naperts # the number of apertures + +begin + call malloc (agr, LEN_AGRSTRUCT, TY_STRUCT) + + call malloc (AGR_DM(agr), naperts, TY_DOUBLE) + call malloc (AGR_DDDR(agr), naperts, TY_DOUBLE) + call malloc (AGR_T(agr), lterms * naperts, TY_DOUBLE) + call malloc (AGR_U(agr), lterms * lterms, TY_DOUBLE) + call malloc (AGR_V(agr), lterms, TY_DOUBLE) + call malloc (AGR_POLD(agr), lterms, TY_DOUBLE) + call malloc (AGR_DP(agr), lterms, TY_DOUBLE) + call malloc (AGR_PARAMS(agr), MAX_MTERMS, TY_DOUBLE) + call malloc (AGR_PERRORS(agr), MAX_MTERMS, TY_DOUBLE) + call malloc (AGR_PCLAMPS(agr), MAX_MTERMS, TY_DOUBLE) + + call malloc (AGR_RBAR(agr), naperts, TY_REAL) + call malloc (AGR_W(agr), naperts, TY_REAL) + call malloc (AGR_THEO(agr), naperts, TY_REAL) + call malloc (AGR_ADOPT(agr), naperts, TY_REAL) + call malloc (AGR_WOBS(agr), naperts, TY_REAL) + call malloc (AGR_OBS(agr), naperts, TY_REAL) + call malloc (AGR_WADO(agr), naperts, TY_REAL) + + call malloc (AGR_CUM(agr), naperts + 1, TY_REAL) + call malloc (AGR_TCUM(agr), naperts + 1, TY_REAL) + call malloc (AGR_WCUM(agr), naperts + 1, TY_REAL) + call malloc (AGR_MAGS(agr), naperts, TY_REAL) + call malloc (AGR_CMAGS(agr), naperts, TY_REAL) + call malloc (AGR_TMAGS(agr), naperts, TY_REAL) + call malloc (AGR_WMAGS(agr), naperts, TY_REAL) + + call calloc (AGR_AVE(agr), naperts, TY_REAL) + call calloc (AGR_RESID(agr), naperts, TY_REAL) + call calloc (AGR_RESSQ(agr), naperts, TY_REAL) + call calloc (AGR_WR(agr), naperts, TY_REAL) + call calloc (AGR_TAVE(agr), naperts, TY_REAL) + call calloc (AGR_TRESID(agr), naperts, TY_REAL) + call calloc (AGR_TRESSQ(agr), naperts, TY_REAL) + call calloc (AGR_TWR(agr), naperts, TY_REAL) +end + + +# PH_AMFIT -- Fit the seeing disk widths and the model parameters. + +int procedure ph_amfit (imtable, agr, nap, rap, mag, merr, naperts, params, + lterms) + +pointer imtable # pointer to the image symbol table +pointer agr # pointer to the fitting structure +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 parameters to be fit +int lterms # number of terms to fit + +double pold, sumd, sumw, sumn, sumr, sumy, old +int niter, i, j, k, nimages, ipow, nterms, iflag, ntimes, ngood +pointer sp, sym, symbol +real gain, rold, rclamp, dr +int stnsymbols() +pointer sthead(), stnext() + +begin + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) + return (ERR) + + # Allocate temporary space + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Order the symbol table. + symbol = sthead (imtable) + do k = nimages, 1, -1 { + Memi[sym+k-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Initialize some variables + ipow = 1 + nterms = 0 + gain = 0.0 + pold = 0.0d0 + call aclrd (Memd[AGR_POLD(agr)], lterms) + call ph_apinit (params, Memd[AGR_PARAMS(agr)], Memd[AGR_PERRORS(agr)], + Memd[AGR_PCLAMPS(agr)]) + + # Compute an improved value for ro and for the stellar profile + # model parameters. + + for (niter = 1; niter <= AGR_ITMAX; niter = niter + 1) { + + sumd = 0.0d0 + sumw = 0.0d0 + sumn = 0.0d0 + + # Compute the number of terms to be fit this iteration. + call ph_anterms (niter, lterms, gain, nterms) + + # Zero the accumulation matrix and vector + call aclrd (Memd[AGR_U(agr)], lterms * lterms) + call aclrd (Memd[AGR_V(agr)], lterms) + + do k = 1, nimages { + + # Get the current symbol. + symbol = Memi[sym+k-1] + + # Check to see if there is any data. + if (IMT_NENTRIES(symbol) <= 0) { + IMT_RO(symbol) = INDEFR + next + } + + # Check to see if there is any good data. + ngood = 0 + do i = 1, IMT_NENTRIES(symbol) { + do j = 2, nap[i] + ngood = ngood + 1 + } + if (ngood <= 0) { + IMT_RO(symbol) = INDEFR + next + } + + # Get better estimage of ro. + if (niter == 1) { + if (k > 1 && ! IS_INDEFR(IMT_RO(Memi[sym+k-2]))) { + IMT_RO(symbol) = IMT_RO(Memi[sym+k-2]) + } else + IMT_RO(symbol) = 0.5 * rap[1,1] + } + + # Set the model derivative vector + call ph_adddr (rap[1,IMT_OFFSET(symbol)], Memd[AGR_DDDR(agr)], + naperts, Memd[AGR_PARAMS(agr)], IMT_RO(symbol), + IMT_NXAIRMASS(symbol)) + + # Get the new estimate of ro + rold = 0.0 + rclamp = IMT_RO(symbol) / 4.0 + ntimes = 0 + repeat { + + sumr = 0.0d0 + sumy = 0.0d0 + + call ph_adddrdm (rap[1,IMT_OFFSET(symbol)], + Memd[AGR_DM(agr)], Memd[AGR_DDDR(agr)], naperts, + Memd[AGR_PARAMS(agr)], IMT_RO(symbol), + IMT_NXAIRMASS(symbol), niter) + + call ph_awsum (nap[IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], + merr[1,IMT_OFFSET(symbol)], Memd[AGR_DM(agr)], + Memr[AGR_W(agr)], Memd[AGR_DDDR(agr)], naperts, + IMT_NENTRIES(symbol), ipow, sumr, sumy) + + if (sumy > 0.0d0) { + dr = sumr / sumy + if (dr * rold < 0.0) + rclamp = 0.5 * rclamp + IMT_RO(symbol) = IMT_RO(symbol) - dr / (1.0 + abs (dr / + rclamp)) + if (IMT_RO(symbol) <= EPSILONR) { + IMT_RO(symbol) = INDEFR + call sfree (sp) + return (ERR) + } + rold = dr + if (abs (dr / IMT_RO(symbol)) <= 3.0e-4) + break + } else { + IMT_RO(symbol) = INDEFR + break + } + + ntimes = ntimes + 1 + if (ntimes >= 100) { + IMT_RO(symbol) = INDEFR + call sfree (sp) + return (ERR) + } + } + + call ph_adp (rap[1,IMT_OFFSET(symbol)], Memd[AGR_DM(agr)], + Memd[AGR_T(agr)], naperts, lterms, nterms, + Memd[AGR_PARAMS(agr)], IMT_RO(symbol), + IMT_NXAIRMASS(symbol)) + + call ph_aaccum (nap[IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], merr[1,IMT_OFFSET(symbol)], + Memd[AGR_DM(agr)], Memr[AGR_W(agr)], naperts, + IMT_NENTRIES(symbol), Memd[AGR_T(agr)], Memd[AGR_U(agr)], + Memd[AGR_V(agr)], lterms, nterms, ipow, sumd, sumw, sumn) + } + + # Invert the matrix. + call dinver (Memd[AGR_U(agr)], lterms, nterms, iflag) + if (iflag != 0) { + call sfree (sp) + return (ERR) + } + + # Get new values for the parameters. + call dmvmul (Memd[AGR_U(agr)], lterms, nterms, Memd[AGR_V(agr)], + Memd[AGR_DP(agr)]) + call ph_apsolve (Memd[AGR_POLD(agr)], Memd[AGR_PARAMS(agr)], + Memd[AGR_DP(agr)], Memd[AGR_PCLAMPS(agr)], nterms, ipow) + + # Test the fit. + #if (sumn > 6.0) + if (sumn > (nterms + 1.0)) + #sumn = sqrt (sumd / (sumn - 6.0)) + sumn = sqrt (sumd / (sumn - (nterms + 1.0))) + else + sumn = 0.0 + sumd = sqrt (sumd / sumw) + gain = 2.0 * abs (old - sumd) / (old + sumd) + old = sumd + if ((nterms < lterms) || (gain > 0.001)) { + next + } else if (ipow <= 2) { + ipow = 3 + call amovkd (0.1d0, Memd[AGR_PCLAMPS(agr)], 4) + next + } else + break + } + + if (niter > AGR_ITMAX) { + call sfree (sp) + return (ERR) + } + + # Compute the errors. + do k = 1, nterms { + if (Memd[AGR_U(agr)+(k-1)*lterms+k-1] > 0.0) + Memd[AGR_PERRORS(agr)+k-1] = + sumn * sqrt (Memd[AGR_U(agr)+(k-1)*lterms+k-1]) + else + Memd[AGR_PERRORS(agr)+k-1] = 0.0d0 + } + + + call sfree (sp) + return (OK) +end + + +# PH_ISHOW -- Show the initial values of the parameters. + +procedure ph_ishow (fd, params, mterms) + +int fd # pointer to the output file descriptor +double params[ARB] # the parameter values +int mterms # the number of terms to be fit + +int i + +begin + # Print out the best fit parameters. + call fprintf (fd, + "\nThe initial cog model parameter values for nparams %d\n") + call pargi (mterms) + do i = 1, MAX_MTERMS { + call fprintf (fd, "\t%10s: %15.7g\n") + switch (i) { + case 1: + call pargstr ("swings") + case 2: + call pargstr ("pwings") + case 3: + call pargstr ("pgauss") + case 4: + call pargstr ("rgescale") + case 5: + call pargstr ("xwings") + } + call pargd (params[i]) + } + call fprintf (fd, "\n") +end + + +# PH_PSHOW -- Show the results of the parameter fitting. + +procedure ph_pshow (fd, agr, mterms) + +int fd # the output file descriptor +pointer agr # the pointer to the fitting descriptor +int mterms # the number of parameters to fit + +begin + call ph_pwrite (fd, Memd[AGR_PARAMS(agr)], Memd[AGR_PERRORS(agr)], + MAX_MTERMS, mterms) +end + + +# PH_RSHOW -- Show the computed seeing parameters as a function of image. + +procedure ph_rshow (fd, imtable) + +int fd # the output file descriptor +pointer imtable # pointer to the symbol table + +int i, nimages +pointer sp, sym, symbol +int stnsymbols() +pointer sthead(), stnext() + +begin + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) + return + + # Allocate temporary space + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Order the symbol table. + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + call fprintf (fd, + "\nThe seeing radius and assumed airmass for each image\n") + do i = 1, nimages { + symbol = Memi[sym+i-1] + call fprintf (fd, "%30s %8.4f %8.3f\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (IMT_RO(symbol)) + call pargr (IMT_XAIRMASS(symbol)) + } + call fprintf (fd, "\n") + + call sfree (sp) +end + + +# PH_AEVAL -- Evaluate the fit for all the images. + +procedure ph_aeval (imtable, agr, apcat, magfd, logfd, mgd, id, x, y, nap, + rap, mag, merr, naperts, smallap, largeap) + +pointer imtable # pointer to the image symbol table +pointer agr # pointer to the fitting structure +int apcat # the aperture correction file descriptor +int magfd # the best magnitudes file descriptor +int logfd # the log file descriptor +pointer mgd # pointer to the plot metacode file +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 +int smallap # the small aperture number +int largeap # the large aperture number + +int k, nimages +pointer sp, sym, symbol +int stnsymbols() +pointer sthead(), stnext() + +begin + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) + return + + # Initialize some vectors. + call aclrr (Memr[AGR_TAVE(agr)], naperts) + call aclrr (Memr[AGR_TRESID(agr)], naperts) + call aclrr (Memr[AGR_TRESSQ(agr)], naperts) + call aclrr (Memr[AGR_TWR(agr)], naperts) + + # Allocate temporary space + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Order the symbol table. + symbol = sthead (imtable) + do k = nimages, 1, -1 { + Memi[sym+k-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Write the banner for the apfile. + call fprintf (apcat, + "# The aperture correction from apertures %d to %d in magnitudes\n\n") + call pargi (smallap) + call pargi (largeap) + + # Write the banner for the magfile. + if (magfd != NULL) { + call fprintf (magfd, "# Magnitudes corrected to aperture %d\n\n") + call pargi (largeap) + call fprintf (magfd, +"#%19tImage%30tFilter%39tExptime%47tAirmass%62tOtime%70tXcenter%80tYcenter\ +%93tMag%101tMerr%107tRadius\n\n") + } + + # Compute the aperture corrections. + do k = 1, nimages { + + symbol = Memi[sym+k-1] + + # Initialize some vectors. + call aclrr (Memr[AGR_AVE(agr)], naperts) + call aclrr (Memr[AGR_RESID(agr)], naperts) + call aclrr (Memr[AGR_RESSQ(agr)], naperts) + call aclrr (Memr[AGR_WR(agr)], naperts) + + call aclrr (Memr[AGR_ADOPT(agr)], naperts) + call aclrr (Memr[AGR_WADO(agr)], naperts) + call aclrr (Memr[AGR_WOBS(agr)], naperts) + + call ph_rinit (rap[1,IMT_OFFSET(symbol)], Memr[AGR_RBAR(agr)], + naperts) + call ph_tinit (rap[1,IMT_OFFSET(symbol)], Memr[AGR_THEO(agr)], + naperts, Memd[AGR_PARAMS(agr)], IMT_RO(symbol), + IMT_NXAIRMASS(symbol)) + + # Accumulate the difference data. + call ph_taccum (nap[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + merr[1,IMT_OFFSET(symbol)], Memr[AGR_THEO(agr)], + Memr[AGR_W(agr)], Memr[AGR_ADOPT(agr)], + Memr[AGR_WOBS(agr)], Memr[AGR_AVE(agr)], + Memr[AGR_RESID(agr)], Memr[AGR_RESSQ(agr)], + Memr[AGR_WR(agr)], Memr[AGR_OBS(agr)], Memr[AGR_WADO(agr)], + naperts, IMT_NENTRIES(symbol)) + + # Compute the cumulative differences. + call ph_tcum (rap[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)], + Memr[AGR_WADO(agr)], Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)], + Memr[AGR_WCUM(agr)], naperts, Memd[AGR_PARAMS(agr)], + IMT_RO(symbol), IMT_NXAIRMASS(symbol)) + + # Write the aperture photometry file. + call ph_wapcat (apcat, IMT_IMNAME(symbol), IMT_RO(symbol), + Memr[AGR_ADOPT(agr)], Memr[AGR_WADO(agr)], naperts, smallap, + largeap) + + if (logfd != NULL) { + call ph_twrite (logfd, IMT_IMNAME(symbol), IMT_RO(symbol), + IMT_XAIRMASS(symbol), Memr[AGR_RBAR(agr)], + Memr[AGR_THEO(agr)], Memr[AGR_OBS(agr)], + Memr[AGR_WOBS(agr)], Memr[AGR_ADOPT(agr)], + Memr[AGR_WADO(agr)], rap[1,IMT_OFFSET(symbol)], + Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)], + Memr[AGR_WCUM(agr)], naperts) + call fprintf (logfd, "\n") + call ph_tmags (logfd, id[IMT_OFFSET(symbol)], + x[IMT_OFFSET(symbol)], y[IMT_OFFSET(symbol)], + nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], merr[1,IMT_OFFSET(symbol)], + Memr[AGR_ADOPT(agr)], Memr[AGR_W(agr)], + Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)], + Memr[AGR_WCUM(agr)], Memr[AGR_MAGS(agr)], + Memr[AGR_CMAGS(agr)], Memr[AGR_TMAGS(agr)], + Memr[AGR_WMAGS(agr)], naperts, IMT_NENTRIES(symbol)) + call ph_papcor (logfd, IMT_IMNAME(symbol), IMT_RO(symbol), + rap[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)], + Memr[AGR_WADO(agr)], naperts, smallap, largeap) + call fprintf (logfd, "\n") + } + + if (magfd != NULL) { + call ph_wmags (magfd, symbol, x[IMT_OFFSET(symbol)], + y[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)], + rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + merr[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)], + Memr[AGR_WADO(agr)], Memr[AGR_W(agr)], Memr[AGR_CUM(agr)], + Memr[AGR_WCUM(agr)], Memr[AGR_MAGS(agr)], + Memr[AGR_CMAGS(agr)], Memr[AGR_WMAGS(agr)], naperts, + IMT_NENTRIES(symbol), largeap) + } + + if (mgd != NULL) { + call ph_gimfit (mgd, 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), YES) + if (! IS_INDEFR(IMT_RO(symbol))) { + call ph_gaimres (mgd, 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), YES) + call ph_gbimres (mgd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)], + mag[1,IMT_OFFSET(symbol)], naperts, + IMT_NENTRIES(symbol), YES) + call ph_gaximres (mgd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)], + x[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + naperts, IMT_NENTRIES(symbol), YES) + call ph_gayimres (mgd, agr, IMT_IMNAME(symbol), + IMT_RO(symbol), IMT_XAIRMASS(symbol), + nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)], + y[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)], + naperts, IMT_NENTRIES(symbol), YES) + call ph_gacum (mgd, 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), smallap, largeap, YES) + } + } + + if (! IS_INDEFR(IMT_RO(symbol))) { + call aaddr (Memr[AGR_AVE(agr)], Memr[AGR_TAVE(agr)], + Memr[AGR_TAVE(agr)], naperts) + call aaddr (Memr[AGR_RESID(agr)], Memr[AGR_TRESID(agr)], + Memr[AGR_TRESID(agr)], naperts) + call aaddr (Memr[AGR_RESSQ(agr)], Memr[AGR_TRESSQ(agr)], + Memr[AGR_TRESSQ(agr)], naperts) + call aaddr (Memr[AGR_WR(agr)], Memr[AGR_TWR(agr)], + Memr[AGR_TWR(agr)], naperts) + } + } + + call sfree (sp) +end + + +# PH_A1EVAL -- Evaluate the fit for all the images. + +procedure ph_a1eval (agr, image, r0, xairmass, nap, rap, mag, merr, + naperts, npts) + +pointer agr # pointer to the fitting structure +char image[ARB] # the image name +real r0 # the seeing radius +real xairmass # the airmass value +int nap[ARB] # the number of apertures array +real rap[naperts,ARB] # the list of aperture radii +real mag[naperts,ARB] # the magnitude difference array +real merr[naperts,ARB] # the magnitude error array +int naperts # the number of apertures +int npts # the number of points + +begin + # Initialize some vectors. + call aclrr (Memr[AGR_AVE(agr)], naperts) + call aclrr (Memr[AGR_RESID(agr)], naperts) + call aclrr (Memr[AGR_RESSQ(agr)], naperts) + call aclrr (Memr[AGR_WR(agr)], naperts) + call aclrr (Memr[AGR_ADOPT(agr)], naperts) + call aclrr (Memr[AGR_WOBS(agr)], naperts) + + # Compute the theoretical curve + call ph_rinit (rap, Memr[AGR_RBAR(agr)], naperts) + call ph_tinit (rap, Memr[AGR_THEO(agr)], naperts, + Memd[AGR_PARAMS(agr)], r0, xairmass) + + # Accumulate the difference data. + call ph_taccum (nap, mag, merr, Memr[AGR_THEO(agr)], Memr[AGR_W(agr)], + Memr[AGR_ADOPT(agr)], Memr[AGR_WOBS(agr)], Memr[AGR_AVE(agr)], + Memr[AGR_RESID(agr)], Memr[AGR_RESSQ(agr)], Memr[AGR_WR(agr)], + Memr[AGR_OBS(agr)], Memr[AGR_WADO(agr)], naperts, npts) + + # Compute the cumulative differences. + call ph_tcum (rap, Memr[AGR_ADOPT(agr)], Memr[AGR_WADO(agr)], + Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)], Memr[AGR_WCUM(agr)], + naperts, Memd[AGR_PARAMS(agr)], r0, xairmass) +end + + +# PH_TFSHOW -- Compute and print the summary statistics. + +procedure ph_tfshow (fd, agr, naperts) + +int fd # the output file descriptor +pointer agr # the pointer to the fit structure +int naperts # the number of apertures + +begin + call ph_rwrite (fd, Memr[AGR_TAVE(agr)], Memr[AGR_TWR(agr)], + Memr[AGR_TRESID(agr)], Memr[AGR_TRESSQ(agr)], naperts) +end + + +# PH_AMEMFREE -- Free the memory used for doing the fit + +procedure ph_amemfree (agr) + +pointer agr # pointer to the fitting structure + +begin + call mfree (AGR_DM(agr), TY_DOUBLE) + call mfree (AGR_DDDR(agr), TY_DOUBLE) + call mfree (AGR_T(agr), TY_DOUBLE) + call mfree (AGR_U(agr), TY_DOUBLE) + call mfree (AGR_V(agr), TY_DOUBLE) + call mfree (AGR_POLD(agr), TY_DOUBLE) + call mfree (AGR_DP(agr), TY_DOUBLE) + call mfree (AGR_PARAMS(agr), TY_DOUBLE) + call mfree (AGR_PERRORS(agr), TY_DOUBLE) + call mfree (AGR_PCLAMPS(agr), TY_DOUBLE) + + call mfree (AGR_RBAR(agr), TY_REAL) + call mfree (AGR_W(agr), TY_REAL) + call mfree (AGR_THEO(agr), TY_REAL) + call mfree (AGR_WR(agr), TY_REAL) + call mfree (AGR_ADOPT(agr), TY_REAL) + call mfree (AGR_WOBS(agr), TY_REAL) + call mfree (AGR_OBS(agr), TY_REAL) + call mfree (AGR_WADO(agr), TY_REAL) + + call mfree (AGR_CUM(agr), TY_REAL) + call mfree (AGR_TCUM(agr), TY_REAL) + call mfree (AGR_WCUM(agr), TY_REAL) + call mfree (AGR_MAGS(agr), TY_REAL) + call mfree (AGR_CMAGS(agr), TY_REAL) + call mfree (AGR_TMAGS(agr), TY_REAL) + call mfree (AGR_WMAGS(agr), TY_REAL) + + call mfree (AGR_AVE(agr), TY_REAL) + call mfree (AGR_RESID(agr), TY_REAL) + call mfree (AGR_RESSQ(agr), TY_REAL) + call mfree (AGR_TAVE(agr), TY_REAL) + call mfree (AGR_TRESID(agr), TY_REAL) + call mfree (AGR_TRESSQ(agr), TY_REAL) + call mfree (AGR_TWR(agr), TY_REAL) + + call mfree (agr, TY_STRUCT) +end + + +# PH_AGETP -- Fetch the initial values of the parameters for the cog model. + +procedure ph_agetp (params) + +double params[ARB] + +real clgetr() + +begin + params[1] = clgetr ("swings") + params[2] = clgetr ("pwings") + params[3] = clgetr ("pgauss") + params[4] = clgetr ("rgescale") + params[5] = clgetr ("xwings") +end + + +# PH_APINIT -- Set the initial values for the curve of growth model parameters. + +procedure ph_apinit (iparams, oparams, perrors, pclamps) + +double iparams[ARB] # the input parameters array +double oparams[ARB] # the output parameters array +double perrors[ARB] # the parameter errors array +double pclamps[ARB] # array of parameter clamps + +begin + call amovd (iparams, oparams, MAX_MTERMS) + call aclrd (perrors, MAX_MTERMS) + call amovkd (0.2d0, pclamps, MAX_MTERMS) + pclamps[4] = 0.5d0 +end + + +# PH_ANTERMS -- Compute the number of terms to use in the fit. + +procedure ph_anterms (niter, lterms, gain, nterms) + +int niter # the current iteration +int lterms # the maximum number of terms to be fit +real gain # the current value of the gain +int nterms # the current number of terms to be fit + +int n + +begin + n = 1 + if (nterms >= 1 && gain <= 0.04) { + n = 2 + if (nterms >= 2 && gain <= 0.016) { + n = 3 + if (nterms >= 3 && gain <= 0.006) { + n = 4 + if (nterms >= 4 && gain <= 0.0025) + n = 5 + } + } + } + + nterms = min (n, lterms) +end + + +# PH_ADDDR -- Compute the derivative wrt ro vector. + +procedure ph_adddr (rap, dddr, naperts, params, r0, xairmass) + +real rap[ARB] # array of aperture radii +double dddr[ARB] # the output derivatives array +int naperts # the number of apertures +double params[ARB] # the input parameter array +real r0 # the current value of r0 +real xairmass # the current value of the airmass + +int i +double ph_dmag() + +begin + do i = 2, naperts { + dddr[i] = 500.0d0 * (ph_dmag (rap[i-1], rap[i], xairmass, + r0 - 0.001, params) - ph_dmag (rap[i-1], rap[i], xairmass, + r0 + 0.001, params)) + } +end + + +# PH_ADDDRDM -- Compute the derivative wrt ro vector. + +procedure ph_adddrdm (rap, dm, dddr, naperts, params, r0, xairmass, niter) + +real rap[ARB] # array of aperture radii +double dm[ARB] # the output model differences array +double dddr[ARB] # the output derivatives array +int naperts # the number of apertures +double params[ARB] # the input parameter array +real r0 # the current value of r0 +real xairmass # the current value of the airmass +int niter # the current iteratiom + +int i +double ph_dmag() + +begin + do i = 2, naperts { + dm[i] = ph_dmag (rap[i-1], rap[i], xairmass, r0, params) + #if (niter == 1) + dddr[i] = 500.0d0 * (ph_dmag (rap[i-1], rap[i], xairmass, + r0 - 0.001, params) - ph_dmag (rap[i-1], rap[i], xairmass, + r0 + 0.001, params)) + } +end + + +# PH_AWSUM -- Accumulate the weighted sums necessary to fit r0. + +procedure ph_awsum (nap, mag, merr, dm, w, dddr, naperts, npts, ipow, + sumr, sumy) + +int nap[ARB] # array of aperture numbers +real mag[naperts,ARB] # array of magnitude difference +real merr[naperts,ARB] # array of magnitude errors +double dm[ARB] # array of model differences +real w[ARB] # array of weights +double dddr[ARB] # array of model derivatives +int naperts # number of apertures +int npts # number of data points +int ipow # weighting power factor +double sumr, sumy # the output sums + +int i, j +real diff, wt + +begin + do i = 1, npts { + w[1] = 1.0d0 + if (nap[i] <= 1) + next + do j = 2, nap[i] { + diff = mag[j,i] - dm[j] + w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** + ipow) + w[j] = min (w[j], w[j-1]) + wt = w[j] / merr[j,i] ** 2 + sumr = sumr + wt * diff * dddr[j] + sumy = sumy + wt * dddr[j] ** 2 + } + } +end + + +# PH_ADP -- Compute the parameter derivative vector with the new +# value of r0. + +procedure ph_adp (rap, dm, t, naperts, lterms, nterms, params, r0, xairmass) + +real rap[ARB] # array of aperture radii +double dm[ARB] # the new model differences +double t[lterms,ARB] # the parameter derivatives matrix +int naperts # the number of apertures +int lterms # the size of the derivatives matrix +int nterms # the number of terms to be fit +double params[ARB] # the current model parameter values +real r0 # the current r0 value +real xairmass # the current airmass + +int i, j +double ph_dmag() + +begin + do j = 2, naperts { + dm[j] = ph_dmag (rap[j-1], rap[j], xairmass, r0, params) + do i = 1, nterms { + params[i] = params[i] - 0.001 + t[i,j] = ph_dmag (rap[j-1], rap[j], xairmass, r0, params) + params[i] = params[i] + 0.002 + t[i,j] = 500.0d0 * (t[i,j] - ph_dmag (rap[j-1], rap[j], + xairmass, r0, params)) + params[i] = params[i] - 0.001 + } + } +end + + +# PH_AACCUM -- Accumulate the matrix and vector required to solve for the +# parameter increments + +procedure ph_aaccum (nap, mag, merr, dm, w, naperts, npts, t, u, v, lterms, + nterms, ipow, sumd, sumw, sumn) + +int nap[ARB] # array of aperture numbers +real mag[naperts,ARB] # array of magnitude difference +real merr[naperts,ARB] # array of magnitude errors +double dm[ARB] # array of model differences +real w[ARB] # array of weights +int naperts # number of apertures +int npts # number of data points +double t[lterms,ARB] # the array of parameter derivatives +double u[lterms,ARB] # the matrix to be accumulated +double v[ARB] # the vector to be accumulated +int lterms # the maximum number of terms to fit +int nterms # the current number of terms to fit +int ipow # power factor for the weights +double sumd # sum of the differences +double sumw # sum of the scatter +double sumn # sum of the weights + +int i, j, l, m +real diff, wt + +begin + do i = 1, npts { + w[1] = 1.0 + if (nap[i] <= 1) + next + do j = 2, nap[i] { + diff = mag[j,i] - dm[j] + w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** ipow) + w[j] = min (w[j], w[j-1]) + wt = w[j] / merr[j,i] ** 2 + sumd = sumd + wt * diff ** 2 + sumw = sumw + wt + sumn = sumn + w[j] + do l = 1, nterms { + v[l] = v[l] + wt * diff * t[l,j] + do m = 1, nterms + u[l,m] = u[l,m] + wt * t[l,j] * t[m,j] + } + } + } +end + + +# PH_APSOLVE -- Solve for new values of the parameters. + +procedure ph_apsolve (pold, params, dp, pclamps, nterms, ipow) + +double pold[ARB] # the previous parameter imcrement values +double params[ARB] # the current values of the parameters +double dp[ARB] # the parameter increment values +double pclamps[ARB] # the parameter clamp values +int nterms # the number of terms that were fit +int ipow # the current weighting factor + +int i + +begin + do i = 1, nterms { + if (dp[i] * pold[i] < 0.0d0) + pclamps[i] = 0.5 * pclamps[i] + pold[i] = dp[i] + } + + params[1] = params[1] - dp[1] / ( 1.0d0 + abs (dp[1] / (pclamps[1] * + (params[1] - 1.0d0)))) + if (nterms >= 2) { + params[2] = params[2] - dp[2] / (1.0d0 + abs (dp[2] / (pclamps[2] * + params[2] * (1.0d0 - params[2])))) + if (nterms >= 3) { + params[3] = params[3] - dp[3] / (1.0d0 + abs (dp[3] / + (pclamps[3] * params[3] * (1.0d0 - params[3])))) + if (ipow == 1) { + ipow = 2 + pclamps[1] = 0.1d0 + pclamps[2] = 0.1d0 + } + if (nterms >= 4) { + params[4] = params[4] - dp[4] / (1.0d0 + abs (dp[4] / + (pclamps[4] * params[4]))) + if (nterms >= 5) + params[5] = params[5] - dp[5] / (1.0d0 + abs (dp[5] / + (pclamps[5] * params[2]))) + } + } + } +end + +# PH_PWRITE -- Write out the theoetical model parameters. + +procedure ph_pwrite (fd, params, perrors, max_nterms, nterms) + +int fd # the output file descriptor +double params[ARB] # the current model parameter values +double perrors[ARB] # the current parameter errors +int max_nterms # the number of terms +int nterms # the number of terms to fit + +int i + +begin + # Print out the best fit parameters. + call fprintf (fd, + "\nThe computed cog model parameters and their errors for nparams %d\n") + call pargi (nterms) + do i = 1, max_nterms { + call fprintf (fd, "\t%10s: %15.7g +/- %15.7g\n") + switch (i) { + case 1: + call pargstr ("swings") + case 2: + call pargstr ("pwings") + case 3: + call pargstr ("pgauss") + case 4: + call pargstr ("rgescale") + case 5: + call pargstr ("xwings") + } + call pargd (params[i]) + call pargd (perrors[i]) + } + call fprintf (fd, "\n") +end + +# PH_RINIT -- Initialize the rbar vector. + +procedure ph_rinit (rap, rbar, naperts) + +real rap[ARB] # the array of aperture radii +real rbar[ARB] # the mean radius estimates +int naperts # the number of aperture radii + +int j + +begin + do j = 2, naperts + rbar[j] = 0.5 * (rap[j-1] + rap[j]) +end + +# PH_TINIT -- Initialize the rbar vector and compute the theoretical estimates. + +procedure ph_tinit (rap, theo, naperts, params, r0, airmass) + +real rap[ARB] # the array of aperture radii +real theo[ARB] # the current model estimates +int naperts # the number of aperture radii +double params[ARB] # the current parameter estimates +real r0 # the seeing radius +real airmass # the airmass value + +int j +double ph_dmag() + +begin + if (IS_INDEFR(r0)) + call amovkr (INDEFR, theo[2], naperts - 1) + else { + do j = 2, naperts + theo[j] = ph_dmag (rap[j-1], rap[j], airmass, r0, params) + } +end + + +# PH_TACCUM -- Accumulate the data. + +procedure ph_taccum (nap, mag, merr, theo, w, adopt, wobs, ave, resid, ressq, + wr, obs, wado, naperts, npts) + +int nap[ARB] # array of aperture numbers +real mag[naperts,ARB] # the array of magnitude differences +real merr[naperts,ARB] # the array of magnitude errors +real theo[ARB] # the current model values +real w[ARB] # the working weight array +real adopt[ARB] # the adopted weight differences +real wobs[ARB] # the sum of the weights +real ave[ARB] # the average model values +real resid[ARB] # the residuals +real ressq[ARB] # the residuals squared +real wr[ARB] # the sum of the weights +real obs[ARB] # the observations +real wado[ARB] # the error in the observations +int naperts # the number of apertures +int npts # the number of npts + +int i, j +real diff, wt, scale + +begin + do i = 1, npts { + w[1] = 1.0 + if (nap[i] <= 1) + next + do j = 2, nap[i] { + if (IS_INDEFR(theo[j])) { + wt = 1.0 / merr[j,i] ** 2 + ave[1] = INDEFR + ave[j] = INDEFR + resid[1] = INDEFR + resid[j] = INDEFR + ressq[1] = INDEFR + ressq[j] = INDEFR + wr[1] = INDEFR + wr[j] = INDEFR + } else { + diff = mag[j,i] - theo[j] + w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i]))) + w[j] = min (w[j], w[j-1]) + wt = w[j] / merr[j,i] ** 2 + ave[1] = ave[1] + wt * theo[j] + ave[j] = ave[j] + wt * theo[j] + resid[1] = resid[1] + wt * diff + resid[j] = resid[j] + wt * diff + ressq[1] = ressq[1] + wt * diff ** 2 + ressq[j] = ressq[j] + wt * diff ** 2 + wr[1] = wr[1] + wt + wr[j] = wr[j] + wt + } + adopt[j] = adopt[j] + wt * mag[j,i] + wobs[j] = wobs[j] + wt + } + } + + do j = 2, naperts { + if (wobs[j] <= 0.0) + next + obs[j] = adopt[j] / wobs[j] + adopt[j] = 0.0 + wobs[j] = 0.0 + } + + do i = 1, npts { + w[1] = 1.0 + if (nap[i] <= 1) + next + do j = 2, nap[i] { + diff = mag[j,i] - obs[j] + w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** 2) + w[j] = min (w[j], w[j-1]) + wt = w[j] / merr[j,i] ** 2 + adopt[j] = adopt[j] + wt * mag[j,i] + wobs[j] = wobs[j] + wt + } + } + + do j = 2, naperts { + if (wobs[j] <= 0.0) + next + obs[j] = adopt[j] / wobs[j] + adopt[j] = 0.0 + wobs[j] = 0.0 + } + + do i = 1, npts { + w[1] = 1.0 + if (nap[i] <= 1) + next + do j = 2, nap[i] { + diff = mag[j,i] - obs[j] + w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** 3) + w[j] = min (w[j], w[j-1]) + wt = w[j] / merr[j,i] ** 2 + adopt[j] = adopt[j] + wt * mag[j,i] + wobs[j] = wobs[j] + wt + } + } + + scale = 0.1 + do j = 2, naperts { + if (wobs[j] > 0.0) { + obs[j] = adopt[j] / wobs[j] + if (IS_INDEFR(theo[j])) + wt = 0.0 + else { + wt = 1.0 / (scale * theo[j]) ** 2 + adopt[j] = adopt[j] + wt * theo[j] + } + wt = 1.0 / (wobs[j] + wt) + adopt[j] = wt * adopt[j] + #wado[j] = sqrt (wt) + wado[j] = sqrt (wt + (scale * adopt[j]) ** 2) + wobs[j] = sqrt (1.0 / wobs[j]) + } else { + adopt[j] = theo[j] + if (IS_INDEFR(theo[j])) { + wado[j] = INDEFR + obs[j] = INDEFR + wobs[j] = INDEFR + } else + wado[j] = 2.0 * scale * abs (theo[j]) + } + } +end + + +# PH_TCUM -- Compute the cumulative differences. + +procedure ph_tcum (rap, adopt, wado, cum, tcum, wcum, naperts, params, + r0, airmass) + +real rap[ARB] # the list of aperture radii +real adopt[ARB] # the adopted differences +real wado[ARB] # the errors in the adopted differences +real cum[ARB] # the accumulated differences +real tcum[ARB] # the total accumulated differences +real wcum[ARB] # the accumulated difference errors +int naperts # the number of aperture radii +double params[ARB] # the current parameter values +real r0 # the seeing radius +real airmass # the airmass + +int i +double ph_dmag() + +begin + if (IS_INDEFR(r0)) { + call amovkr (INDEFR, cum, naperts + 1) + call amovkr (INDEFR, tcum, naperts + 1) + call amovkr (INDEFR, wcum, naperts + 1) + } else { + cum[naperts+1] = 0.0 + tcum[naperts+1] = ph_dmag (rap[naperts], 2.0 * rap[naperts], + airmass, r0, params) + wcum[naperts+1] = 0.0 + do i = naperts, 2, -1 { + cum[i] = adopt[i] + cum[i+1] + tcum[i] = adopt[i] + tcum[i+1] + wcum[i] = wado[i] ** 2 + wcum[i+1] + } + } +end + + +# PH_WAPCAT -- Write the image entry to the aperture correction catalog. + +procedure ph_wapcat (apcat, image, r0, adopt, wado, naperts, smallap, largeap) + +int apcat # the aperture correction catalog descriptor +char image[ARB] # the image name +real r0 # the seeing radius +real adopt[ARB] # the adopted difference values +real wado[ARB] # the adopted difference errors +int naperts # the number of apertures +int smallap # the small aperture number +int largeap # the large aperture number + +int i +real cum, wcum + +begin + if (IS_INDEFR(r0)) { + cum = INDEFR + wcum = INDEFR + } else { + cum = 0.0 + wcum = 0.0 + do i = largeap, smallap + 1, -1 { + cum = cum + adopt[i] + wcum = wcum + wado[i] ** 2 + } + wcum = sqrt (wcum) + } + call fprintf (apcat, "%s %g %g\n") + call pargstr (image) + call pargr (cum) + call pargr (wcum) +end + + +# PH_PAPCOR -- Print the aperture correction on the standard output. + +procedure ph_papcor (fd, image, r0, rap, adopt, wado, naperts, smallap, largeap) + +int fd # the output file descriptor +char image[ARB] # the image name +real r0 # the seeing radius +real rap[ARB] # the aperture radii +real adopt[ARB] # the adopted difference values +real wado[ARB] # the adopted difference errors +int naperts # the number of apertures +int smallap # the small aperture number +int largeap # the large aperture number + +int i +real cum, wcum + +begin + if (IS_INDEFR(r0)) { + cum = INDEFR + wcum = INDEFR + } else { + cum = 0.0 + wcum = 0.0 + do i = largeap, smallap + 1, -1 { + cum = cum + adopt[i] + wcum = wcum + wado[i] ** 2 + } + wcum = sqrt (wcum) + } + + call fprintf (fd, + "Image: %s rin=%.2f rout=%.2f apercor=%.3f +/- %.4f\n") + call pargstr (image) + call pargr (rap[smallap]) + call pargr (rap[largeap]) + call pargr (cum) + call pargr (wcum) +end + + +define NPERLINE 7 + +# PH_TWRITE -- Write the results to the output file. + +procedure ph_twrite (fd, image, r0, xairmass, rbar, theo, obs, wobs, adopt, + wado, rap, cum, tcum, wcum, naperts) + +int fd # the file descriptor +char image # the iamge name +real r0 # the seeing disk +real xairmass # the assumed airmass +real rbar[ARB] # the list of mean aperture radii +real theo[ARB] # the theoretical model differences +real obs[ARB] # the observed differences +real wobs[ARB] # the observed difference errors +real adopt[ARB] # the adopted differences +real wado[ARB] # the adopted difference errors +real rap[ARB] # the list of aperture radii +real cum[ARB] # the cumulative differences +real tcum[ARB] # the total cumulative differences +real wcum[ARB] # the errors in the cumulative differences +int naperts # the number of aperture + +int j +real sqwcum + +begin + # Print the title. + call fprintf (fd, "\nThe cog for image %s from radius %.2f to %.2f\n") + call pargstr (image) + call pargr (rbar[2]) + call pargr (rbar[naperts]) + + # Print the mean aperture radius. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " radius %9.4f") + call pargr (rbar[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (rbar[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (rbar[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (rbar[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the theoretical model. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " model %9.4f") + call pargr (theo[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (theo[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (theo[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (theo[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the observed cog. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, "observed %9.4f") + call pargr (obs[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (obs[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (obs[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (obs[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the errors in the observed cog. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " sigma %9.4f") + call pargr (wobs[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (wobs[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (wobs[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (wobs[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the adopted cog. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " adopted %9.4f") + call pargr (adopt[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (adopt[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (adopt[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (adopt[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the errors in the adopted cog. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " sigma %9.4f") + call pargr (wado[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (wado[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (wado[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (wado[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the title. + call fprintf (fd, + "\nThe aperture correction for image %s from radius %.2f to %.2f\n") + call pargstr (image) + call pargr (rap[1]) + call pargr (rap[naperts]) + + # Print the aperture radii. + do j = 1, naperts { + if (j == 1) { + call fprintf (fd, " radius %9.4f") + call pargr (rap[j]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (fd, "%9.4f\n") + call pargr (rap[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, " %9.4f") + call pargr (rap[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (rap[j]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (fd, "\n") + + # Print the aperture correction. + do j = 2, naperts { + if (j == 2) { + call fprintf (fd, " apercor %9.4f") + call pargr (cum[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (cum[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (cum[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (cum[j]) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + + # Print the error in the aperture correction. + do j = 2, naperts { + sqwcum = wcum[j] + if (j == 2) { + call fprintf (fd, " sigma %9.4f") + call pargr (sqwcum) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (sqwcum) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (sqwcum) + } else { + call fprintf (fd, "%9.4f") + call pargr (sqwcum) + } + } + if (mod (naperts, NPERLINE) != 1) + call fprintf (fd, "\n") + + # Print the title. + call fprintf (fd, + "\nThe aperture correction for image %s from radius %.2f to %.2f\n") + call pargstr (image) + call pargr (rap[1]) + call pargr (2.0 * rap[naperts]) + + # Print the total aperture correction. + do j = 2, naperts + 1 { + if (j == 2) { + call fprintf (fd, "tapercor %9.4f") + call pargr (tcum[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (fd, "%9.4f\n") + call pargr (tcum[j]) + } else if (mod (j, NPERLINE) == 2) { + call fprintf (fd, " %9.4f") + call pargr (tcum[j]) + } else { + call fprintf (fd, "%9.4f") + call pargr (tcum[j]) + } + } + if (mod (naperts + 1, NPERLINE) != 1) + call fprintf (fd, "\n") + + call fprintf (fd, "\n") +end + + +# PH_TMAGS -- Compute the correction to the last computed magnitude for +# each star and the total magnitude as a function of aperture. + +procedure ph_tmags (logfd, id, x, y, nap, rap, mag, merr, adopt, w, cum, + tcum, wcum, tmpmags, obs, tobs, wobs, naperts, npts) + +int logfd # the output file descriptor +int id[ARB] # stellar id numbers +real x[ARB] # stellar x coordinates +real y[ARB] # stellar y coordinates +int nap[ARB] # array of aperture numbers +real rap[naperts,ARB] # the list of aperture radii +real mag[naperts,ARB] # the input magnitude difference array +real merr[naperts,ARB] # the input magnitude error array +real adopt[ARB] # the adopted difference values +real w[ARB] # the working weight array +real cum[ARB] # the accumulated difference array +real tcum[ARB] # the total accumulated difference array +real wcum[ARB] # the errors in the accumulated diff array +real tmpmags[ARB] # temporary magnitude array +real obs[ARB] # the accumulated magnitude array +real tobs[ARB] # the total accumulated magnitude array +real wobs[ARB] # the observations array +int naperts # number of apertures +int npts # number of points + +int i, j +real diff + +begin + # Print the logfile banner. + if (logfd != NULL) { + call fprintf (logfd, "\nThe observed, and corrected to radii %.2f ") + call pargr (rap[naperts,1]) + call fprintf (logfd, "and %.2f, magnitudes\n") + call pargr (2.0 * rap[naperts,1]) + } + + do i = 1, npts { + + # Compute the observed, correctged and estimated total + # magnitudes. + tmpmags[1] = mag[1,i] + do j = 1, nap[i] { + if (j == 1) + w[j] = 1.0 + else { + diff = mag[j,i] - adopt[j] + tmpmags[j] = tmpmags[j-1] + mag[j,i] + w[j] = 1.0 / (1.0 + (diff / (2.0 * merr[j,i])) ** 2) + w[j] = min (w[j], w[j-1]) + } + if (IS_INDEFR(cum[j+1])) { + obs[j] = INDEFR + tobs[j] = INDEFR + wobs[j] = INDEFR + } else { + obs[j] = tmpmags[j] + cum[j+1] + tobs[j] = tmpmags[j] + tcum[j+1] + wobs[j] = sqrt (wcum[j+1] + merr[j,i] ** 2 / w[j]) + } + } + do j = nap[i] + 1, naperts { + obs[j] = INDEFR + tobs[j] = INDEFR + wobs[j] = INDEFR + } + + # Write out the results for the star to the log file. + + # Print the banner of the star. + call fprintf (logfd, "Star: %d x: %.3f y: %.3f\n") + call pargi (id[i]) + call pargr (x[i]) + call pargr (y[i]) + + # Print the aperture radii. + do j = 1, naperts { + if (j == 1) { + call fprintf (logfd, " radius %9.4f") + call pargr (rap[j,i]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (logfd, "%9.4f\n") + call pargr (rap[j,i]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (logfd, " %9.4f") + call pargr (rap[j,i]) + } else { + call fprintf (logfd, "%9.4f") + call pargr (rap[j,i]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (logfd, "\n") + + # Print the observed magnitudes. + do j = 1, naperts { + if (j == 1) { + call fprintf (logfd, " mags %9.4f") + call pargr (tmpmags[j]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (logfd, "%9.4f\n") + call pargr (tmpmags[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (logfd, " %9.4f") + call pargr (tmpmags[j]) + } else { + call fprintf (logfd, "%9.4f") + call pargr (tmpmags[j]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (logfd, "\n") + + # Print the corrected magnitudes. + do j = 1, naperts { + if (j == 1) { + call fprintf (logfd, " cmags %9.4f") + call pargr (obs[j]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (logfd, "%9.4f\n") + call pargr (obs[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (logfd, " %9.4f") + call pargr (obs[j]) + } else { + call fprintf (logfd, "%9.4f") + call pargr (obs[j]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (logfd, "\n") + + # Print the estimated total corrected magnitudes. + do j = 1, naperts { + if (j == 1) { + call fprintf (logfd, " tcmags %9.4f") + call pargr (tobs[j]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (logfd, "%9.4f\n") + call pargr (tobs[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (logfd, " %9.4f") + call pargr (tobs[j]) + } else { + call fprintf (logfd, "%9.4f") + call pargr (tobs[j]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (logfd, "\n") + + # Print the errors in the total magnitudes + do j = 1, naperts { + if (j == 1) { + call fprintf (logfd, " sigma %9.4f") + call pargr (wobs[j]) + } else if (mod (j, NPERLINE) == 0) { + call fprintf (logfd, "%9.4f\n") + call pargr (wobs[j]) + } else if (mod (j, NPERLINE) == 1) { + call fprintf (logfd, " %9.4f") + call pargr (wobs[j]) + } else { + call fprintf (logfd, "%9.4f") + call pargr (wobs[j]) + } + } + if (mod (naperts, NPERLINE) != 0) + call fprintf (logfd, "\n") + } + + call fprintf (logfd, "\n") +end + + +# PH_WMAGS -- Compute the correction to the last computed magnitude for +# each star and the total magnitude as a function of aperture. + +procedure ph_wmags (magfd, imsymbol, x, y, nap, rap, mag, merr, adopt, wado, + w, cum, wcum, tmpmags, obs, wobs, naperts, npts, largeap) + +int magfd # the best magnitudes file descriptor +pointer imsymbol # the image symbol +real x[ARB] # stellar x coordinates +real y[ARB] # stellar y coordinates +int nap[ARB] # array of aperture numbers +real rap[naperts,ARB] # the input aperture radii +real mag[naperts,ARB] # the input magnitude difference array +real merr[naperts,ARB] # the input magnitude error array +real adopt[ARB] # the adopted difference values +real wado[ARB] # the adopted difference errors +real w[ARB] # the working weight array +real cum[ARB] # the accumulated difference array +real wcum[ARB] # the errors in the accumulated diff array +real tmpmags[ARB] # temporary magnitude array +real obs[ARB] # the accumulated magnitude array +real wobs[ARB] # the observations array +int naperts # number of apertures +int npts # number of points +int largeap # the largest aperture + +int i, j, jfinal +real diff, sigmin, sigcor +real assqr() + +begin + do i = 1, npts { + + # Compute the observed, correctged and estimated total magnitudes. + tmpmags[1] = mag[1,i] + sigmin = MAX_REAL + #jfinal = min (nap[i], largeap) + jfinal = largeap + if (largeap < nap[i]) + sigcor = assqr (wado[largeap+1], nap[i] - largeap) + else + sigcor = 0.0 + do j = 1, min (nap[i], largeap) { + if (j == 1) + w[j] = 1.0 + else { + diff = mag[j,i] - adopt[j] + tmpmags[j] = tmpmags[j-1] + mag[j,i] + w[j] = 1.0 / (1.0 + (diff / (2.0 * merr[j,i])) ** 2) + w[j] = min (w[j], w[j-1]) + } + if (IS_INDEFR(cum[j+1])) { + obs[j] = INDEFR + wobs[j] = INDEFR + } else { + obs[j] = tmpmags[j] + cum[j+1] - cum[largeap+1] + wobs[j] = sqrt (wcum[j+1] - sigcor + merr[j,i] ** 2 / w[j]) + if (wobs[j] < sigmin) { + jfinal = j + sigmin = wobs[j] + } + } + } + + # Write out the best magnitude. + if (magfd != NULL) { + call fprintf (magfd, +"%23.23s %10.10s %8.2f %6.3f %11.1h %8.2f %8.2f %7.3f %7.3f %6.3f\n") + call pargstr (IMT_IMNAME(imsymbol)) + call pargstr (IMT_IFILTER(imsymbol)) + call pargr (IMT_ITIME(imsymbol)) + call pargr (IMT_XAIRMASS(imsymbol)) + call pargr (IMT_OTIME(imsymbol)) + call pargr (x[i]) + call pargr (y[i]) + call pargr (obs[jfinal]) + call pargr (wobs[jfinal]) + call pargr (rap[jfinal,i]) + } + } +end + + +# PH_ACHI -- Compute the chi statistic. + +real procedure ph_achi (wr, resid, ressq, naperts) + +real wr[ARB] # the weights +real resid[ARB] # the residuals +real ressq[ARB] # the residuals +int naperts # the number of apertures + +double sumwr, sumres, sumressq +real chi +real asumr() + +begin + sumwr = asumr (wr, naperts) + sumres = asumr (resid, naperts) + sumressq = asumr (ressq, naperts) + if (sumwr <= 0.0d0) + chi = 0.0 + else + chi = sumressq / sumwr - (sumres / sumwr) ** 2 + if (chi <= 0.0 || chi > MAX_REAL) + return (INDEFR) + else + return (chi) +end + + +# PH_RWRITE -- Write out the fit statistics for each aperture. + +procedure ph_rwrite (fd, ave, wr, resid, ressq, naperts) + +int fd # the output file descriptor +real ave[ARB] # the average values +real wr[ARB] # the weights +real resid[ARB] # the residuals +real ressq[ARB] # the residuals +int naperts # the number of apertures + +int j +real rtmp + +begin + call fprintf (fd, "\nAverage model cog, residual, and rms residual ") + call fprintf (fd, "for each aperture all images\n") + do j = 1, naperts { + if (wr[j] <= 0.0) { + ave[j] = INDEFR + resid[j] = INDEFR + ressq[j] = INDEFR + } else { + ave[j] = ave[j] / wr[j] + resid[j] = resid[j] / wr[j] + rtmp = ressq[j] / wr[j] - resid[j] ** 2 + if (rtmp <= 0.0) + ressq[j] = 0.0 + else + ressq[j] = sqrt (rtmp) + } + call fprintf (fd, "\t%3d %9.5f %9.5f %9.5f\n") + call pargi (j) + call pargr (ave[j]) + call pargr (resid[j]) + call pargr (ressq[j]) + } + call fprintf (fd, "\n") +end + + +# PH_DMAG -- Compute the integral of the magnitude between two radii +# assuming that the stellar profile can be approximated by a circular +# Moffat function. + +double procedure ph_dmag (r1, r2, x, r0, params) + +real r1 # beginning radius for the integration +real r2 # ending radius for the integration +real x # the airmass value +real r0 # the hwhm of the psf +double params[ARB] # the model parameter array + +double bpex, pm1, x1, x2, x1sq, x2sq, d1, d2, i1, i2, dmag + +begin + bpex = params[2] + params[5] * x + pm1 = params[1] - 1.0d0 + x1 = (r1 / r0) + x2 = (r2 / r0) + x1sq = x1 ** 2 + x2sq = x2 ** 2 + x1 = x1 / params[4] + x2 = x2 / params[4] + d1 = 1.0d0 + r1 ** 2 + d2 = 1.0d0 + r2 ** 2 + i1 = bpex * (1.0d0 - 1.0d0 / d1 ** pm1) + (1.0d0 - bpex) * + (params[3] * (1.0d0 - exp (-0.5d0 * x1sq)) + (1.0d0 - params[3]) * + (1.0d0 - (1.0d0 + x1) * exp (-x1))) + i2 = bpex * (1.0d0 - 1.0d0 / d2 ** pm1) + (1.0d0 - bpex) * + (params[3] * (1.0d0 - exp (-0.5d0 * x2sq)) + (1.0d0 - params[3]) * + (1.0d0 - (1.0d0 + x2) * exp (-x2))) + dmag = -2.5d0 * log10 (i2 / i1) + + return (dmag) +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/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 diff --git a/noao/digiphot/photcal/mkobsfile/phaimtable.x b/noao/digiphot/photcal/mkobsfile/phaimtable.x new file mode 100644 index 00000000..a1e8cd94 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phaimtable.x @@ -0,0 +1,409 @@ +include <fset.h> +include "../lib/apfile.h" + +# PH_EAOBSIMTABLE -- Enter the correct values of the filterid, exposure time, +# and airmass into the image symbol table from the standard input. + +procedure ph_eaobsimtable (imtable, verify) + +pointer imtable # pointer to the symbol table +int verify # verify the user input + +int nimages, i +pointer sp, sym, filterid, symbol +real itime, xairmass, otime +int stnsymbols(), scan(), nscan(), strncmp() +pointer sthead(), stnext() + +begin + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) + return + + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + # Reverse the order of the symbols in the symbol table. + + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + call printf ("\n") + call printf ("Enter filterid, exposure time and airmass\n") + + # Loop over the images. + do i = 1, nimages { + + symbol = Memi[sym+i-1] + + # Issue the prompt for each image. + call printf ( + " Image %s (f t X T, <CR>=4 X INDEF, <EOF>=quit entry): " ) + call pargstr (IMT_IMNAME(symbol)) + call flush (STDOUT) + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + break + } + + # Read in the airmass. + call gargwrd (Memc[filterid], SZ_FNAME) + call gargr (itime) + call gargr (xairmass) + call gargr (otime) + if (nscan() < 1) { + call strcpy ("INDEF", Memc[filterid], SZ_FNAME) + itime = INDEFR + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 2) { + itime = INDEFR + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 3) { + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 4) { + otime = INDEFR + } + + # Update the symbol table. + if (strncmp (Memc[filterid], "INDEF", 5) != 0) + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + if (! IS_INDEFR(itime)) + IMT_ITIME(symbol) = itime + if (! IS_INDEFR(xairmass)) { + IMT_XAIRMASS(symbol) = xairmass + IMT_NXAIRMASS(symbol) = xairmass - DEF_AIROFFSET + } else { + IMT_XAIRMASS(symbol) = DEF_AIROFFSET + IMT_NXAIRMASS(symbol) = 0.0 + call printf ( + " Warning: Setting airmass for image %s to %g\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (DEF_AIROFFSET) + } + if (! IS_INDEFR(otime)) + IMT_OTIME(symbol) = otime + + # Verify the input. + if (verify == NO) + next + + # Issue the verify prompt. + call printf (" Verify (f t X T, <CR>=%s %g %g %0.1h): ") + call pargstr (IMT_IFILTER(symbol)) + call pargr (IMT_ITIME(symbol)) + call pargr (IMT_XAIRMASS(symbol)) + call pargr (IMT_OTIME(symbol)) + call flush (STDOUT) + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + next + } + + # Read in the airmass. + call gargwrd (Memc[filterid], SZ_FNAME) + call gargr (itime) + call gargr (xairmass) + call pargr (otime) + if (nscan() == 4) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + if (! IS_INDEFR(xairmass)) { + IMT_XAIRMASS(symbol) = xairmass + IMT_NXAIRMASS(symbol) = xairmass - DEF_AIROFFSET + } else { + IMT_XAIRMASS(symbol) = DEF_AIROFFSET + IMT_NXAIRMASS(symbol) = 0.0 + call printf ( + " Warning: Setting airmass for image %s to %g\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (DEF_AIROFFSET) + } + IMT_OTIME(symbol) = otime + } else if (nscan() == 3) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + if (! IS_INDEFR(xairmass)) { + IMT_XAIRMASS(symbol) = xairmass + IMT_NXAIRMASS(symbol) = xairmass - DEF_AIROFFSET + } else { + IMT_XAIRMASS(symbol) = DEF_AIROFFSET + IMT_NXAIRMASS(symbol) = 0.0 + call printf ( + " Warning: Setting airmass for image %s to %g\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (DEF_AIROFFSET) + } + } else if (nscan() == 2) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + } else if (nscan() == 1) + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + } + + call printf ("\n") + + call sfree (sp) +end + + +# PH_AOBSIMTABLE -- Enter the correct values of the filterid, exposure time, +# and airmass into the image table. + +procedure ph_aobsimtable (imtable, fd, columns) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns + +int stat, lbufsize +pointer sp, image, fname, filterid, line, sym +real itime, airmass, otime +bool streq() +int ph_anxtimage() +pointer stfind() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + Memc[image] = EOS + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + line = NULL + lbufsize = 0 + + repeat { + + # Get the image name. + stat = ph_anxtimage (fd, columns, Memc[image], line, lbufsize) + if (stat == EOF) + break + + # Locate the image name in the symbol table. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + call printf ( + "Warning: File: %s image: %s ") + call pargstr (Memc[fname]) + call pargstr (Memc[image]) + call printf ("is not in the input image list\n") + next + } + + # Decode the data. + call ph_aobsdata (Memc[line], columns, Memc[filterid], SZ_FNAME, + itime, airmass, otime) + + # Enter the data in the symbol table. + if (! streq (Memc[filterid], "INDEF")) + call strcpy (Memc[filterid], IMT_IFILTER(sym), SZ_FNAME) + if (! IS_INDEFR(itime)) + IMT_ITIME(sym) = itime + if (! IS_INDEFR(airmass)) { + IMT_XAIRMASS(sym) = airmass + IMT_NXAIRMASS(sym) = airmass - DEF_AIROFFSET + } else { + IMT_XAIRMASS(sym) = DEF_AIROFFSET + IMT_NXAIRMASS(sym) = 0.0 + call printf ( + " Warning: Setting airmass for image %s to %g\n") + call pargstr (IMT_IMNAME(sym)) + call pargr (DEF_AIROFFSET) + } + if (! IS_INDEFR(otime)) + IMT_OTIME(sym) = otime + } + + if (line != NULL) { + call mfree (line, TY_CHAR) + line = NULL + lbufsize = 0 + } + + call sfree (sp) +end + + +# PH_AMKIMTABLE -- Create the image table using the input text file(s) +# and the list of column numbers defining the fields required by the +# program. + +int procedure ph_amkimtable (imtable, fd, columns, naperts, imid, id, x, y, + nap, rap, mag, merr, nptr, sortimid, magerr) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +int naperts # the number of apertures +pointer imid # pointer to the image id +pointer id # pointer to the image ids +pointer x # pointer to the x coordinate array +pointer y # pointer to the y coordinate array +pointer nap # pointer to the number of apertures array +pointer rap # pointer to the aperture radii array +pointer mag # pointer to the magnitude array +pointer merr # pointer to the magnitude error array +int nptr # pointer to the current data point +int sortimid # does data need to be sorted on imid or id +real magerr # the maximum magnitude error + +int stat, dbufsize, lbufsize +pointer sp, fname, image, imname, filterid, line, sym +real itime, airmass, otime + +bool streq() +int ph_anxtimage(), ph_astardata(), ph_agetimage(), stnsymbols() +pointer stfind(), stenter() + +begin + # Allocate working space. + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + # Allocate some initial buffer space. + if (nptr == 0) { + dbufsize = DEF_BUFSIZE + call malloc (imid, dbufsize, TY_INT) + call malloc (id, dbufsize, TY_INT) + call malloc (x, dbufsize, TY_REAL) + call malloc (y, dbufsize, TY_REAL) + call malloc (nap, dbufsize, TY_INT) + call malloc (rap, naperts * dbufsize, TY_REAL) + call malloc (mag, naperts * dbufsize, TY_REAL) + call malloc (merr, naperts * dbufsize, TY_REAL) + } + + # Initialize. + Memc[image] = EOS + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + # Find the first image. + lbufsize = 0 + line = NULL + stat = ph_anxtimage (fd, columns, Memc[image], line, lbufsize) + + # Read the data. + while (stat != EOF) { + + # Check to see if the image name is already in the symbol table + # and if not enter it. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + sym = stenter (imtable, Memc[image], LEN_IMT_STRUCT) + IMT_IMNO(sym) = stnsymbols (imtable, 0) + IMT_NENTRIES(sym) = 0 + IMT_RO(sym) = INDEFR + IMT_OFFSET(sym) = 0 + call strcpy ("INDEF", IMT_IFILTER(sym), SZ_FNAME) + IMT_ITIME(sym) = INDEFR + IMT_XAIRMASS(sym) = INDEFR + IMT_NXAIRMASS(sym) = INDEFR + IMT_OTIME(sym) = INDEFR + call strcpy (Memc[image], IMT_IMNAME(sym), SZ_FNAME) + } + + # Decode the airmass information. + call ph_aimdata (Memc[line], columns, Memc[filterid], SZ_FNAME, + itime, airmass, otime) + + # Enter the new filterid, itime, and airmass only if there are no + # previous entries for that image. + + if (streq (IMT_IFILTER(sym), "INDEF")) + call strcpy (Memc[filterid], IMT_IFILTER(sym), SZ_FNAME) + if (IS_INDEFR(IMT_ITIME(sym))) { + if (IS_INDEFR(itime)) + IMT_ITIME(sym) = 1.0 + else + IMT_ITIME(sym) = itime + } + if (IS_INDEFR(IMT_XAIRMASS(sym))) { + if (! IS_INDEFR(airmass)) { + IMT_XAIRMASS(sym) = airmass + IMT_NXAIRMASS(sym) = airmass - DEF_AIROFFSET + } else { + IMT_XAIRMASS(sym) = DEF_AIROFFSET + IMT_NXAIRMASS(sym) = 0.0 + call printf ( + " Warning: Setting airmass for image %s to %g\n") + call pargstr (IMT_IMNAME(sym)) + call pargr (DEF_AIROFFSET) + } + } + if (IS_INDEFR(IMT_OTIME(sym))) + IMT_OTIME(sym) = otime + + # Does the data need to be sorted on image or id. This should + # usually not be necessary for most DAOPHOT or APPHOT data + # but will be necessary if data from the same images is contained + # in more than 1 file. + + if (IMT_NENTRIES(sym) > 0) + sortimid = YES + else + IMT_OFFSET(sym) = nptr + 1 + + # Get the data. + repeat { + + # Decode x, y, magnitude and error. + if (ph_astardata (Memc[line], columns, Memr[x+nptr], + Memr[y+nptr], Memi[nap+nptr], Memr[rap+naperts*nptr], + Memr[mag+naperts*nptr], Memr[merr+naperts*nptr], + naperts, magerr) == OK) { + Memi[imid+nptr] = IMT_IMNO(sym) + IMT_NENTRIES(sym) = IMT_NENTRIES(sym) + 1 + if (IMT_NENTRIES(sym) == 1) + IMT_RO(sym) = 0.5 * Memr[rap] + Memi[id+nptr] = IMT_NENTRIES(sym) + nptr = nptr + 1 + } + + # Allocate more buffer space if necessary. + if (nptr >= dbufsize) { + dbufsize = dbufsize + DEF_BUFSIZE + call realloc (imid, dbufsize, TY_INT) + call realloc (id, dbufsize, TY_INT) + call realloc (x, dbufsize, TY_REAL) + call realloc (y, dbufsize, TY_REAL) + call realloc (nap, dbufsize, TY_INT) + call realloc (rap, naperts * dbufsize, TY_REAL) + call realloc (mag, naperts * dbufsize, TY_REAL) + call realloc (merr, naperts * dbufsize, TY_REAL) + } + + # Decode the next data. + stat = ph_agetimage (fd, columns, Memc[imname], line, lbufsize) + + } until ((stat == EOF) || (! streq (Memc[imname], Memc[image]))) + + call strcpy (Memc[imname], Memc[image], SZ_FNAME) + + } + + if (line != NULL) { + call mfree (line, TY_CHAR) + line = NULL + lbufsize = 0 + } + + call sfree (sp) + + return (nptr) +end diff --git a/noao/digiphot/photcal/mkobsfile/phimtable.x b/noao/digiphot/photcal/mkobsfile/phimtable.x new file mode 100644 index 00000000..dbbadb9e --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phimtable.x @@ -0,0 +1,971 @@ +include <fset.h> +include "../lib/obsfile.h" + +# PH_ESETIMTABLE -- Read in the image set name and image names from the +# standard input and initialize the table fields. + +int procedure ph_esetimtable (imtable, nfilters, verify) + +pointer imtable # pointer to the image symbol table +int nfilters # the number of filters in the filter set +int verify # verify the user input + +int setno, nsymbols, field, stat +pointer sp, name, image, str, sym +bool streq() +int scan(), nscan(), strlen() +pointer stenter() + +begin + # Allocate working space. + call smark (sp) + call salloc (name, SZ_FNAME, TY_CHAR) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Initialize the set and image counters. + setno = 0 + nsymbols = 0 + field = 1 + call printf ("\n") + + # Loop over the user defined image sets. + repeat { + + if (field == 1) { + + # Print the prompt for the image set name. + call printf ( + "Enter name of image set %d (name, <EOF>=quit entry): ") + call pargi (setno + 1) + call flush (STDOUT) + + # Read the image set name. + stat = scan() + if (stat == EOF) { + call printf ("\n") + break + } else { + call gargwrd (Memc[name], SZ_FNAME) + if (nscan() != 1) + next + } + + # Verify the image set name. + if (verify == YES) { + call printf (" Verify name (name, <CR>=%s): ") + call pargstr (Memc[name]) + call flush (STDOUT) + stat = scan() + if (stat == EOF) + call printf ("\n") + else { + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() == 1) + call strcpy (Memc[str], Memc[name], SZ_FNAME) + } + } + + # Prepare for the next image set field. + field = field + 1 + setno = setno + 1 + next + + } else if (field <= (nfilters + 1)) { + + # Prompt for the image name. + call printf ( + " Enter image name %d (name, <CR>=INDEF): ") + call pargi (field - 1) + call flush (STDOUT) + + # Get the next image name. + stat = scan () + if (stat == EOF) { + call printf ("\n") + next + } else { + call gargwrd (Memc[image], SZ_FNAME) + if (nscan() != 1) + call strcpy ("INDEF", Memc[image], SZ_FNAME) + } + + # Verify + if (verify == YES) { + call printf ( + " Verify name (name, <CR>=%s): ") + call pargstr (Memc[image]) + call flush (STDOUT) + stat = scan () + if (stat == EOF) + call printf ("\n") + else { + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() == 1) + call strcpy (Memc[str], Memc[image], SZ_FNAME) + } + } + + nsymbols = nsymbols + 1 + field = field + 1 + + # Enter the new symbol name. + if (streq (Memc[image], "INDEF")) { + call sprintf (Memc[image+strlen(Memc[image])], SZ_FNAME, + "%d") + call pargi (nsymbols) + } + sym = stenter (imtable, Memc[image], LEN_IMT_STRUCT) + + IMT_IMSETNO(sym) = setno + IMT_IMNO(sym) = nsymbols + IMT_OFFSET(sym) = 0 + IMT_NENTRIES(sym) = 0 + + IMT_XSHIFT(sym) = 0.0 + IMT_YSHIFT(sym) = 0.0 + IMT_APERCOR(sym) = 0.0 + IMT_ITIME(sym) = INDEFR + IMT_OTIME(sym) = INDEFR + IMT_XAIRMASS(sym) = INDEFR + call strcpy ("INDEF", IMT_IFILTER(sym), SZ_FNAME) + + call strcpy (Memc[name], IMT_LABEL(sym), SZ_FNAME) + call strcpy (Memc[image], IMT_IMNAME(sym), SZ_FNAME) + + next + + } else { + field = 1 + } + } + + # Print a newline. + call printf ("\n") + + call sfree (sp) + + return (nsymbols) +end + + +# PH_SETIMTABLE -- Read in the image names from the image sets file +# and initialize the fields. + +int procedure ph_setimtable (imtable, fd, nfilters, verbose) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int nfilters # the number of filters in the set +int verbose # verbose mode + +char colon +int field, setno, nsymbols, nimages +pointer sp, name, image, sym +bool streq() +int fscan(), nscan(), strlen(), itoc() +pointer stenter() + +begin + call smark (sp) + call salloc (name, SZ_FNAME, TY_CHAR) + call salloc (image, SZ_FNAME, TY_CHAR) + + setno = 0 + nsymbols = 0 + + while (fscan (fd) != EOF) { + + # Initialize. + nimages = 0 + field = 1 + colon = ' ' + + # Decode the image set name and the image names for each set. + # Skip blank lines and lines beginning with a pound sign. + + repeat { + + if (field == 1) { + call gargwrd (Memc[name], SZ_FNAME) + if (Memc[name] == EOS || Memc[name] == '#') + break + field = field + 1 + if (Memc[name] == ':') { + Memc[name] = EOS + colon = ':' + } else + call gargc (colon) + if ((nscan() != field) || (colon != ':')) + break + field = field + 1 + setno = setno + 1 + if (Memc[name] == EOS) { + if (itoc (setno, Memc[name], SZ_FNAME) <= 0) + Memc[name] = EOS + } + next + } else { + call gargwrd (Memc[image], SZ_FNAME) + if (nscan() != field) + break + field = field + 1 + } + + nimages = nimages + 1 + nsymbols = nsymbols + 1 + + if (streq (Memc[image], "INDEF")) { + call sprintf (Memc[image+strlen(Memc[image])], SZ_FNAME, + "%d") + call pargi (nsymbols) + } + + # Enter the new symbol name. + sym = stenter (imtable, Memc[image], LEN_IMT_STRUCT) + + IMT_IMSETNO(sym) = setno + IMT_IMNO(sym) = nsymbols + IMT_OFFSET(sym) = 0 + IMT_NENTRIES(sym) = 0 + + IMT_XSHIFT(sym) = 0.0 + IMT_YSHIFT(sym) = 0.0 + IMT_APERCOR(sym) = 0.0 + IMT_ITIME(sym) = INDEFR + IMT_XAIRMASS(sym) = INDEFR + IMT_OTIME(sym) = INDEFR + call strcpy ("INDEF", IMT_IFILTER(sym), SZ_FNAME) + + call strcpy (Memc[name], IMT_LABEL(sym), SZ_FNAME) + call strcpy (Memc[image], IMT_IMNAME(sym), SZ_FNAME) + } + + # Test for error conditions. + if (verbose == NO) + next + if (Memc[name] == EOS || Memc[name] == '#') { + # blank lines and comment lines are skipped + } else if (colon != ':') { + call eprintf ("Warning: Error decoding image set %s\n") + call pargstr (Memc[name]) + } else if (nimages == 0) { + call eprintf ("Warning: Image set %d name %s is empty\n") + call pargi (setno) + call pargstr (Memc[name]) + } else if (nimages < nfilters) { + call eprintf ("Warning: Image set %d name %s is incomplete\n") + call pargi (setno) + call pargstr (Memc[name]) + } else if (nimages > nfilters) { + call eprintf ("Warning: Image set %d name %s ") + call pargi (setno) + call pargstr (Memc[name]) + call eprintf ("has more than %d images\n") + call pargi (nimages) + } + } + + call sfree (sp) + + return (nsymbols) +end + + +# PH_ESHIMTABLE -- Enter the correct values of the x-y coordinate shifts +# the image symbol table from the standard input. + +procedure ph_eshimtable (imtable, nimages, verify) + +pointer imtable # pointer to the symbol table +int nimages # number of images +int verify # verify the user input + +char comma +int i, osetno, setno +pointer sp, sym, symbol +real xshift, yshift +int scan(), nscan() +pointer sthead(), stnext() + +begin + if (nimages <= 0) + return + + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Reverse the order of the symbols in the symbol table. + + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + call printf ("\n") + + # Loop over the images. + osetno = 0 + do i = 1, nimages { + + symbol = Memi[sym+i-1] + setno = IMT_IMSETNO(symbol) + + # Issue the prompt for each set. + if (setno != osetno) { + call printf ("Image set %d (%s): ") + call pargi (setno) + call pargstr (IMT_LABEL(symbol)) + call printf ("Enter the shift in x and y\n") + } + + # Issue the prompt for each image. + call printf ( + " Image %s (xshift yshift, <CR>=0.0 0.0, <EOF>=quit entry): " ) + call pargstr (IMT_IMNAME(symbol)) + call flush (STDOUT) + + # Initialize + osetno = setno + + # Read a line from STDIN. + if (scan() == EOF) { + call printf ("\n") + break + } + + # Decode the x and y shifts and update the symbol table. + call gargr (xshift) + call gargr (yshift) + if (nscan() < 1) { + xshift = 0.0 + yshift = 0.0 + } else if (nscan() < 2) { + call reset_scan() + call gargr (xshift) + call gargc (comma) + call gargr (yshift) + if (nscan() < 3 || comma != ',') + yshift = 0.0 + } + IMT_XSHIFT(symbol) = xshift + IMT_YSHIFT(symbol) = yshift + + # Verify the input. + if (verify == NO) + next + + # Issue the verify prompt for each image. + call printf ( + " Verify (xshift yshift, <CR>=%g %g): " ) + call pargr (IMT_XSHIFT(symbol)) + call pargr (IMT_YSHIFT(symbol)) + call flush (STDOUT) + + # Read a line from STDIN. + if (scan() == EOF) { + call printf ("\n") + next + } + + # Decode the x and y shifts. + call gargr (xshift) + call gargr (yshift) + if (nscan() == 2) { + IMT_XSHIFT(symbol) = xshift + IMT_YSHIFT(symbol) = yshift + } else if (nscan() == 1) { + call reset_scan() + call gargr (xshift) + call gargc (comma) + call gargr (yshift) + IMT_XSHIFT(symbol) = xshift + if (nscan() == 3 && comma == ',') + IMT_YSHIFT(symbol) = yshift + } + } + + call printf ("\n") + + call sfree (sp) +end + + +# PH_SHIMTABLE -- Read in the x and y shifts file and enter the correct +# x and y shifts in the symbol table. + +procedure ph_shimtable (imtable, fd, verbose) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int verbose # verbose mode + +pointer sp, image, fname, sym +int fscan(), nscan() +pointer stfind() +real xoffset, yoffset + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + while (fscan (fd) != EOF) { + + # Get the image name and the x and y shifts data. + call gargwrd (Memc[image], SZ_FNAME) + call gargr (xoffset) + call gargr (yoffset) + if (nscan () != 3) + next + + # Locate the image name in the symbol table and enter the shifts. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + if (verbose == YES) { + call eprintf ( + "Warning: File: %s image: %s is ") + call pargstr (Memc[fname]) + call pargstr (Memc[image]) + call eprintf ("not in the image sets file\n") + } + next + } + + IMT_XSHIFT(sym) = xoffset + IMT_YSHIFT(sym) = yoffset + } + + call sfree (sp) +end + + +# PH_EAPIMTABLE -- Enter the correct values of the aperture correction in +# the image symbol table from the standard input. + +procedure ph_eapimtable (imtable, nimages, verify) + +pointer imtable # pointer to the symbol table +int nimages # number of images +int verify # verify the user input + +int i, osetno, setno +pointer sp, sym, symbol +real apercor +int scan(), nscan() +pointer sthead(), stnext() + +begin + if (nimages <= 0) + return + + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Reverse the order of the symbols in the symbol table. + + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + call printf ("\n") + + # Loop over the images. + osetno = 0 + do i = 1, nimages { + + symbol = Memi[sym+i-1] + setno = IMT_IMSETNO(symbol) + + # Issue the prompt for each set. + if (setno != osetno) { + call printf ("Image set %d (%s): ") + call pargi (setno) + call pargstr (IMT_LABEL(symbol)) + call printf ("Enter the aperture correction\n") + } + + # Issue the prompt for each image. + call printf ( + " Image %s (magnitude, <CR>=0.0, <EOF>=quit entry): " ) + call pargstr (IMT_IMNAME(symbol)) + call flush (STDOUT) + + # Initialize + osetno = setno + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + break + } + + # Decode the aperture correction and update the symbl table. + call gargr (apercor) + if (nscan() < 1) + apercor = 0.0 + IMT_APERCOR(symbol) = apercor + + # Optionally verify the input. + if (verify == NO) + next + + # Issue the verify prompt. + call printf ( + " Verify (magnitude, <CR>=%g): " ) + call pargr (IMT_APERCOR(symbol)) + call flush (STDOUT) + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + next + } + + # Decode the aperture correction. + call gargr (apercor) + if (nscan() == 1) + IMT_APERCOR(symbol) = apercor + } + + call printf ("\n") + + call sfree (sp) +end + + +# PH_APIMTABLE -- Read in the aperture corrections file and enter the correct +# aperture correction in the symbol table. + +procedure ph_apimtable (imtable, fd, verbose) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int verbose # print status, warning and error messages + +pointer sp, image, fname, sym +int fscan(), nscan() +pointer stfind() +real apercor + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + while (fscan (fd) != EOF) { + + # Get the image name and the x and y shifts data. + call gargwrd (Memc[image], SZ_FNAME) + call gargr (apercor) + if (nscan () != 2) + next + + # Locate the image name in the symbol table and enter the shifts. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + if (verbose == YES) { + call eprintf ( + "Warning: File: %s image: %s ") + call pargstr (Memc[fname]) + call pargstr (Memc[image]) + call eprintf ("is not in the image sets file\n") + } + next + } + + IMT_APERCOR(sym) = apercor + } + + call sfree (sp) +end + + +# PH_EOBSIMTABLE -- Enter the correct values of the filterid, exposure time and +# and airmass into the image symbol table from the standard input. + +procedure ph_eobsimtable (imtable, nimages, verify) + +pointer imtable # pointer to the symbol table +int nimages # number of images +int verify # verify the user input + +int i, osetno, setno +pointer sp, sym, filterid, symbol +real itime, xairmass, otime +int scan(), nscan(), strncmp() +pointer sthead(), stnext() + +begin + if (nimages <= 0) + return + + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + # Reverse the order of the symbols in the symbol table. + + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + call printf ("\n") + + # Loop over the images. + osetno = 0 + do i = 1, nimages { + + symbol = Memi[sym+i-1] + setno = IMT_IMSETNO(symbol) + + # Issue the prompt for each set. + if (setno != osetno) { + call printf ("Image set %d (%s): ") + call pargi (setno) + call pargstr (IMT_LABEL(symbol)) + call printf ( + "Enter filter id, exposure time, airmass, observation time\n") + } + + # Issue the prompt for each image. + call printf ( + " Image %s (f t X T, <CR>=4 X (INDEF), <EOF>=quit entry): " ) + call pargstr (IMT_IMNAME(symbol)) + call flush (STDOUT) + + # Initialize. + osetno = setno + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + break + } + + # Read in the filterid, exposure time and airmass. + call gargwrd (Memc[filterid], SZ_FNAME) + call gargr (itime) + call gargr (xairmass) + call gargr (otime) + if (nscan() < 1) { + call strcpy ("INDEF", Memc[filterid], SZ_FNAME) + itime = INDEFR + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 2) { + itime = INDEFR + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 3) { + xairmass = INDEFR + otime = INDEFR + } else if (nscan() < 4) { + otime = INDEFR + } + + # Update the symbol table. + if (strncmp (IMT_IFILTER(symbol), "INDEF", 5) == 0) + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + if (! IS_INDEFR(itime)) + IMT_ITIME(symbol) = itime + if (! IS_INDEFR(xairmass)) + IMT_XAIRMASS(symbol) = xairmass + if (! IS_INDEFR(otime)) + IMT_OTIME(symbol) = otime + + # Verify the input. + if (verify == NO) + next + + # Issue the verify prompt. + call printf ( + " Verify (f t X T, <CR>=%s %g %g %0.1h): ") + call pargstr (IMT_IFILTER(symbol)) + call pargr (IMT_ITIME(symbol)) + call pargr (IMT_XAIRMASS(symbol)) + call pargr (IMT_OTIME(symbol)) + call flush (STDOUT) + + # Scan the standard input. + if (scan() == EOF) { + call printf ("\n") + next + } + + # Read in the filterid, exposure time and airmass. + call gargwrd (Memc[filterid], SZ_FNAME) + call gargr (itime) + call gargr (xairmass) + call gargr (otime) + if (nscan() == 4) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + IMT_XAIRMASS(symbol) = xairmass + IMT_OTIME(symbol) = otime + } else if (nscan() == 3) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + IMT_XAIRMASS(symbol) = xairmass + } else if (nscan() == 2) { + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + IMT_ITIME(symbol) = itime + } else if (nscan() == 1) + call strcpy (Memc[filterid], IMT_IFILTER(symbol), SZ_FNAME) + } + + call printf ("\n") + + call sfree (sp) +end + + +# PH_OBSIMTABLE -- Enter the correct values of the exposure time, airmass and +# and filter id into the image table. The exposure time entered into the +# image table is the effective exposure time required to convert the +# computed magnitudes to the new exposure time. + +procedure ph_obsimtable (imtable, fd, columns, verbose) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +int verbose # print status, warning and error messages + +int stat, lbufsize +pointer sp, image, fname, filterid, line, sym +real itime, airmass, otime +bool streq() +int ph_nxtimage() +pointer stfind() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + Memc[image] = EOS + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + line = NULL + lbufsize = 0 + repeat { + + # Get the image name. + stat = ph_nxtimage (fd, columns, Memc[image], line, lbufsize) + if (stat == EOF) + break + + # Locate the image name in the symbol table. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + if (verbose == YES) { + call eprintf ( + "Warning: File: %s image: %s ") + call pargstr (Memc[fname]) + call pargstr (Memc[image]) + call eprintf ("is not in the image sets file\n") + } + next + } + + # Decode the data. + call ph_obsdata (Memc[line], columns, Memc[filterid], SZ_FNAME, + itime, airmass, otime) + + # Enter the data in the symbol table. + if (! IS_INDEFR(itime)) + IMT_ITIME(sym) = itime + if (! IS_INDEFR(airmass)) + IMT_XAIRMASS(sym) = airmass + if (! IS_INDEFR(otime)) + IMT_OTIME(sym) = otime + if (! streq (Memc[filterid], "INDEF")) + call strcpy (Memc[filterid], IMT_IFILTER(sym), SZ_FNAME) + } + + if (line != NULL) { + call mfree (line, TY_CHAR) + line = NULL + lbufsize = 0 + } + + call sfree (sp) +end + + +# PH_MKIMTABLE -- Create the image table using the input text file(s) +# and the list of column numbers defining the fields required by the +# program. + +int procedure ph_mkimtable (imtable, fd, columns, objid, x, y, mag, merr, + imid, id, nptr, normtime, sortimid, verbose) + +pointer imtable # pointer to the symbol table +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +pointer objid # pointer to the id array +pointer x # pointer to the x coordinate array +pointer y # pointer to the y coordinate array +pointer mag # pointer to the magnitude array +pointer merr # pointer to the magnitude error array +pointer imid # pointer to the image id +pointer id # pointer to the image ids +int nptr # pointer to the current data point +int normtime # normalize exposure times +int sortimid # does data need to be sorted on imid or id +int verbose # print status, warning and error messages + +int idcol, stat, dbufsize, lbufsize +pointer sp, fname, image, imname, filterid, objname, line, sym +real itime, airmass, otime + +bool streq() +int ph_nxtimage(), ph_stardata(), ph_getimage() +pointer stfind() + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (filterid, SZ_FNAME, TY_CHAR) + call salloc (objname, DEF_LENLABEL, TY_CHAR) + + # Is there id information. + idcol = columns[CAT_ID] + + # Allocate some initial buffer space. + if (nptr == 0) { + dbufsize = DEF_BUFSIZE + if (idcol > 0) + call malloc (objid, dbufsize * (DEF_LENLABEL+ 1), TY_CHAR) + call malloc (x, dbufsize, TY_REAL) + call malloc (y, dbufsize, TY_REAL) + call malloc (mag, dbufsize, TY_REAL) + call malloc (merr, dbufsize, TY_REAL) + call malloc (imid, dbufsize, TY_INT) + call malloc (id, dbufsize, TY_INT) + } + + # Initialize. + Memc[image] = EOS + call fstats (fd, F_FILENAME, Memc[fname], SZ_FNAME) + + # Find the first image. + lbufsize = 0 + line = NULL + stat = ph_nxtimage (fd, columns, Memc[image], line, lbufsize) + + # Read the data. + while (stat != EOF) { + + # Locate the image name in the symbol table. + sym = stfind (imtable, Memc[image]) + if (sym == NULL) { + if (Memc[image] != EOS && verbose == YES) { + call eprintf ( + "Warning: File: %s Image: %s ") + call pargstr (Memc[fname]) + call pargstr (Memc[image]) + call eprintf ("is not in the image sets file\n") + } + stat = ph_nxtimage (fd, columns, Memc[image], line, lbufsize) + next + } + + # Decode the filter, exposure time and airmass information. + call ph_imdata (Memc[line], columns, Memc[filterid], SZ_FNAME, + itime, airmass, otime) + + # Enter the new values in the symbol table. + if (IS_INDEFR(IMT_ITIME(sym))) { + if (IS_INDEFR(itime)) + IMT_ITIME(sym) = 1.0 + else if (normtime == NO) + IMT_ITIME(sym) = 1.0 + else + IMT_ITIME(sym) = itime + } else if (IMT_NENTRIES(sym) > 0) { + # do nothing + } else if (IS_INDEFR(itime)) { + if (normtime == NO) + IMT_ITIME(sym) = 1.0 + } else { + IMT_ITIME(sym) = IMT_ITIME(sym) / itime + } + if (IS_INDEFR(IMT_XAIRMASS(sym))) + IMT_XAIRMASS(sym) = airmass + if (IS_INDEFR(IMT_OTIME(sym))) + IMT_OTIME(sym) = otime + if (streq (IMT_IFILTER(sym), "INDEF")) + call strcpy (Memc[filterid], IMT_IFILTER(sym), SZ_FNAME) + + # Does the data need to be sorted on image or id. This should + # usually not be necessary for most DAOPHOT or APPHOT data + # but will be necessary if data from the same images is contained + # in more than 1 file. + + if (IMT_NENTRIES(sym) > 0) + sortimid = YES + else + IMT_OFFSET(sym) = nptr + 1 + + # Get the data. + repeat { + + # Decode x, y, magnitude and error. + if (ph_stardata (Memc[line], columns, Memc[objname], + Memr[x+nptr], Memr[y+nptr], Memr[mag+nptr], + Memr[merr+nptr]) == OK) { + if (idcol > 0) + call strcpy (Memc[objname], + Memc[objid+nptr*(DEF_LENLABEL+1)], + DEF_LENLABEL) + Memi[imid+nptr] = IMT_IMNO(sym) + IMT_NENTRIES(sym) = IMT_NENTRIES(sym) + 1 + Memi[id+nptr] = IMT_NENTRIES(sym) + nptr = nptr + 1 + } + + # Allocate more buffer space if necessary. + if (nptr >= dbufsize) { + dbufsize = dbufsize + DEF_BUFSIZE + if (idcol > 0) + call realloc (objid, dbufsize * (DEF_LENLABEL + 1), + TY_CHAR) + call realloc (x, dbufsize, TY_REAL) + call realloc (y, dbufsize, TY_REAL) + call realloc (mag, dbufsize, TY_REAL) + call realloc (merr, dbufsize, TY_REAL) + call realloc (imid, dbufsize, TY_INT) + call realloc (id, dbufsize, TY_INT) + } + + # Decode the next data. + stat = ph_getimage (fd, columns, Memc[imname], line, lbufsize) + + } until ((stat == EOF) || (! streq (Memc[imname], Memc[image]))) + + call strcpy (Memc[imname], Memc[image], SZ_FNAME) + + } + + if (line != NULL) { + call mfree (line, TY_CHAR) + line = NULL + lbufsize = 0 + } + + call sfree (sp) + + return (nptr) +end diff --git a/noao/digiphot/photcal/mkobsfile/phmatch.x b/noao/digiphot/photcal/mkobsfile/phmatch.x new file mode 100644 index 00000000..59fadb60 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phmatch.x @@ -0,0 +1,487 @@ +include <mach.h> +include "../lib/obsfile.h" + +# PH_JOIN -- For each individual field join the data from each image +# line for line. + +procedure ph_join (out, label, sym, filters, nfilters, objid, x, y, mag, merr, + allfilters, wrap, verbose) + +int out # the output file descriptor +char label[ARB] # the label string +pointer sym[ARB] # list of pointers to the image set data +char filters[ARB] # the filter id string +int nfilters # number of filters in an image set +char objid[DEF_LENLABEL,ARB] # array of object ids +real x[ARB] # array of shifted x coordinates +real y[ARB] # array of shifted y coordinates +real mag[ARB] # array of corrected magnitudes +real merr[ARB] # array of magnitude errors +int allfilters # match objects in all filters ? +int wrap # format the output file for easy reading ? +int verbose # print status, warning and error messages ? + +int i, j, nobjs, dptr, len_label, useid +pointer sp, outfilter, idlabel, newlabel, bptrs, eptrs +bool streq() +int ph_nthlabel(), strlen() + +begin + # Find the maximum length of the output. Return if there is no data + # in any of the frames, or if there is no data in one frame and + # the match in all filters switch is on. + + if (allfilters == YES) { + nobjs = MAX_INT + do i = 1, nfilters { + if (sym[i] == NULL) { + nobjs = 0 + break + } + if (IMT_NENTRIES(sym[i]) == 0) { + nobjs = 0 + break + } + nobjs = min (nobjs, IMT_NENTRIES(sym[i])) + } + } else { + nobjs = 0 + do i = 1, nfilters { + if (sym[i] == NULL) + next + nobjs = max (nobjs, IMT_NENTRIES(sym[i])) + } + } + + + # Return if there is no data. + if (nobjs <= 0) { + if (verbose == YES) { + call printf ( + "\tImage set: %s 0 stars written to the observations file\n") + call pargstr (label) + } + return + } + + # Allocate working space. + call smark (sp) + call salloc (idlabel, SZ_FNAME, TY_CHAR) + call salloc (newlabel, SZ_FNAME, TY_CHAR) + call salloc (outfilter, SZ_FNAME, TY_CHAR) + call salloc (bptrs, nfilters, TY_INT) + call salloc (eptrs, nfilters, TY_INT) + + # Determine whether or not to use the object id as the object label + # and determine the length of the label. + + if (objid[1,1] == EOS) { + useid = NO + if (nobjs == 1) + len_label = max (DEF_LENLABEL, strlen (label)) + else + len_label = max (DEF_LENLABEL, strlen (label)+DEF_LENINDEX) + } else { + useid = YES + len_label = DEF_LENLABEL + } + + # Initialize the arrays which point to the beginning and ending of + # the data for each image in the set. + + do j = 1, nfilters { + if (sym[j] == NULL) { + Memi[bptrs+j-1] = 0 + Memi[eptrs+j-1] = 0 + } else { + Memi[bptrs+j-1] = IMT_OFFSET(sym[j]) + Memi[eptrs+j-1] = IMT_OFFSET(sym[j]) + IMT_NENTRIES(sym[j]) - 1 + } + } + + # Write out the output. + do i = 1, nobjs { + + # Determine the object id. This is the object id of the first + # defined filter. + if (useid == YES) { + Memc[idlabel] = EOS + do j = 1, nfilters { + if (Memi[bptrs+j-1] == 0) + next + call strcpy (objid[1,Memi[bptrs+j-1]], Memc[idlabel], + DEF_LENLABEL) + break + } + } + + do j = 1, nfilters { + + # Write the label. + if (j == 1) { + if (useid == YES && Memc[idlabel] != EOS) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (Memc[idlabel]) + } else if (nobjs == 1) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (label) + } else { + call sprintf (Memc[newlabel], SZ_FNAME, "%s-%d") + call pargstr (label) + call pargi (i) + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (Memc[newlabel]) + } + } else if (wrap == YES) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr ("*") + } else { + call fprintf (out, " ") + } + + # Write the results. + if (ph_nthlabel (filters, j, Memc[outfilter], SZ_FNAME) != j) + Memc[outfilter] = EOS + dptr = Memi[bptrs+j-1] + if (dptr == 0) { + call fprintf (out, + "%-10.10s %11.1h %6.3f %9.3f %9.3f %7.3f %6.3f") + call pargstr (Memc[outfilter]) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + } else { + if (IMT_OTIME(sym[j]) >= 0.0 && IMT_OTIME(sym[j]) <= 24.0) + call fprintf (out, + "%-10.10s %11.1h %6.3f %9.3f %9.3f %7.3f %6.3f") + else + call fprintf (out, + "%-10.10s %11g %6.3f %9.3f %9.3f %7.3f %6.3f") + if (streq (IMT_IFILTER(sym[j]), "INDEF")) + call pargstr (Memc[outfilter]) + else + call pargstr (IMT_IFILTER(sym[j])) + call pargr (IMT_OTIME(sym[j])) + call pargr (IMT_XAIRMASS(sym[j])) + call pargr (x[dptr] - IMT_XSHIFT(sym[j])) + call pargr (y[dptr] - IMT_YSHIFT(sym[j])) + call pargr (mag[dptr]) + call pargr (merr[dptr]) + } + if (wrap == YES) + call fprintf (out, "\n") + else if (j == nfilters) + call fprintf (out, "\n") + + # Increment the data pointers making sure to check for non- + # existent data and unequal length object lists. + + if (dptr == 0) + next + if (dptr >= Memi[eptrs+j-1]) + dptr = 0 + else + Memi[bptrs+j-1] = dptr + 1 + } + } + + if (verbose == YES) { + call printf ( + "\tImage set: %s %d stars written to the observations file\n") + call pargstr (label) + call pargi (nobjs) + } + + call sfree (sp) +end + + +# PH_MERGE -- For an individual field join the data from each image (filter) in +# the image (filter) set by matching the x-y positions. + +procedure ph_merge (out, label, sym, filters, nfilters, objid, x, y, mag, merr, + ysort, match, index, ndata, tolerance, allfilters, wrap, verbose) + +int out # the output file descriptor +char label[ARB] # label string +pointer sym[ARB] # list of pointers to the image set data +char filters[ARB] # filter string +int nfilters # number of images an image set +char objid[DEF_LENLABEL,ARB] # array of object ids +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real merr[ARB] # array of magnitude errors +int ysort[ARB] # the y sort index +int match[ARB] # the matching array +int index[ndata,ARB] # the matching index array +int ndata # the number of data points +real tolerance # tolerance in pixels +int allfilters # match objects in all filters +int wrap # format the output file for easy reading ? +int verbose # print status, warning and error messages + +int i, j, k, l, biptr, eiptr, bjptr, ejptr, len_label +int ntimes, nframe, nobjs, nmatch, starok, lsort, useid +pointer sp, idlabel, newlabel, outfilter +real tol2, dx, dy, maxrsq, distsq +bool streq() +int ph_nthlabel(), strlen() + +begin + # Allocate working space. + call smark (sp) + call salloc (idlabel, SZ_FNAME, TY_CHAR) + call salloc (newlabel, SZ_FNAME, TY_CHAR) + call salloc (outfilter, SZ_FNAME, TY_CHAR) + + # Initialize. + nframe = 0 + do i = 1, nfilters { + if (sym[i] != NULL) { + call amovki (NO, match[IMT_OFFSET(sym[i])], + IMT_NENTRIES(sym[i])) + if (IMT_NENTRIES(sym[i]) > 0) + nframe = nframe + 1 + } + } + tol2 = tolerance ** 2 + nobjs = 0 + if (allfilters == YES) + ntimes = 1 + else + ntimes = nfilters + + # Loop over the filters. + do i = 1, ntimes { + + # Find the data pointers for the first frame. + if (sym[i] == NULL) + next + biptr = IMT_OFFSET(sym[i]) + if (biptr == 0) + next + eiptr = biptr + IMT_NENTRIES(sym[i]) - 1 + + # Initialize. + nmatch = 0 + + # Clear the index array. + do j = 1, nfilters + call aclri (index[1,j], ndata) + + # Loop over the other frames matching stars. + do j = i + 1, nfilters { + + # Find the data pointers for the second frame. + if (sym[j] == NULL) + next + bjptr = IMT_OFFSET(sym[j]) + if (bjptr == 0) + next + ejptr = bjptr + IMT_NENTRIES(sym[j]) - 1 + + # Loop over the stars in the first frame. + nmatch = 0 + do k = biptr, eiptr { + + if (match[ysort[k]] == YES) + next + if (IS_INDEFR(x[ysort[k]]) || IS_INDEFR(y[ysort[k]])) + next + nmatch = nmatch + 1 + index[ysort[k]-biptr+1,i] = ysort[k] + maxrsq = tol2 + lsort = 0 + + do l = bjptr, ejptr { + if (match[ysort[l]] == YES) + next + if (IS_INDEFR(x[ysort[l]]) || IS_INDEFR(y[ysort[l]])) + next + dy = y[ysort[l]] - y[ysort[k]] + if (dy > tolerance) + break + dx = x[ysort[l]] - x[ysort[k]] + distsq = dx * dx + dy * dy + if (distsq > maxrsq) + next + if (lsort > 0) + match[ysort[lsort]] = NO + match[ysort[l]] = YES + index[ysort[k]-biptr+1,j] = ysort[l] + maxrsq = distsq + lsort = l + } + } + + # If all the stars have been matched quit the loop. + if (nmatch == 0) + break + } + + # Search for unmatched stars in the last frame. + if ((nframe == 1) || (allfilters == NO && i == nfilters)) { + do k = biptr, eiptr { + if (match[ysort[k]] == YES) + next + if (IS_INDEFR(x[ysort[k]]) || IS_INDEFR(y[ysort[k]])) + next + nmatch = nmatch + 1 + index[ysort[k]-biptr+1,i] = ysort[k] + } + } + + # Use the individual star id. + if (objid[1,1] == EOS) + useid = NO + else + useid = YES + + # Write the results. + do j = 1, eiptr - biptr + 1 { + + # Break if no stars were matched. + if (nmatch <= 0) + break + #next + + # If the allfilters switch is on only print points with + # data for all the filters, otherwise print all the objects. + + if (allfilters == YES) { + starok = YES + do k = 1, nfilters { + if (index[j,k] > 0) + next + starok = NO + break + } + } else { + starok = NO + do k = 1, nfilters { + if (index[j,k] <= 0) + next + starok = YES + break + } + } + if (starok == NO) + next + + if (useid == YES) + len_label = DEF_LENLABEL + else if (nobjs == 1) + len_label = max (DEF_LENLABEL, strlen (label)) + else + len_label = max (DEF_LENLABEL, strlen (label)+DEF_LENINDEX) + + # Find the label. + if (useid == YES) { + Memc[idlabel] = EOS + do k = 1, nfilters { + if (index[j,k] == 0) + next + call strcpy (objid[1,index[j,k]], Memc[idlabel], + DEF_LENLABEL) + break + } + } + + # Loop over the filters, writing the data for each point. + do k = 1, nfilters { + + # Write the label. + if (k == 1) { + if (useid == YES && Memc[idlabel] != EOS) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (Memc[idlabel]) + } else if ((nmatch == 1) && (nobjs == 0)) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (label) + } else { + call sprintf (Memc[newlabel], SZ_FNAME, "%s-%d") + call pargstr (label) + call pargi (nobjs + 1) + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr (Memc[newlabel]) + } + } else if (wrap == YES) { + call fprintf (out, "%*.*s ") + call pargi (-len_label) + call pargi (len_label) + call pargstr ("*") + } else { + call fprintf (out, " ") + } + + # Write the data. + if (ph_nthlabel (filters, k, Memc[outfilter], + SZ_FNAME) != k) + Memc[outfilter] = EOS + if (index[j,k] == 0) { + call fprintf (out, + "%-10.10s %11.1h %6.3f %9.3f %9.3f %7.3f %6.3f") + call pargstr (Memc[outfilter]) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + } else { + if (IMT_OTIME(sym[k]) >= 0.0 && + IMT_OTIME(sym[k]) <= 24.0) + call fprintf (out, + "%-10.10s %11.1h %6.3f %9.3f %9.3f %7.3f %6.3f") + else + call fprintf (out, + "%-10.10s %11g %6.3f %9.3f %9.3f %7.3f %6.3f") + if (streq (IMT_IFILTER(sym[k]), "INDEF")) + call pargstr (Memc[outfilter]) + else + call pargstr (IMT_IFILTER(sym[k])) + call pargr (IMT_OTIME(sym[k])) + call pargr (IMT_XAIRMASS(sym[k])) + call pargr (x[index[j,k]] - IMT_XSHIFT(sym[k])) + call pargr (y[index[j,k]] - IMT_YSHIFT(sym[k])) + call pargr (mag[index[j,k]]) + call pargr (merr[index[j,k]]) + } + if (wrap == YES) + call fprintf (out, "\n") + else if (k == nfilters) + call fprintf (out, "\n") + + } + + nobjs = nobjs + 1 + } + } + + if (verbose == YES) { + call printf ( + "\tImage set: %s %d stars written to the observations file\n") + call pargstr (label) + call pargi (nobjs) + } + + call sfree (sp) +end diff --git a/noao/digiphot/photcal/mkobsfile/phsort.x b/noao/digiphot/photcal/mkobsfile/phsort.x new file mode 100644 index 00000000..a29920b0 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/phsort.x @@ -0,0 +1,528 @@ +include "../lib/obsfile.h" + +define LOGPTR 32 # log2(maxpts) (4e9) + +# PH_1C4R2ISORT -- Vector quicksort on the first and second indices arrays, +# where the second index is used to resolve ambiguities in the first index. +# An additional 4 input arrays are sorted as well. + +procedure ph_1c4r2isort (label, d1, d2, d3, d4, findex, sindex, npix) + +char label[DEF_LENLABEL,ARB] # the char array +real d1[ARB] # the first input data array +real d2[ARB] # the second input data array +real d3[ARB] # the third input data array +real d4[ARB] # the fourth input data array +int findex[ARB] # first index array which is sorted on +int sindex[ARB] # second index array which is sorted on +int npix # number of pixels + +real tempr +int fpivot, spivot, tempi +int i, j, k, p, lv[LOGPTR], uv[LOGPTR] +char tempc[DEF_LENLABEL] +int ph_2icompare() +define swapi {tempi=$1;$1=$2;$2=tempi} +define swapr {tempr=$1;$1=$2;$2=tempr} +define swapc {call strcpy ($1, tempc, DEF_LENLABEL);call strcpy ($2, $1, DEF_LENLABEL);call strcpy (tempc, $2, DEF_LENLABEL)} + +begin + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already sorted + # array. + + k = (lv[p] + uv[p]) / 2 + swapr (d1[j], d1[k]) + swapr (d2[j], d2[k]) + swapr (d3[j], d3[k]) + swapr (d4[j], d4[k]) + swapc (label[1,j], label[1,k]) + swapi (findex[j], findex[k]) + swapi (sindex[j], sindex[k]) + fpivot = findex[j] # pivot line + spivot = sindex[j] + + while (i < j) { + for (i=i+1; ph_2icompare (findex[i], sindex[i], fpivot, + spivot) < 0; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (ph_2icompare (findex[j], sindex[j], fpivot, + spivot) <= 0) + break + if (i < j) { # switch elements + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (d4[i], d4[j]) + swapc (label[1,i], label[1,j]) + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + } + } + + j = uv[p] # move pivot to position i + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (d4[i], d4[j]) + swapc (label[1,i], label[1,j]) + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } +end + + +# PH_4R2ISORT -- Vector quicksort on the first and second indices arrays, +# where the second index is used to resolve ambiguities in the first index. +# An additional 4 input arrays are sorted as well. + +procedure ph_4r2isort (d1, d2, d3, d4, findex, sindex, npix) + +real d1[ARB] # the first input data array +real d2[ARB] # the second input data array +real d3[ARB] # the third input data array +real d4[ARB] # the fourth input data array +int findex[ARB] # first index array which is sorted on +int sindex[ARB] # second index array which is sorted on +int npix # number of pixels + +int fpivot, spivot, tempi +int i, j, k, p, lv[LOGPTR], uv[LOGPTR] +real tempr +int ph_2icompare() +define swapi {tempi=$1;$1=$2;$2=tempi} +define swapr {tempr=$1;$1=$2;$2=tempr} + +begin + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already sorted + # array. + + k = (lv[p] + uv[p]) / 2 + swapr (d1[j], d1[k]) + swapr (d2[j], d2[k]) + swapr (d3[j], d3[k]) + swapr (d4[j], d4[k]) + swapi (findex[j], findex[k]) + swapi (sindex[j], sindex[k]) + fpivot = findex[j] # pivot line + spivot = sindex[j] + + while (i < j) { + for (i=i+1; ph_2icompare (findex[i], sindex[i], fpivot, + spivot) < 0; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (ph_2icompare (findex[j], sindex[j], fpivot, + spivot) <= 0) + break + if (i < j) { # switch elements + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (d4[i], d4[j]) + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + } + } + + j = uv[p] # move pivot to position i + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (d4[i], d4[j]) + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } +end + + +# PH_3RIRSORT -- Vector quicksort on the index and the rindex array, +# where the rindex array is used to resolve ambiguities in the index array. +# The three real arrays are sorted as well. + +procedure ph_3rirsort (d1, d2, d3, index, rindex, npix) + +real d1[ARB] # the first input data array +real d2[ARB] # the second input data array +real d3[ARB] # the third input data array +int index[ARB] # first index array which is sorted on +real rindex[ARB] # second index array which is sorted on +int npix # number of pixels + +int fpivot, tempi +int i, j, k, p, lv[LOGPTR], uv[LOGPTR] +real rpivot, tempr +int ph_ircompare() +define swapi {tempi=$1;$1=$2;$2=tempi} +define swapr {tempr=$1;$1=$2;$2=tempr} + +begin + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already sorted + # array. + + k = (lv[p] + uv[p]) / 2 + swapr (d1[j], d1[k]) + swapr (d2[j], d2[k]) + swapr (d3[j], d3[k]) + swapi (index[j], index[k]) + swapr (rindex[j], rindex[k]) + fpivot = index[j] # pivot line + rpivot = rindex[j] + + while (i < j) { + for (i=i+1; ph_ircompare (index[i], rindex[i], fpivot, + rpivot) < 0; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (ph_ircompare (index[j], rindex[j], fpivot, + rpivot) <= 0) + break + if (i < j) { # switch elements + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (rindex[i], rindex[j]) + swapi (index[i], index[j]) # interchange elements + } + } + + j = uv[p] # move pivot to position i + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + swapr (d3[i], d3[j]) + swapr (rindex[i], rindex[j]) + swapi (index[i], index[j]) # interchange elements + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } +end + + +# PH_QSORT -- Vector quicksort. In this version the index array is +# sorted. + +procedure ph_qsort (data, index, npix, offset) + +real data[ARB] # data array +int index[ARB] # index array +int npix # number of pixels +int offset # the index offset + +int i, j, lv[LOGPTR], p, uv[LOGPTR], temp +real pivot + +begin + # Initialize the indices for an inplace sort. + do i = 1, npix + index[i] = i + + p = 1 + lv[1] = 1 + uv[1] = npix + while (p > 0) { + + # If only one elem in subset pop stack otherwise pivot line. + if (lv[p] >= uv[p]) + p = p - 1 + else { + i = lv[p] - 1 + j = uv[p] + pivot = data[index[j]] + + while (i < j) { + for (i=i+1; data[index[i]] < pivot; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (data[index[j]] <= pivot) + break + if (i < j) { # out of order pair + temp = index[j] # interchange elements + index[j] = index[i] + index[i] = temp + } + } + + j = uv[p] # move pivot to position i + temp = index[j] # interchange elements + index[j] = index[i] + index[i] = temp + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + p = p + 1 # push onto stack + } + } + + do i = 1, npix + index[i] = index[i] + offset - 1 +end + + +# PH_2ICOMPARE -- Comparison routine for PH_4R2ISORT. + +int procedure ph_2icompare (findex, sindex, fpivot, spivot) + +int findex # the first index value +int sindex # the second index value +int fpivot # the first pivot value +int spivot # the second pivot value + +begin + if (findex < fpivot) + return (-1) + else if (findex > fpivot) + return (1) + else if (sindex < spivot) + return (-1) + else if (sindex > spivot) + return (1) + else + return (0) +end + + +# PH_IRCOMPARE -- Comparison routine for PH_3RIRSORT. + +int procedure ph_ircompare (findex, sindex, fpivot, spivot) + +int findex # the first index value +real sindex # the second index value +int fpivot # the first pivot value +real spivot # the second pivot value + +begin + if (findex < fpivot) + return (-1) + else if (findex > fpivot) + return (1) + else if (sindex < spivot) + return (-1) + else if (sindex > spivot) + return (1) + else + return (0) +end + + +# PH_5R3ISORT -- Vector quicksort on the first and second indices arrays, +# where the second index is used to resolve ambiguities in the first index. +# An additional 5 input arrays are sorted as well. + +procedure ph_5r3isort (findex, sindex, i1, d1, d2, d3, d4, d5, naperts, npix) + +int findex[ARB] # first index array which is sorted on +int sindex[ARB] # second index array which is sorted on +int i1[ARB] # the 3 integer array +real d1[ARB] # the first input data array +real d2[ARB] # the second input data array +real d3[naperts,ARB] # the third input data array +real d4[naperts,ARB] # the fourth input data array +real d5[naperts,ARB] # the fifth input data array +int naperts # number of apertures +int npix # number of pixels + +int fpivot, spivot, tempi +int i, j, k, l,p, lv[LOGPTR], uv[LOGPTR] +real tempr +int ph_2aicompare() +define swapi {tempi=$1;$1=$2;$2=tempi} +define swapr {tempr=$1;$1=$2;$2=tempr} + +begin + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already sorted + # array. + + k = (lv[p] + uv[p]) / 2 + swapi (findex[j], findex[k]) + swapi (sindex[j], sindex[k]) + swapi (i1[j], i1[k]) + swapr (d1[j], d1[k]) + swapr (d2[j], d2[k]) + do l = 1, naperts { + swapr (d3[l,j], d3[l,k]) + swapr (d4[l,j], d4[l,k]) + swapr (d5[l,j], d5[l,k]) + } + fpivot = findex[j] # pivot line + spivot = sindex[j] + + while (i < j) { + for (i=i+1; ph_2aicompare (findex[i], sindex[i], fpivot, + spivot) < 0; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (ph_2aicompare (findex[j], sindex[j], fpivot, + spivot) <= 0) + break + if (i < j) { # switch elements + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + swapi (i1[i], i1[j]) + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + do l = 1, naperts { + swapr (d3[l,i], d3[l,j]) + swapr (d4[l,i], d5[l,j]) + swapr (d5[l,i], d5[l,j]) + } + } + } + + j = uv[p] # move pivot to position i + swapi (sindex[i], sindex[j]) + swapi (findex[i], findex[j]) # interchange elements + swapi (i1[i], i1[j]) + swapr (d1[i], d1[j]) + swapr (d2[i], d2[j]) + do l = 1, naperts { + swapr (d3[l,i], d3[l,j]) + swapr (d4[l,i], d4[l,j]) + swapr (d5[l,i], d5[l,j]) + } + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } +end + + +# PH_2AICOMPARE -- Comparison routine for PH_5R3ISORT. + +int procedure ph_2aicompare (findex, sindex, fpivot, spivot) + +int findex # the first index value +int sindex # the second index value +int fpivot # the first pivot value +int spivot # the second pivot value + +begin + if (findex < fpivot) + return (-1) + else if (findex > fpivot) + return (1) + else if (sindex < spivot) + return (-1) + else if (sindex > spivot) + return (1) + else + return (0) +end diff --git a/noao/digiphot/photcal/mkobsfile/t_apfile.x b/noao/digiphot/photcal/mkobsfile/t_apfile.x new file mode 100644 index 00000000..59bfb7ec --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/t_apfile.x @@ -0,0 +1,727 @@ +include <fset.h> +include <ctotok.h> +include <ctype.h> +include "../lib/apfile.h" + +# T_APFILE -- Compute apperture corrections using the data extracted +# from a simple text file written by the user and an optional corrections +# file. APFILE expects to see the following nine quantities in the input data +# file(s): image name, xcenter, ycenter, filterid, exposure time, xairmass, +# list of apertures, list of magnitudes and list of magnitude errors, and +# uses a simple list of columns to determine which data is in which column. +# The obsparams file is decoded in the same way. + +procedure t_apfile () + +int photfiles # pointer to the list of photometry files +pointer columnstr # pointer to the list of columns +int naperts # number of apertures to extract +int mode # the file mode for the log and plot files +int smallap # the index of the smallest aperture to use +int largeap # the index of the largest aperture to use +real maglim # the maximum magnitude error +int mterms # the maximum number of terms in the cog model to fit +int interactive # interactive mode +pointer graphics # the default graphics device +int verify # verify interactive user input + +int incat, apcat, magfd, logfd, plotfd, obs, csp, cp, nincolumns, npts +pointer sp, column, incolumns, obscolumns, fname, params +pointer imtable, imid, id, x, y, nap, rap, mag, merr, sortimid, gd, mgd +bool clgetb() +int clpopnu(), clplen(), ph_agrange(), open(), btoi(), clgeti() +int ctoi(), clgfil(), ph_amkimtable(), fstati() +pointer stopen(), gopen() +real clgetr() + +begin + # Setup to flush on a newline. + if (fstati (STDOUT, F_REDIR) == NO) + call fseti (STDOUT, F_FLUSHNL, YES) + + # Allocate working memory. + call smark (sp) + call salloc (columnstr, SZ_LINE, TY_CHAR) + call salloc (column, SZ_FNAME, TY_CHAR) + call salloc (incolumns, CAT_NFIELDS, TY_INT) + call salloc (obscolumns, OBS_NFIELDS, TY_INT) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (params, MAX_MTERMS, TY_DOUBLE) + + # Get the input file list. + photfiles = clpopnu ("photfiles") + if (clplen (photfiles) <= 0) + call error (0, "The input file list is empty") + + # Decode the input columns list. There must be at least CAT_RAPERT + # range specifiers in the incolumns list. The number 0 can be used as + # a place holder for missing columns. + + call clgstr ("incolumns", Memc[columnstr], SZ_LINE) + call amovki (0, Memi[incolumns], CAT_NFIELDS) + nincolumns = 0 + csp = 1 + while (ph_agrange (Memc[columnstr], csp, Memc[column], SZ_FNAME) != + EOF) { + + # Decode the individual ranges. + if (nincolumns >= CAT_NFIELDS) + call error (0, + "Too many fields specified in the incolumns parameter") + cp = 1 + if (ctoi (Memc[column], cp, Memi[incolumns+nincolumns]) > 0) + nincolumns = nincolumns + 1 + else + break + } + + # Check that sufficient fields have been defined. + if (nincolumns < CAT_MERR) + call error (0, "Too few fields defined in the incolumns parameter") + + naperts = clgeti ("naperts") + smallap = max (1, min (naperts, clgeti ("smallap"))) + largeap = clgeti ("largeap") + if (IS_INDEFI(largeap) || largeap <= 0) + largeap = naperts + else + largeap = max (1, min (naperts, largeap)) + if (largeap <= smallap) + call error (0, + "The large aperture is less than or equal to the large aperture\n") + + # Open the output catalog file. + call clgstr ("apercors", Memc[fname], SZ_FNAME) + apcat = open (Memc[fname], NEW_FILE, TEXT_FILE) + + # Determine the access mode for the plot and log files. + if (clgetb ("append")) + mode = APPEND + else + mode = NEW_FILE + + # Open the best magnitudes file if any. + call clgstr ("magfile", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + magfd = NULL + else + magfd = open (Memc[fname], NEW_FILE, TEXT_FILE) + + # Open the logfile if any. + call clgstr ("logfile", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + logfd = NULL + else + logfd = open (Memc[fname], mode, TEXT_FILE) + + # Open the plot file if any. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("plotfile", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + plotfd = NULL + else + plotfd = open (Memc[fname], mode, BINARY_FILE) + if (plotfd != NULL) + mgd = gopen (Memc[graphics], NEW_FILE, plotfd) + else + mgd = NULL + + # Open the observing parameters file and decode the column list. + + call clgstr ("obsparams", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + obs = NULL + else { + obs = open (Memc[fname], READ_ONLY, TEXT_FILE) + call amovki (0, Memi[obscolumns], OBS_NFIELDS) + if (obs != STDIN) { + call clgstr ("obscolumns", Memc[columnstr], SZ_LINE) + csp = 1 + nincolumns = 1 + Memi[obscolumns] = 1 + while (ph_agrange (Memc[columnstr], csp, Memc[column], + SZ_FNAME) != EOF && (nincolumns < OBS_NFIELDS)) { + cp = 1 + if (ctoi (Memc[column], cp, + Memi[obscolumns+nincolumns]) <= 0) + Memi[obscolumns+nincolumns] = 0 + else if (Memi[obscolumns+nincolumns] == 1) + nincolumns = nincolumns - 1 + nincolumns = nincolumns + 1 + } + } + } + + # Verify interactive user input. + maglim = clgetr ("maglim") + mterms = max (1, min (MAX_MTERMS, clgeti ("nparams"))) + call ph_agetp (Memd[params]) + interactive = btoi (clgetb ("interactive")) + verify = btoi (clgetb ("verify")) + + # Open a symbol table to hold a list of image charactersitics. + + imtable = stopen ("imtable", 2 * LEN_IMTABLE, LEN_IMTABLE, 10 * + LEN_IMTABLE) + + # Read in the data. Extract each unique image name and accompanying + # airmass and add it to the symbol table. For each star in each image + # extract the x and y centers, and the list of aperture radii, + # magnitudes, and magnitude errors. + + imid = NULL + id = NULL + x = NULL + y = NULL + nap = NULL + rap = NULL + mag = NULL + merr = NULL + sortimid = NO + + npts = 0 + while (clgfil (photfiles, Memc[fname], SZ_FNAME) != EOF) { + incat = open (Memc[fname], READ_ONLY, TEXT_FILE) + npts = ph_amkimtable (imtable, incat, Memi[incolumns], naperts, + imid, id, x, y, nap, rap, mag, merr, npts, sortimid, maglim) + call close (incat) + } + + call realloc (imid, npts, TY_INT) + call realloc (id, npts, TY_INT) + call realloc (x, npts, TY_REAL) + call realloc (y, npts, TY_REAL) + call realloc (nap, npts, TY_INT) + call realloc (rap, naperts * npts, TY_REAL) + call realloc (mag, naperts * npts, TY_REAL) + call realloc (merr, naperts * npts, TY_REAL) + + # For each image referenced in the airmass file extract + # the image name, and the airmass if defined. Enter the new value of + # airmass in the symbol table and compute the effective + # exposure time required to correct the instrumental magnitudes to + # the new exposure time. + + if (obs != NULL) { + if (obs == STDIN) + call ph_eaobsimtable (imtable, verify) + else + call ph_aobsimtable (imtable, obs, Memi[obscolumns]) + } + + # Sort the data by image id and object id if required. + + if (sortimid == YES) + call ph_asort (imtable, Memi[imid], Memi[id], Memr[x], Memr[y], + Memi[nap], Memr[rap], Memr[mag], Memr[merr], naperts, npts) + + # Free the image id sorting index arrays. + + if (imid != NULL) + call mfree (imid, TY_INT) + + # Compute the curves of growth. + + if (interactive == YES) { + gd = gopen (Memc[graphics], NEW_FILE, STDGRAPH) + call ph_aigrow (gd, apcat, magfd, logfd, mgd, imtable, Memi[id], + Memr[x], Memr[y], Memi[nap], Memr[rap], Memr[mag], Memr[merr], + naperts, Memd[params], mterms, smallap, largeap) + call gclose (gd) + } else + call ph_agrow (apcat, magfd, logfd, mgd, imtable, Memi[id], Memr[x], + Memr[y], Memi[nap], Memr[rap], Memr[mag], Memr[merr], naperts, + Memd[params], mterms, smallap, largeap) + + # Free the data arrays. + + if (id != NULL) + call mfree (id, TY_INT) + if (x != NULL) + call mfree (x, TY_REAL) + if (y != NULL) + call mfree (y, TY_REAL) + if (nap != NULL) + call mfree (nap, TY_INT) + if (rap != NULL) + call mfree (rap, TY_REAL) + if (mag != NULL) + call mfree (mag, TY_REAL) + if (merr != NULL) + call mfree (merr, TY_REAL) + + # Close the symbol table. + call stclose (imtable) + + # Close the files and file lists. + call close (apcat) + if (magfd != NULL) + call close (magfd) + if (logfd != NULL) + call close (logfd) + if (mgd != NULL) + call gclose (mgd) + if (plotfd != NULL) + call close (plotfd) + if (obs != NULL) + call close (obs) + call clpcls (photfiles) + + call sfree (sp) +end + + +# PH_ASORT -- Sort the data. + +procedure ph_asort (imtable, imid, id, x, y, nap, rap, mag, merr, naperts, npts) + +pointer imtable # pointer to the symbol table +int imid[ARB] # array of image ids +int id[ARB] # array of object ids +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +int nap[ARB] # array of aperture numbers +real rap[naperts,ARB] # the aperture radii list +real mag[naperts,ARB] # array of magnitudes +real merr[naperts,ARB] # array of magnitude errors +int naperts # the number of apertures +int npts # the number of points + +int nimages, i, nptr +pointer sp, sym, symbol +int stnsymbols() +pointer sthead(), stnext() + +begin + nimages = stnsymbols (imtable, 0) + if (nimages <= 0) + return + + # Allocate working space. + call smark (sp) + call salloc (sym, nimages, TY_INT) + + # Read in the symbol list in the last in first out sense. + + symbol = sthead (imtable) + do i = 1, nimages { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Do the image/object id sort if necessary. This is a real sort + # where the actual contents of the arrays are rearcolumnd. + # After the image/object id sort recompute the offsets specifying the + # location of the image data in the symbol table. The first + # pass through the symbol table to picks up the information for the + # images that have only one entry. The second pass picks up + # information for images that have more than one entry. + + call ph_5r3isort (imid, id, nap, x, y, rap, mag, merr, naperts, npts) + nptr = npts + 1 + do i = 1, nimages { + symbol = Memi[sym+i-1] + if (IMT_NENTRIES(symbol) <= 0) + next + nptr = nptr - IMT_NENTRIES(symbol) + IMT_OFFSET(symbol) = nptr + } + + call sfree (sp) +end + + +# PH_ANXTIMAGE -- Find the first line in the input catalog containing the next +# image name, given the current image name. Return the new image nam and line. + +int procedure ph_anxtimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # current maximum line buffer size + +int i, op, colno +long stat +pointer sp, str +bool streq() +int getline(), nscan() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + Memc[str] = EOS + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() != colno) + next + + } until (! streq (image, Memc[str])) + call strcpy (Memc[str], image, SZ_FNAME) + + call sfree (sp) + + return (stat) +end + + +# PH_AGETIMAGE -- Read the next line in the input catalog and return the image +# name and the line. + +int procedure ph_agetimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # size of the output line + +int op, i, colno +int stat +int getline(), nscan() + +begin + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (image, SZ_FNAME) + if (nscan() != colno) + next + + break + } + + return (stat) +end + + +# PH_AIMDATA -- Decode the airmass from a a line of the input catalog. + +procedure ph_aimdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the filter id +int sz_filterid # the maximum size of the filter id +real itime # the integration time +real airmass # the airmass of the observation +real otime # the time of observation + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[CAT_IFILTER] + icol = columns[CAT_ITIME] + acol = columns[CAT_XAIRMASS] + ocol = columns[CAT_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + airmass = INDEFR + + call sscan (line) + do i = 1, maxcol { + if (i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_AOBSDATA -- Decode the filterid, exposure time, and airmass from a +# line of the airmass file. + +procedure ph_aobsdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the filter id +int sz_filterid # maximum size if filter id string +real itime # the exposure time +real airmass # the airmass of the observation +real otime # the time of observation + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[OBS_IFILTER] + icol = columns[OBS_ITIME] + acol = columns[OBS_XAIRMASS] + ocol = columns[OBS_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + call strcpy ("INDEF", filterid, sz_filterid) + itime = INDEFR + airmass = INDEFR + + call sscan (line) + do i = 1, maxcol { + if (i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_ASTARDATA -- Decode the image name, x and y coordinates, magnitude +# and magnitude error from the input catalog. + +int procedure ph_astardata (line, columns, x, y, nap, rap, mag, merr, naperts, + magerr) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +real x, y # the output x and y coordinates +int nap # the number of apertures +real rap[ARB] # the output aperture radii +real mag[ARB] # the output magnitude array +real merr[ARB] # the output magnitude error array +int naperts # the number of apertures +real magerr # the maximum magnitude error + +int i, xcol, ycol, rcol1, rcol2, mcol1, mcol2, ecol1, ecol2, maxcol, stat +pointer sp, str +int nscan() +real rval + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + xcol= columns[CAT_XCENTER] + ycol= columns[CAT_YCENTER] + rcol1 = columns[CAT_RAPERT] + rcol2 = columns[CAT_RAPERT] + naperts - 1 + mcol1 = columns[CAT_MAG] + mcol2 = columns[CAT_MAG] + naperts - 1 + ecol1 = columns[CAT_MERR] + ecol2 = columns[CAT_MERR] + naperts - 1 + maxcol = max (xcol, ycol, rcol1, rcol2, mcol1, mcol2, ecol1, ecol2) + + x = INDEFR + y = INDEFR + nap = 0 + call amovkr (INDEFR, rap, naperts) + call amovkr (INDEFR, mag, naperts) + call amovkr (INDEFR, merr, naperts) + stat = OK + + call sscan (line) + do i = 1, maxcol { + + if (i == xcol) { + call gargr (rval) + if (nscan() != xcol) + stat = ERR + else + x = rval + } else if (i == ycol) { + call gargr (rval) + if (nscan() != ycol) + stat = ERR + else + y = rval + } else if (i >= rcol1 && i <= rcol2) { + call gargr (rval) + if (nscan() != i) + stat = ERR + else + rap[i-rcol1+1] = rval + } else if (i >= mcol1 && i <= mcol2) { + call gargr (rval) + if (nscan() != i) + stat = ERR + #else if (IS_INDEFR(rval)) + #stat = ERR + else { + mag[i-mcol1+1] = rval + nap = nap + 1 + } + } else if (i >= ecol1 && i <= ecol2) { + call gargr (rval) + if (nscan() != i) + stat = ERR + #else if (IS_INDEFR(rval)) + #stat = ERR + #else if (rval > magerr) + #stat = ERR + else + #merr[i-ecol1+1] = sqrt (1.0e-6 + rval ** 2) + merr[i-ecol1+1] = rval + } else { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() != i) + stat = ERR + } + + if (stat == ERR) + break + nap = nap + 1 + } + + call sfree (sp) + + # Compute the magnitude differences. + if (stat == ERR) { + return (ERR) + } else if (nap > 1) { + nap = 0 + do i = 1, naperts { + if (IS_INDEFR(mag[i]) || IS_INDEFR(merr[i]) || merr[i] > + magerr) + break + nap = nap + 1 + } + #if (nap <= 1) + if (nap <= 0) + return (ERR) + if (nap > 1) { + do i = nap, 2, -1 { + mag[i] = mag[i] - mag[i-1] + merr[i] = sqrt (1.0e-6 + merr[i] ** 2) + } + } + return (OK) + } else + return (ERR) +end + + +# PH_AGRANGE -- Get the next column from a list of columns. + +int procedure ph_agrange (list, ip, label, maxch) + +char list[ARB] # list of labels +int ip # pointer in to the list of labels +char label[ARB] # the output label +int maxch # maximum length of a column name + +int op, token +int ctotok(), strlen() + +begin + # Decode the column labels. + op = 1 + while (list[ip] != EOS) { + + token = ctotok (list, ip, label[op], maxch) + if (label[op] == EOS) + next + if ((token == TOK_UNKNOWN) || (token == TOK_CHARCON)) + break + if ((token == TOK_PUNCTUATION) && (label[op] == ',')) { + if (op == 1) + next + else + break + } + + op = op + strlen (label[op]) + if (IS_WHITE(list[ip])) + break + } + + label[op] = EOS + if ((list[ip] == EOS) && (op == 1)) + return (EOF) + else + return (op - 1) +end diff --git a/noao/digiphot/photcal/mkobsfile/t_mkphotcors.x b/noao/digiphot/photcal/mkobsfile/t_mkphotcors.x new file mode 100644 index 00000000..c3e22646 --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/t_mkphotcors.x @@ -0,0 +1,273 @@ +include "../lib/obsfile.h" + +# T_MKPHOTCORS - Query the user for the image sets, observing parameters, +# positional shifts and aperture corrections required by the preprocessor +# task MKOBSFILE. Input to OBSQUERY is the names of the output image set, +# observing parameters, shifts and aperture corrections files, and a list +# of filter ids. + +procedure t_mkphotcors () + +pointer infilters # pointer to the input list of filter ids +pointer imsets # name of the output image set file +pointer obsparams # pointer to the output obsparams file +pointer shifts # pointer to the output shifts file +pointer apercors # pointer to the output aperture corrections file +pointer obscolumnstr # pointer to format of obsparams +int verify # verify interactive user input +int verbose # print status, warning and error messages + +int nfilters, csp, cp, ncolumns, nimages +int insets, outsets, inobs, outobs, insh, outsh, inap, outap +pointer sp, outfilters, imtable, obscolname, obscolumns +bool clgetb() +int btoi(), open(), ph_filters(), ph_esetimtable(), ph_setimtable() +int ctoi(), access(), ph_getlabels() +pointer stopen() + +begin + # Allocate working memory. + call smark (sp) + call salloc (imsets, SZ_FNAME, TY_CHAR) + call salloc (infilters, SZ_LINE, TY_CHAR) + call salloc (outfilters, SZ_LINE, TY_CHAR) + call salloc (obsparams, SZ_FNAME, TY_CHAR) + call salloc (shifts, SZ_FNAME, TY_CHAR) + call salloc (apercors, SZ_FNAME, TY_CHAR) + call salloc (obscolumnstr, SZ_LINE, TY_CHAR) + call salloc (obscolname, SZ_FNAME, TY_CHAR) + call salloc (obscolumns, OBS_NFIELDS, TY_INT) + + # Get the parameters. + call clgstr ("imsets", Memc[imsets], SZ_FNAME) + call clgstr ("idfilters", Memc[infilters], SZ_LINE) + call clgstr ("obsparams", Memc[obsparams], SZ_FNAME) + call clgstr ("shifts", Memc[shifts], SZ_FNAME) + call clgstr ("apercors", Memc[apercors], SZ_FNAME) + call clgstr ("obscolumns", Memc[obscolumnstr], SZ_LINE) + verify = btoi (clgetb ("verify")) + verbose = btoi (clgetb ("verbose")) + + # Open the idfilters string. + nfilters = ph_filters (Memc[infilters], Memc[outfilters], SZ_LINE) + + # Open the input and output image sets file. + if (Memc[imsets] == EOS) { + outsets = NULL + insets = open ("STDIN", READ_ONLY, TEXT_FILE) + } else if (access (Memc[imsets], READ_ONLY, TEXT_FILE) == YES) { + outsets = NULL + insets = open (Memc[imsets], READ_ONLY, TEXT_FILE) + } else { + outsets = open (Memc[imsets], NEW_FILE, TEXT_FILE) + insets = open ("STDIN", READ_ONLY, TEXT_FILE) + } + + # Open the input and output observing parameters file. + if (Memc[obsparams] == EOS) { + outobs = NULL + inobs = NULL + } else if (access (Memc[obsparams], READ_ONLY, TEXT_FILE) == YES) { + outobs = NULL + inobs = open (Memc[obsparams], READ_ONLY, TEXT_FILE) + call amovki (0, Memi[obscolumns], OBS_NFIELDS) + csp = 1 + ncolumns = 0 + while (ph_getlabels (Memc[obscolumnstr], csp, Memc[obscolname], + SZ_FNAME) != EOF) { + cp = 1 + if (ctoi (Memc[obscolname], cp, Memi[obscolumns+ncolumns]) <= 0) + Memi[obscolumns+ncolumns] = 0 + ncolumns = ncolumns + 1 + } + } else { + outobs = open (Memc[obsparams], NEW_FILE, TEXT_FILE) + inobs = open ("STDIN", READ_ONLY, TEXT_FILE) + } + + # Open the input and output shifts file. + if (Memc[shifts] == EOS) { + outsh = NULL + insh = NULL + } else if (access (Memc[shifts], READ_ONLY, TEXT_FILE) == YES) { + outsh = NULL + insh = open (Memc[shifts], READ_ONLY, TEXT_FILE) + } else { + outsh = open (Memc[shifts], NEW_FILE, TEXT_FILE) + insh = open ("STDIN", READ_ONLY, TEXT_FILE) + } + + # Open the aperture corrections file. + if (Memc[apercors] == EOS) { + outap = NULL + inap = NULL + } else if (access (Memc[apercors], READ_ONLY, TEXT_FILE) == YES) { + outap = NULL + inap = open (Memc[apercors], READ_ONLY, TEXT_FILE) + } else { + outap = open (Memc[apercors], NEW_FILE, TEXT_FILE) + inap = open ("STDIN", READ_ONLY, TEXT_FILE) + } + + # Open a symbol table to hold a list of image charactersitics. + + imtable = stopen ("imtable", 2 * LEN_IMTABLE, LEN_IMTABLE, 10 * + LEN_IMTABLE) + + # Read in the image set file and store each unique image name in + # the symbol table. Initialize the records fields. + + if (insets == STDIN) + nimages = ph_esetimtable (imtable, nfilters, verify) + else + nimages = ph_setimtable (imtable, insets, nfilters, verbose) + + # For each image referenced in the observing parameters file extract + # the image name, the filter id if defined, the exposure time if + # defined and the airmass if defined. Enter the new values of filter + # id and airmass in the symbol table and compute the effective + # exposure time required to correct the instrumental magnitudes to + # the new exposure time. + + if (inobs == STDIN) + call ph_eobsimtable (imtable, nimages, verify) + else if (inobs != NULL) + call ph_obsimtable (imtable, inobs, Memi[obscolumns], verbose) + + # Read in the x and y shifts for the image. Images for which no + # shift is defined are assigned an x-y shift of 0.0. + + if (insh == STDIN) + call ph_eshimtable (imtable, nimages, verify) + else if (insh != NULL) + call ph_shimtable (imtable, insh, verbose) + + # Read in the aperture corrections. Images for which no aperture + # correction is defined are assigned an aperture correction of 0.0. + + if (inap == STDIN) + call ph_eapimtable (imtable, nimages, verify) + else if (inap != NULL) + call ph_apimtable (imtable, inap, verbose) + + # Write the results to the various output files. + if (nimages > 0) + call ph_wimtable (imtable, outsets, outobs, outsh, outap, nimages) + + # Close the symbol table. + call stclose (imtable) + + # Close the input files. + if (insets != NULL) + call close (insets) + if (inobs != NULL) + call close (inobs) + if (insh != NULL) + call close (insh) + if (inap != NULL) + call close (inap) + + # Close the output files. + if (outsets != NULL) { + call close (outsets) + if (nimages <= 0) + call delete (Memc[imsets]) + } + if (outobs != NULL) { + call close (outobs) + if (nimages <= 0) + call delete (Memc[obsparams]) + } + if (outsh != NULL) { + call close (outsh) + if (nimages <= 0) + call delete (Memc[shifts]) + } + if (outap != NULL) { + call close (outap) + if (nimages <= 0) + call delete (Memc[apercors]) + } + + call sfree (sp) +end + + +# PH_WIMTABLE -- Write the contents of the symbol table to the output image +# set file, the output observing parameters file, the output shifts file, +# and the output aperture corrections file. The only assumption here is +# that all images in a given image set are adjacent to each other in the +# symbol table, an assumption that is always true. + +procedure ph_wimtable (imtable, outsets, outobs, outsh, outap, nimages) + +pointer imtable # pointer to the the input symbol table +int outsets # the output image set file descriptor +int outobs # the output observing parameters file descriptor +int outsh # the output positional shifts file descriptor +int outap # the output aperture corrections file descriptor +int nimages # number of images in the symbol table + +int i, osetno, setno +pointer sp, sym, symbol +pointer sthead(), stnext() + +begin + call smark (sp) + call salloc (sym, nimages, TY_POINTER) + + # Reorder the list of symbols. + symbol = sthead (imtable) + do i = nimages, 1, -1 { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Write the output files. + osetno = 0 + do i = 1, nimages { + + symbol = Memi[sym+i-1] + setno = IMT_IMSETNO(symbol) + + if (outsets != NULL) { + if (setno != osetno) { + if (setno == 1) + call fprintf (outsets, "%s :") + else + call fprintf (outsets, "\n%s :") + call pargstr (IMT_LABEL(symbol)) + } + call fprintf (outsets, " %s") + call pargstr (IMT_IMNAME(symbol)) + if (i == nimages) + call fprintf (outsets, "\n") + } + + if (outobs != NULL) { + call fprintf (outobs, "%s %s %g %g %0.1h\n") + call pargstr (IMT_IMNAME(symbol)) + call pargstr (IMT_IFILTER(symbol)) + call pargr (IMT_ITIME(symbol)) + call pargr (IMT_XAIRMASS(symbol)) + call pargr (IMT_OTIME(symbol)) + } + + if (outsh != NULL) { + call fprintf (outsh, "%s %g %g\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (IMT_XSHIFT(symbol)) + call pargr (IMT_YSHIFT(symbol)) + } + + if (outap != NULL) { + call fprintf (outap, "%s %g\n") + call pargstr (IMT_IMNAME(symbol)) + call pargr (IMT_APERCOR(symbol)) + } + + osetno = setno + } + + call sfree (sp) +end diff --git a/noao/digiphot/photcal/mkobsfile/t_obsfile.x b/noao/digiphot/photcal/mkobsfile/t_obsfile.x new file mode 100644 index 00000000..6c997cea --- /dev/null +++ b/noao/digiphot/photcal/mkobsfile/t_obsfile.x @@ -0,0 +1,1131 @@ +include <fset.h> +include <ctotok.h> +include <ctype.h> +include "../lib/obsfile.h" + +# T_OBSFILE -- Make an observations catalog using the data extracted +# from the APPHOT/DAOPHOT catalogs by TXDUMP/TDUMP or supplied in a simple +# text file by the user, the observing parameters file, the shifts and +# aperture corrections file and the image sets file. OBSFILE expects to +# see the following 8 fields in the input data file: image name, xcenter, +# ycenter, exposure time, filter id, airmass, magnitude and magnitude error +# and uses a simple string containing a list of column numbers to determine +# which data is in which column. The observing parameters file is decoded +# in the same way. + +procedure t_obsfile () + +int photfiles # pointer to the list of photometry files +pointer columnstr # pointer to the list of input columns +pointer infilters # pointer to the input list of filter ids +int normtime # normalize to an exposure time of 1 time unit +int allfilters # only output objects with data in all filters +real tolerance # tolerance in pixels for position matching +int verify # verify interactive user input +int verbose # print status, warning and errors messages + +int incat, sets, sh, ap, obs, outcat, fmt, wrap +int csp, cp, nincolumns, nfilters, nimages, npts, sortimid +pointer sp, outfilters, fname, tfname, column, incolumns, obscolumns, imtable +pointer objid, x, y, mag, merr, imid, id +real minmagerr + +bool clgetb() +int clpopnu(), clgfil(), open(), clplen(), ph_filters(), ph_getlabels() +int ctoi(), ph_setimtable(), ph_esetimtable(), ph_mkimtable(), btoi() +int access() +pointer stopen() +real clgetr() + +begin + # Allocate working memory. + call smark (sp) + + call salloc (infilters, SZ_LINE, TY_CHAR) + call salloc (outfilters, SZ_LINE, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (tfname, SZ_FNAME, TY_CHAR) + + call salloc (columnstr, SZ_LINE, TY_CHAR) + call salloc (column, SZ_FNAME, TY_CHAR) + call salloc (incolumns, CAT_NFIELDS, TY_INT) + call salloc (obscolumns, OBS_NFIELDS, TY_INT) + + # Get the input file list. + photfiles = clpopnu ("photfiles") + if (clplen (photfiles) <= 0) + call error (0, "The input file list is empty") + + # Decode the input columns list. There must be at least CAT_NFIELDS + # column numbers in the incolumns list. The number 0 can be used as + # a place holder for missing columns. + + call clgstr ("incolumns", Memc[columnstr], SZ_LINE) + csp = 1 + nincolumns = 0 + while (ph_getlabels (Memc[columnstr], csp, Memc[column], SZ_FNAME) != + EOF) { + cp = 1 + if (ctoi (Memc[column], cp, Memi[incolumns+nincolumns]) <= 0) + Memi[incolumns+nincolumns] = 0 + nincolumns = nincolumns + 1 + } + if (nincolumns < CAT_NFIELDS) + call error (0, + "There are too few columns defined in the incolumns string") + + # Open the idfilters string. + call clgstr ("idfilters", Memc[infilters], SZ_LINE) + nfilters = ph_filters (Memc[infilters], Memc[outfilters], SZ_LINE) + + # Open the image sets file. + call clgstr ("imsets", Memc[fname], SZ_FNAME) + sets = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the output catalog file. + call clgstr ("observations", Memc[fname], SZ_FNAME) + outcat = open (Memc[fname], NEW_FILE, TEXT_FILE) + call sprintf (Memc[tfname], SZ_FNAME, "f%s.dat") + call pargstr (Memc[fname]) + if (access (Memc[tfname], 0, 0) == YES) + call delete (Memc[tfname]) + fmt = open (Memc[tfname], NEW_FILE, TEXT_FILE) + wrap = btoi (clgetb ("wrap")) + + # Get the minimum magnitude error. + minmagerr = clgetr ("minmagerr") + + # Normalize the exposure times? + normtime = btoi (clgetb ("normtime")) + + # Get the position matching parameters. + tolerance = clgetr ("tolerance") + allfilters = btoi (clgetb ("allfilters")) + + # Verify interactive user input. + verify = btoi (clgetb ("verify")) + verbose = btoi (clgetb ("verbose")) + + # Open the shifts file. + call clgstr ("shifts", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + sh = NULL + else + sh = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the aperture corrections file. + call clgstr ("apercors", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + ap = NULL + else + ap = open (Memc[fname], READ_ONLY, TEXT_FILE) + + # Open the observing parameters file and decode the column list. + + call clgstr ("obsparams", Memc[fname], SZ_FNAME) + if (Memc[fname] == EOS) + obs = NULL + else { + obs = open (Memc[fname], READ_ONLY, TEXT_FILE) + call amovki (0, Memi[obscolumns], OBS_NFIELDS) + if (obs != STDIN) { + call clgstr ("obscolumns", Memc[columnstr], SZ_LINE) + csp = 1 + nincolumns = 1 + Memi[obscolumns] = 1 + while (ph_getlabels (Memc[columnstr], csp, Memc[column], + SZ_FNAME) != EOF && (nincolumns < OBS_NFIELDS)) { + cp = 1 + if (ctoi (Memc[column], cp, + Memi[obscolumns+nincolumns]) <= 0) + Memi[obscolumns+nincolumns] = 0 + else if (Memi[obscolumns+nincolumns] == 1) + nincolumns = nincolumns - 1 + nincolumns = nincolumns + 1 + } + } + } + + # Open a symbol table to hold a list of image charactersitics. + + imtable = stopen ("imtable", 2 * LEN_IMTABLE, LEN_IMTABLE, 10 * + LEN_IMTABLE) + + # Read in the image set file and store each unique image name in + # the symbol table. Initialize the records fields. + + if (sets == NULL) + nimages = 0 + else if (sets == STDIN) + nimages = ph_esetimtable (imtable, nfilters, verify) + else + nimages = ph_setimtable (imtable, sets, nfilters, verbose) + + # For each image referenced in the observing parameters file extract + # the image name, the filter id if defined, the exposure time if + # defined and the airmass if defined. Enter the new values of filter + # id and airmass in the symbol table and compute the effective + # exposure time required to correct the instrumental magnitudes to + # the new exposure time. + + if (obs != NULL) { + if (obs == STDIN) + call ph_eobsimtable (imtable, nimages, verify) + else { + if (sets == STDIN && verbose == YES) { + call fstats (obs, F_FILENAME, Memc[fname], SZ_FNAME) + call eprintf ("Warning: Filter ids in %s ") + call pargstr (Memc[fname]) + call eprintf ("will replace those listed in idfilters\n") + } + call ph_obsimtable (imtable, obs, Memi[obscolumns], verbose) + } + } + + # Read in the x and y shifts for the image. Images for which no + # shift is defined are assigned an x-y shift of 0.0. + + if (sh != NULL) { + if (sh == STDIN) + call ph_eshimtable (imtable, nimages, verify) + else + call ph_shimtable (imtable, sh, verbose) + } + + # Read in the aperture corrections. Images for which no aperture + # correction is defined are assigned an aperture correction of 0.0. + + if (ap != NULL) { + if (ap == STDIN) + call ph_eapimtable (imtable, nimages, verify) + else + call ph_apimtable (imtable, ap, verbose) + } + + # Read in the data. Extract each unique image name and match it + # with the corresponding image name in the symbol table. Skip data + # for images which are not in the symbol table. For each unique image + # name which does match a symbol table image name, extract the exposure + # time, filter id, and airmass and enter them in the symbol table. + # Read the x and y and magnitude and magnitude error data into + # separate data arrays. + + objid = NULL + x = NULL + y = NULL + mag = NULL + merr = NULL + imid = NULL + id = NULL + sortimid = NO + + npts = 0 + while (clgfil (photfiles, Memc[fname], SZ_FNAME) != EOF) { + incat = open (Memc[fname], READ_ONLY, TEXT_FILE) + npts = ph_mkimtable (imtable, incat, Memi[incolumns], objid, x, + y, mag, merr, imid, id, npts, normtime, sortimid, verbose) + call close (incat) + } + + if (objid != NULL) + call realloc (objid, npts * (DEF_LENLABEL+1), TY_CHAR) + call realloc (x, npts, TY_REAL) + call realloc (y, npts, TY_REAL) + call realloc (mag, npts, TY_REAL) + call realloc (merr, npts, TY_REAL) + call realloc (imid, npts, TY_INT) + call realloc (id, npts, TY_INT) + + # Sort the data by image id and object id if position matching is + # turned off or by image id and y coordinate if position matching is + # turned on. + + if (objid == NULL) + call ph_sort (imtable, nimages, "", Memr[x], Memr[y], Memr[mag], + Memr[merr], Memi[imid], Memi[id], npts, sortimid, tolerance) + else + call ph_sort (imtable, nimages, Memc[objid], Memr[x], Memr[y], + Memr[mag], Memr[merr], Memi[imid], Memi[id], npts, sortimid, + tolerance) + # Free the image id sorting index arrays. + + if (imid != NULL) + call mfree (imid, TY_INT) + + # Apply the exposure time corrections, x-y shifts, and the aperture + # corrections. + + call ph_correct (imtable, Memr[x], Memr[y], Memr[mag], Memr[merr], + minmagerr) + + if (verbose == YES) { + call fstats (outcat, F_FILENAME, Memc[fname], SZ_FNAME) + call printf ("\nObservations file: %s\n") + call pargstr (Memc[fname]) + } + + # Match by order in the catalog if tolerance is less than or equal to + # zero or by x-y position if tolerance is greater than zero and output + # the data. + + if (objid == NULL) + call ph_match (imtable, nimages, Memc[outfilters], nfilters, outcat, + "", Memr[x], Memr[y], Memr[mag], Memr[merr], Memi[id], npts, + tolerance, allfilters, wrap, verbose) + else + call ph_match (imtable, nimages, Memc[outfilters], nfilters, outcat, + Memc[objid], Memr[x], Memr[y], Memr[mag], Memr[merr], + Memi[id], npts, tolerance, allfilters, wrap, verbose) + + if (verbose == YES) + call printf ("\n") + + # Write the format file for the output catalog. + + call ph_format (fmt, Memc[outfilters], nfilters) + + # Free the data arrays. + + if (objid != NULL) + call mfree (objid, TY_CHAR) + if (x != NULL) + call mfree (x, TY_REAL) + if (y != NULL) + call mfree (y, TY_REAL) + if (mag != NULL) + call mfree (mag, TY_REAL) + if (merr != NULL) + call mfree (merr, TY_REAL) + if (id != NULL) + call mfree (id, TY_INT) + + # Close the symbol table. + call stclose (imtable) + + # Close the files and file lists. + if (sets != NULL) + call close (sets) + if (sh != NULL) + call close (sh) + if (ap != NULL) + call close (ap) + if (obs != NULL) + call close (obs) + call close (outcat) + call close (fmt) + call clpcls (photfiles) + + call sfree (sp) +end + + +# PH_SORT -- Sort the data. + +procedure ph_sort (imtable, nimages, objid, x, y, mag, merr, imid, id, npts, + sortimid, tolerance) + +pointer imtable # pointer to the symbol table +int nimages # the number of images in the symbol table +char objid[ARB] # the array of object ids objid[1] = EOS if none +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real merr[ARB] # array of magnitude errors +int imid[ARB] # array of image ids +int id[ARB] # array of object ids +int npts # the number of points +int sortimid # sort by image and id +real tolerance # the tolerance in pixels + +int i, nptr +pointer sp, sym, symbol +pointer sthead(), stnext() + +begin + # Allocate working space. + call smark (sp) + call salloc (sym, nimages, TY_INT) + + # Read in the symbol list in the last in first out sense. + + symbol = sthead (imtable) + do i = 1, nimages { + Memi[sym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Do the image/object id sort if necessary. This is a real sort + # where the actual contents of the arrays are rearranged. + # After the image/object id sort recompute the offsets specifying the + # location of the image data in the symbol table. The first + # pass through the symbol table to picks up the information for the + # images that have only one entry. The second pass picks up + # information for images that have more than one entry. + + if (sortimid == YES) { + if (objid[1] == EOS) + call ph_4r2isort (x, y, mag, merr, imid, id, npts) + else + call ph_1c4r2isort (objid, x, y, mag, merr, imid, id, npts) + nptr = npts + 1 + do i = 1, nimages { + symbol = Memi[sym+i-1] + if (IMT_NENTRIES(symbol) <= 0) + next + nptr = nptr - IMT_NENTRIES(symbol) + IMT_OFFSET(symbol) = nptr + } + } + + # Sort on y if necessary. This is an index sort and does not alter + # the position of the data. + + if (tolerance > 0) { + do i = 1, nimages { + symbol = Memi[sym+i-1] + if (IMT_NENTRIES(symbol) <= 0) + next + nptr = IMT_OFFSET(symbol) + call ph_qsort (y[nptr], id[nptr], IMT_NENTRIES(symbol), nptr) + } + } + + call sfree (sp) +end + + +# PH_CORRECT -- Correct the x, y and mag values using the x-y shifts, aperture +# corrections, and exposure times. + +procedure ph_correct (imtable, x, y, mag, merr, minmagerr) + +pointer imtable # pointer to the symbol table +real x[ARB] # array of x coordinates +real y[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real merr[ARB] # array of magnitude errors +real minmagerr # the minimum magnitude error + +int i, nptr, ndata +pointer sp, imname, symbol, osymbol +real magshift +pointer sthead(), stnext(), stfind() + +begin + # Allocate working space. + + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + # Loop through the images adding the x, y and magnitude shifts to + # the data. + + symbol = sthead (imtable) + while (symbol != NULL) { + + # Get the correct symbol. + ndata = IMT_NENTRIES(symbol) + osymbol = stfind (imtable, IMT_IMNAME(symbol)) + + # If data is present make the corrections once for each image. + + if ((ndata > 0) && (osymbol == symbol)) { + nptr = IMT_OFFSET(symbol) + + # Correct the positions. + call aaddkr (x[nptr], IMT_XSHIFT(symbol), x[nptr], ndata) + call aaddkr (y[nptr], IMT_YSHIFT(symbol), y[nptr], ndata) + + # Correct the magnitudes. + magshift = 2.5 * log10 (IMT_ITIME(symbol)) + IMT_APERCOR(symbol) + do i = nptr, nptr + ndata - 1 { + if (IS_INDEFR(mag[i])) + next + mag[i] = mag[i] + magshift + if (IS_INDEFR(merr[i])) + next + if (merr[i] < minmagerr) + merr[i] = minmagerr + } + } + + # Get the next symbol. + symbol = stnext (imtable, symbol) + } + + call sfree (sp) +end + + +# PH_MATCH -- Match up the photometry lists by order in the input catalog +# or x-y position. + +procedure ph_match (imtable, nimages, filters, nfilters, out, objid, x, y, + mag, merr, ysort, npts, tolerance, allfilters, wrap, verbose) + +pointer imtable # pointer to the image name symbol table +int nimages # the number of images in the symbol table +char filters[ARB] # the filters string +int nfilters # the number of filters +int out # the output file descriptor +char objid[ARB] # the object id array +real x[ARB] # the x input array +real y[ARB] # the y input array +real mag[ARB] # the magnitude array +real merr[ARB] # the magnitude error array +int ysort[ARB] # the y sort index +int npts # the number of points +real tolerance # the x-y matching tolerance in pixels +int allfilters # match in all filters +int wrap # format the output file for easy reading +int verbose # print status, warning and error messages + +int i, k, nsets, imptr, filterno, ndata +pointer sp, imsym, filterid, label, sym, match, osymbol, symbol, index +bool streq() +int strdic(), ph_nthlabel() +pointer sthead(), stnext(), stfind() + +begin + # Write out the output file column headers. + call fprintf (out, "\n") + if (wrap == YES) { + call fprintf (out, +"# FIELD%17tFILTER%34tOTIME%40tAIRMASS%49tXCENTER%59tYCENTER%71tMAG%77tMERR\n") + } else { + do i = 1, nfilters { + if (i == 1) + call fprintf (out, +"# FIELD%17tFILTER%34tOTIME%40tAIRMASS%49tXCENTER%59tYCENTER%71tMAG%77tMERR ") + else + call fprintf (out, +"%2tFILTER%19tOTIME%25tAIRMASS%34tXCENTER%44tYCENTER%56tMAG%62tMERR ") + } + call fprintf (out,"\n") + } + call fprintf (out, "\n") + + # Allocate and or initialize working space. + call smark (sp) + call salloc (imsym, nimages, TY_INT) + call salloc (filterid, SZ_FNAME, TY_CHAR) + call salloc (label, SZ_FNAME, TY_CHAR) + call salloc (sym, nfilters, TY_INT) + if (tolerance > 0.0) + call malloc (match, npts, TY_INT) + else + match = NULL + index = NULL + + + # Read in the symbol table pointers and rearrange them so that the + # image name symbols are ordered in the same way as the image set + # numbers. Compute the total number of image sets at the same time. + + nsets = 0 + symbol = sthead (imtable) + do i = nimages, 1, -1 { + nsets = max (nsets, IMT_IMSETNO(symbol)) + Memi[imsym+i-1] = symbol + symbol = stnext (imtable, symbol) + } + + # Loop over the image sets. + imptr = 1 + do i = 1, nsets { + + # Clear the array of symbols and set the label. + ndata = 0 + call amovki (NULL, Memi[sym], nfilters) + call strcpy (IMT_LABEL(Memi[imsym+imptr-1]), Memc[label], SZ_FNAME) + + # Loop over the images in the set. + repeat { + + # Quit if the image pointer exceeds the number of images. + if (imptr > nimages) + break + + # Get the next symbol for a given image set. + symbol = Memi[imsym+imptr-1] + if (IMT_IMSETNO(symbol) != i) + break + + # Check for duplicate image names. + osymbol = stfind (imtable, IMT_IMNAME(symbol)) + if (osymbol != symbol) { + IMT_OFFSET(symbol) = IMT_OFFSET(osymbol) + IMT_NENTRIES(symbol) = IMT_NENTRIES(osymbol) + IMT_XAIRMASS(symbol) = IMT_XAIRMASS(osymbol) + IMT_OTIME(symbol) = IMT_OTIME(osymbol) + call strcpy (IMT_IFILTER(osymbol), IMT_IFILTER(symbol), + SZ_FNAME) + } + + # Search for the filter id, if one is not found then + # use the first unoccupied space. + + filterno = strdic (IMT_IFILTER(symbol), Memc[filterid], + SZ_FNAME, filters) + if (streq (IMT_IFILTER(symbol), "INDEF")) { + do k = 1, nfilters { + if (Memi[sym+k-1] != NULL) + next + if (ph_nthlabel (filters, k, Memc[filterid], + SZ_FNAME) != k) + Memc[filterid] = EOS + filterno = k + break + } + } + + # Load the data pointers. + if (filterno <= 0) { + + if (verbose == YES) { + call eprintf ( + "Warning: Image set: %d label: %s image: %s\n") + call pargi (i) + call pargstr (IMT_LABEL(symbol)) + call pargstr (IMT_IMNAME(symbol)) + call eprintf ( + "\tWarning: Filter %s ") + call pargstr (IMT_IFILTER(symbol)) + call eprintf ("does not belong to filter set %s\n") + call pargstr (filters) + } + + } else if (Memi[sym+filterno-1] != NULL) { + + if (verbose == YES) { + call eprintf ( + "Warning: Image set: %d label: %s image: %s\n") + call pargi (i) + call pargstr (IMT_LABEL(symbol)) + call pargstr (IMT_IMNAME(symbol)) + call eprintf ("\tData for filter %s is redundant\n") + call pargstr (Memc[filterid]) + } + + } else { + + Memi[sym+filterno-1] = symbol + ndata = max (ndata, IMT_NENTRIES(symbol)) + } + + imptr = imptr + 1 + + } + + # Skip matching if there are no points. + if (ndata <= 0) { + if (verbose == YES) { + call eprintf ("Image set: %d Label: %s has no data\n") + call pargi (i) + call pargstr (Memc[label]) + } + next + } + + # Do the matching. + if (tolerance <= 0.0) { + if (objid[1] == EOS) + call ph_join (out, Memc[label], Memi[sym], filters, + nfilters, "", x, y, mag, merr, allfilters, wrap, + verbose) + else + call ph_join (out, Memc[label], Memi[sym], filters, + nfilters, objid, x, y, mag, merr, allfilters, wrap, + verbose) + } else { + if (index == NULL) + call malloc (index, nfilters * ndata, TY_INT) + else + call realloc (index, nfilters * ndata, TY_INT) + if (objid[1] == EOS) + call ph_merge (out, Memc[label], Memi[sym], filters, + nfilters, "", x, y, mag, merr, ysort, Memi[match], + Memi[index], ndata, tolerance, allfilters, wrap, + verbose) + else + call ph_merge (out, Memc[label], Memi[sym], filters, + nfilters, objid, x, y, mag, merr, ysort, Memi[match], + Memi[index], ndata, tolerance, allfilters, wrap, + verbose) + } + } + + # Free working space. + if (match != NULL) + call mfree (match, TY_INT) + if (index != NULL) + call mfree (index, TY_INT) + call sfree (sp) +end + + +# PH_FORMAT -- Write the format file for the output catalog. + +procedure ph_format (fd, filters, nfilters) + +int fd # file descripter for the format file +char filters[ARB] # list of filter ids +int nfilters # number of filters + +int i, col +pointer sp, filterid +int ph_nthlabel() + +begin + call smark (sp) + call salloc (filterid, SZ_FNAME, TY_CHAR) + + call fprintf (fd, "# Declare the observations file variables\n\n") + call fprintf (fd, "observations\n\n") + + col = FIRST_COLUMN + do i = 1, nfilters { + if (ph_nthlabel (filters, i, Memc[filterid], SZ_FNAME) != i) + Memc[filterid] = EOS + call fprintf (fd, + "T%s%15t%d%30t# time of observation in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "X%s%15t%d%30t# airmass in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "x%s%15t%d%30t# x coordinate in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, "y%s%15t%d%30t# y coordinate in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, + "m%s%15t%d%30t# instrumental magnitude in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + col = col + 1 + call fprintf (fd, + "error(m%s)%15t%d%30t# magnitude error in filter %s\n") + call pargstr (Memc[filterid]) + call pargi (col) + call pargstr (Memc[filterid]) + call fprintf (fd, "\n") + col = col + DELTA_COLUMN + } + + call sfree (sp) +end + + +# PH_NXTIMAGE -- Find the first line in the input catalog containing the next +# image name, given the current image name. Return the new image name and line. + +int procedure ph_nxtimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # the line buffer size + +int i, op, colno +long stat +pointer sp, str +bool streq() +int getline(), nscan() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + Memc[str] = EOS + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() != colno) + next + + } until (! streq (image, Memc[str])) + call strcpy (Memc[str], image, SZ_FNAME) + + call sfree (sp) + + return (stat) +end + + +# PH_GETIMAGE -- Read the next line in the input catalog and return the image +# name and the line. + +int procedure ph_getimage (fd, columns, image, line, lbufsize) + +int fd # file descriptor of the input text file +int columns[ARB] # the list of input columns +char image[ARB] # the name of the current image +pointer line # pointer to the output line +int lbufsize # the size of the output buffer + +int op, i, colno +int stat +int getline(), nscan() + +begin + if (line == NULL || lbufsize <= 0) { + lbufsize = lbufsize + SZ_LINE + call malloc (line, lbufsize, TY_CHAR) + } + + colno = columns[CAT_IMAGE] + repeat { + + op = 1 + repeat { + stat = getline (fd, Memc[line+op-1]) + if (stat == EOF) + break + op = op + stat + if (op > lbufsize) { + lbufsize = lbufsize + SZ_LINE + call realloc (line, lbufsize, TY_CHAR) + } + } until (Memc[line+op-2] == '\n') + if (stat == EOF) + break + else + stat = op - 1 + + call sscan (Memc[line]) + do i = 1, colno + call gargwrd (image, SZ_FNAME) + if (nscan() != colno) + next + + break + } + + return (stat) +end + + +# PH_IMDATA -- Decode the filter id, integration time and airmass from +# a line of the input catalog. + +procedure ph_imdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the name of the filter +int sz_filterid # maximum length of the filter id +real itime # the exposure time of the observation +real airmass # the airmass of the observation +real otime # the time of observations + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[CAT_IFILTER] + icol = columns[CAT_ITIME] + acol = columns[CAT_XAIRMASS] + ocol = columns[CAT_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + call strcpy ("INDEF", filterid, sz_filterid) + itime = INDEFR + airmass = INDEFR + otime = INDEFR + + call sscan (line) + do i = 1, maxcol { + if ( i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_OBSDATA -- Decode the filter id, integration time and airmass from +# a line of the observing parameters file. + +procedure ph_obsdata (line, columns, filterid, sz_filterid, itime, airmass, + otime) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char filterid[ARB] # the name of the filter +int sz_filterid # maximum length of the filter id +real itime # the exposure time of the observation +real airmass # the airmass of the observation +real otime # the time of observation + +int i, fcol, icol, acol, ocol, maxcol +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + fcol = columns[OBS_IFILTER] + icol = columns[OBS_ITIME] + acol = columns[OBS_XAIRMASS] + ocol = columns[OBS_OTIME] + maxcol = max (fcol, icol, acol, ocol) + + call strcpy ("INDEF", filterid, sz_filterid) + itime = INDEFR + airmass = INDEFR + otime = INDEFR + + call sscan (line) + do i = 1, maxcol { + if ( i == fcol) { + call gargwrd (filterid, sz_filterid) + if (nscan() != fcol) + call strcpy ("INDEF", filterid, sz_filterid) + } else if (i == icol) { + call gargr (itime) + if (nscan() != icol) + itime = INDEFR + } else if (i == acol) { + call gargr (airmass) + if (nscan() != acol) + airmass = INDEFR + } else if (i == ocol) { + call gargr (otime) + if (nscan() != ocol) + otime = INDEFR + } else { + call gargwrd (Memc[str], SZ_LINE) + } + } + + call sfree (sp) +end + + +# PH_STARDATA -- Decode the image name, x and y coordinates, magnitude +# and magnitude error from the input catalog. + +int procedure ph_stardata (line, columns, objid, x, y, mag, merr) + +char line[ARB] # input line to be scanned +int columns[ARB] # the list of input columns +char objid[ARB] # the object id +real x, y # the output x and y coordinates +real mag, merr # the output magnitude and magnitude error + +int i, icol, xcol, ycol, mcol, ecol, maxcol, stat +pointer sp, str +int nscan() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + icol = columns[CAT_ID] + xcol= columns[CAT_XCENTER] + ycol= columns[CAT_YCENTER] + mcol = columns[CAT_MAG] + ecol = columns[CAT_MERR] + maxcol = max (icol, xcol, ycol, mcol, ecol) + + objid[1] = EOS + x = INDEFR + y = INDEFR + mag = INDEFR + merr = INDEFR + stat = OK + + call sscan (line) + do i = 1, maxcol { + + if (i == icol) { + call gargwrd (objid, DEF_LENLABEL) + if (nscan() != icol) + stat = ERR + } else if (i == xcol) { + call gargr (x) + if (nscan() != xcol) + stat = ERR + } else if (i == ycol) { + call gargr (y) + if (nscan() != ycol) + stat = ERR + } else if (i == mcol) { + call gargr (mag) + if (nscan() != mcol) + stat = ERR + } else if (i == ecol) { + call gargr (merr) + if (nscan() != ecol) + stat = ERR + } else { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() != i) + stat = ERR + } + } + + call sfree (sp) + + return (stat) +end + + +# PH_FILTERS -- Reformat the filter string so that it is a suitable string +# dictionary for the STRDIC routine. + +int procedure ph_filters (infilters, outfilters, max_lenfilters) + +char infilters[ARB] # input list of filters +char outfilters[ARB] # output list of filters +int max_lenfilters # maximum length of the filter string + +int ip, nfilters +pointer sp, filter +int ph_getlabels() + +begin + call smark (sp) + call salloc (filter, max_lenfilters, TY_CHAR) + + ip = 1 + nfilters = 0 + outfilters[1] = EOS + while (ph_getlabels (infilters, ip, Memc[filter], max_lenfilters) != + EOF) { + call strcat (",", outfilters, max_lenfilters) + call strcat (Memc[filter], outfilters, max_lenfilters) + nfilters = nfilters + 1 + } + + call sfree (sp) + + return (nfilters) +end + + +# PH_NTHLABEL -- Get the nth label out of a list of labels. + +int procedure ph_nthlabel (list, item, label, maxch) + +char list[ARB] # input list +int item # item to be extracted +char label[ARB] # extracted label +int maxch # maximum length of the extracted label + +int ip, nitems +int ph_getlabels() + +begin + nitems = 0 + + ip = 1 + while (ph_getlabels (list, ip, label, maxch) != EOF) { + nitems = nitems + 1 + if (nitems == item) + break + } + + return (nitems) +end + + +# PH_GETLABELS -- Get the next label from a list of labels. + +int procedure ph_getlabels (list, ip, label, maxch) + +char list[ARB] # list of labels +int ip # pointer in to the list of labels +char label[ARB] # the output label +int maxch # maximum length of a column name + +int op, token +int ctotok(), strlen() + +begin + # Decode the column labels. + op = 1 + while (list[ip] != EOS) { + + token = ctotok (list, ip, label[op], maxch) + if (label[op] == EOS) + next + if ((token == TOK_UNKNOWN) || (token == TOK_CHARCON)) + break + if ((token == TOK_PUNCTUATION) && (label[op] == ',')) { + if (op == 1) + next + else + break + } + + op = op + strlen (label[op]) + if (IS_WHITE(list[ip])) + break + } + + label[op] = EOS + if ((list[ip] == EOS) && (op == 1)) + return (EOF) + else + return (op - 1) +end |