From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/astutil/pdm/README | 9 + noao/astutil/pdm/TODO | 4 + noao/astutil/pdm/mkpkg | 40 +++++ noao/astutil/pdm/pdm.h | 77 +++++++++ noao/astutil/pdm/pdmalltheta.x | 104 +++++++++++ noao/astutil/pdm/pdmampep.x | 38 ++++ noao/astutil/pdm/pdmautorang.x | 101 +++++++++++ noao/astutil/pdm/pdmbatch.x | 49 ++++++ noao/astutil/pdm/pdmclose.x | 63 +++++++ noao/astutil/pdm/pdmcolon.x | 292 +++++++++++++++++++++++++++++++ noao/astutil/pdm/pdmcompare.x | 25 +++ noao/astutil/pdm/pdmcursor.x | 383 +++++++++++++++++++++++++++++++++++++++++ noao/astutil/pdm/pdmdelete.x | 103 +++++++++++ noao/astutil/pdm/pdmdplot.x | 101 +++++++++++ noao/astutil/pdm/pdmfindmin.x | 57 ++++++ noao/astutil/pdm/pdmfitphase.x | 43 +++++ noao/astutil/pdm/pdmgdata.x | 136 +++++++++++++++ noao/astutil/pdm/pdmminmaxp.x | 43 +++++ noao/astutil/pdm/pdmopen.x | 47 +++++ noao/astutil/pdm/pdmphase.x | 72 ++++++++ noao/astutil/pdm/pdmpplot.x | 120 +++++++++++++ noao/astutil/pdm/pdmranperm.x | 56 ++++++ noao/astutil/pdm/pdmshow.x | 56 ++++++ noao/astutil/pdm/pdmsignif.x | 61 +++++++ noao/astutil/pdm/pdmsort.x | 20 +++ noao/astutil/pdm/pdmstats.x | 37 ++++ noao/astutil/pdm/pdmtheta.x | 120 +++++++++++++ noao/astutil/pdm/pdmthetaran.x | 118 +++++++++++++ noao/astutil/pdm/pdmtplot.x | 139 +++++++++++++++ noao/astutil/pdm/pdmundelete.x | 124 +++++++++++++ noao/astutil/pdm/t_pdm.x | 72 ++++++++ 31 files changed, 2710 insertions(+) create mode 100644 noao/astutil/pdm/README create mode 100644 noao/astutil/pdm/TODO create mode 100644 noao/astutil/pdm/mkpkg create mode 100644 noao/astutil/pdm/pdm.h create mode 100644 noao/astutil/pdm/pdmalltheta.x create mode 100644 noao/astutil/pdm/pdmampep.x create mode 100644 noao/astutil/pdm/pdmautorang.x create mode 100644 noao/astutil/pdm/pdmbatch.x create mode 100644 noao/astutil/pdm/pdmclose.x create mode 100644 noao/astutil/pdm/pdmcolon.x create mode 100644 noao/astutil/pdm/pdmcompare.x create mode 100644 noao/astutil/pdm/pdmcursor.x create mode 100644 noao/astutil/pdm/pdmdelete.x create mode 100644 noao/astutil/pdm/pdmdplot.x create mode 100644 noao/astutil/pdm/pdmfindmin.x create mode 100644 noao/astutil/pdm/pdmfitphase.x create mode 100644 noao/astutil/pdm/pdmgdata.x create mode 100644 noao/astutil/pdm/pdmminmaxp.x create mode 100644 noao/astutil/pdm/pdmopen.x create mode 100644 noao/astutil/pdm/pdmphase.x create mode 100644 noao/astutil/pdm/pdmpplot.x create mode 100644 noao/astutil/pdm/pdmranperm.x create mode 100644 noao/astutil/pdm/pdmshow.x create mode 100644 noao/astutil/pdm/pdmsignif.x create mode 100644 noao/astutil/pdm/pdmsort.x create mode 100644 noao/astutil/pdm/pdmstats.x create mode 100644 noao/astutil/pdm/pdmtheta.x create mode 100644 noao/astutil/pdm/pdmthetaran.x create mode 100644 noao/astutil/pdm/pdmtplot.x create mode 100644 noao/astutil/pdm/pdmundelete.x create mode 100644 noao/astutil/pdm/t_pdm.x (limited to 'noao/astutil/pdm') diff --git a/noao/astutil/pdm/README b/noao/astutil/pdm/README new file mode 100644 index 00000000..b6102f37 --- /dev/null +++ b/noao/astutil/pdm/README @@ -0,0 +1,9 @@ +ASTUTIL.PDM -- This directory contains the sources for the PDM program. + +PDM finds periods in light curve data by a method called +Phase Dispersion Minimization. + +Reference: + +Stellingwerf, R. F., 1978, "Period Determination by Phase Dispersion + Minimization", The Astrophysical Journal, 224, pp. 953-960. diff --git a/noao/astutil/pdm/TODO b/noao/astutil/pdm/TODO new file mode 100644 index 00000000..4e668a13 --- /dev/null +++ b/noao/astutil/pdm/TODO @@ -0,0 +1,4 @@ +The following are things which would be nice to add. + +- Compute RMS in theta plots +- Text output of data points for all plots diff --git a/noao/astutil/pdm/mkpkg b/noao/astutil/pdm/mkpkg new file mode 100644 index 00000000..aeb97b97 --- /dev/null +++ b/noao/astutil/pdm/mkpkg @@ -0,0 +1,40 @@ +# PDM -- Make the Phase Dispersion Minimization package. + +update: + $checkout libpkg.a ../ + $update libpkg.a + $checkin libpkg.a ../ + ; + + +libpkg.a: + pdmalltheta.x pdm.h + pdmampep.x pdm.h + pdmautorang.x pdm.h + pdmbatch.x pdm.h + pdmclose.x pdm.h + pdmcolon.x pdm.h + pdmcompare.x pdm.h + pdmcursor.x pdm.h \ + + pdmdelete.x pdm.h + pdmdplot.x pdm.h + pdmfindmin.x pdm.h + pdmfitphase.x pdm.h + pdmgdata.x pdm.h + pdmminmaxp.x pdm.h + pdmopen.x pdm.h + pdmphase.x pdm.h + pdmpplot.x pdm.h + pdmranperm.x pdm.h + pdmshow.x pdm.h + pdmsignif.x pdm.h + pdmsort.x pdm.h + pdmstats.x pdm.h + pdmtheta.x pdm.h + pdmthetaran.x pdm.h + pdmtplot.x pdm.h + pdmundelete.x pdm.h + t_pdm.x pdm.h + ; + diff --git a/noao/astutil/pdm/pdm.h b/noao/astutil/pdm/pdm.h new file mode 100644 index 00000000..7ad38d1e --- /dev/null +++ b/noao/astutil/pdm/pdm.h @@ -0,0 +1,77 @@ +# The PDM data structure and other definitions + +define PDM_LENSTRUCT 49 # Length of PDM structure + +# Double precision. +define PDM_PMIN Memd[P2D($1)] # Minimum period to search +define PDM_PMAX Memd[P2D($1+2)] # Maximum period to search +define PDM_FMIN Memd[P2D($1+4)] # Minimum frequency to search +define PDM_FMAX Memd[P2D($1+6)] # Maximum frequency to search +define PDM_MINR Memd[P2D($1+8)] # Period to remember (min) +define PDM_NSIGMA Memd[P2D($1+10)] # num std dev. for range break +define PDM_SUMSQ Memd[P2D($1+12)] # Sum of squares of the data +define PDM_DVAR Memd[P2D($1+14)] # Variance (s ** 2) of the data +define PDM_AMPL Memd[P2D($1+16)] # Amplitude of light curve +define PDM_EPOCH Memd[P2D($1+18)] # Epoch of first maxima in data + +# Pointers +define PDM_ICD Memi[$1+20] # ICFIT pointer for data fit +define PDM_ICP Memi[$1+21] # ICFIT pointer for phasecurve fit +define PDM_CVD Memi[$1+22] # CURFIT pointer for data fit +define PDM_CVP Memi[$1+23] # CURFIT pointer for phasecurve fit +define PDM_GP Memi[$1+24] # PDM graphics GIO pointer +define PDM_LFD Memi[$1+25] # Log file descriptor +define PDM_PFD Memi[$1+26] # Plot file descriptor +define PDM_GT Memi[$1+27] # PDM gtools pointer +define PDM_XP Memi[$1+28] # Pointer to data ordinates +define PDM_ODYP Memi[$1+29] # Pointer to original data abscissas +define PDM_DYP Memi[$1+30] # Pointer to working data abscissas +define PDM_INUSEP Memi[$1+31] # Pointer to in-use array +define PDM_ERRP Memi[$1+32] # Pointer to error bar array +define PDM_XTHP Memi[$1+33] # Pointer to theta plot ordinates +define PDM_YTHP Memi[$1+34] # Pointer to theta plot abscissas +define PDM_XPHP Memi[$1+35] # Pointer to phasecurve plot ordinates +define PDM_YPHP Memi[$1+36] # Pointer to phasecurve plot abscissas +define PDM_PHERRP Memi[$1+37] # Pointer to phasecurve plot errors +define PDM_SORTP Memi[$1+38] # Pointer to array defining sort +define PDM_SAMPLEP Memi[$1+39] # Pointer to sample (range) string +define PDM_RG Memi[$1+40] # Pointer to range (sample) structure + +# Constants +define PDM_NPT Memi[$1+41] # Number of data points +define PDM_NTHPT Memi[$1+42] # Number of theta points +define PDM_NRANGE Memi[$1+43] # Number of ranges. + +# Other +define PDM_RESID Memi[$1+44] # Using residuals? flag +define PDM_RANGE Memi[$1+45] # Using ranges? flag +define PDM_DEBUG Memb[$1+46] # Debug? flag +define PDM_PLUSPOINT Memi[$1+47] # Threshold data number to use plus +define PDM_EB Memi[$1+48] # Use error bars? flag + +# Macro definitions +define PDM_X Memd[PDM_XP($1)+$2-1] # data ordinates +define PDM_ODY Memd[PDM_ODYP($1)+$2-1] # orig data abscissas +define PDM_DY Memd[PDM_DYP($1)+$2-1] # working data abscissas +define PDM_ERR Memr[PDM_ERRP($1)+$2-1] # error bars +define PDM_INUSE Memi[PDM_INUSEP($1)+$2-1] # in-use array +define PDM_XTH Memd[PDM_XTHP($1)+$2-1] # theta plot ordinates +define PDM_YTH Memd[PDM_YTHP($1)+$2-1] # theta plot abscissas +define PDM_XPH Memd[PDM_XPHP($1)+$2-1] # phase plot ordinates +define PDM_YPH Memd[PDM_YPHP($1)+$2-1] # phase plot abscissas +define PDM_PHERR Memr[PDM_PHERRP($1)+$2-1] # phase plot errors +define PDM_SORT Memi[PDM_SORTP($1)+$2-1] # array defining sort +define PDM_SAMPLE Memc[PDM_SAMPLEP($1)] # sample (range) string + + +# Plot types. (ptype) +define DATAPLOT 0 # data plot +define THETAPPLOT 1 # theta period plot +define THETAFPLOT 2 # theta frequency plot +define PHASEPLOT 3 # phase plot + +define MAX_RANGES 20 # maximum number of range segments +define PDM_SZ_TITLE (4*SZ_LINE) # size of pdm plot title buffer +define BIN10 100 # numpts for bin 5/10 split + # otherwise, plot pluses +define HELP "noao$lib/scr/pdm.key" # where help file is diff --git a/noao/astutil/pdm/pdmalltheta.x b/noao/astutil/pdm/pdmalltheta.x new file mode 100644 index 00000000..b7a48889 --- /dev/null +++ b/noao/astutil/pdm/pdmalltheta.x @@ -0,0 +1,104 @@ +include +include +include +include "pdm.h" + +# PDM_ALLTHETA -- Calculate the theta statistic for the period window. + +procedure pdm_alltheta (pdmp, porf) + +pointer pdmp # structure pointer +int porf # period or frequency flag + +int i +double factor, period, pdm_theta(), mmean, mmin, mmax, scale +pointer rg, rg_xrangesd() +errchk calloc, realloc, rg_xrangesd + +begin + # Check to see that maxp > minp. + if (PDM_PMAX(pdmp) <= PDM_PMIN(pdmp)) + call error (0,"maxp smaller or equal minp in alltheta") + + # Calculate the mean of the Y data, remove it. + mmean = 0.0d+0 + do i = 1, PDM_NPT(pdmp) { + mmean = mmean + PDM_DY(pdmp,i) + } + mmean = mmean/PDM_NPT(pdmp) + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = PDM_DY(pdmp,i) - mmean + } + + # Find the min and max of the data, scale the data 0 to 1. + mmin = PDM_DY(pdmp,1) + mmax = PDM_DY(pdmp,1) + do i = 1, PDM_NPT(pdmp) { + if (PDM_DY(pdmp,i) > mmax) + mmax = PDM_DY(pdmp,i) + if (PDM_DY(pdmp,i) < mmin) + mmin = PDM_DY(pdmp,i) + } + scale = mmax - mmin + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = (PDM_DY(pdmp,i) - mmin)/scale + } + + # Bring the statistics up to date. + iferr (call pdm_statistics (pdmp)) + call error (0, "alltheta: error in statistics") + + # Allocate space for the ordinate and abscissa vectors for theta + # in the pdm data structure. + + if (PDM_XTHP(pdmp) == NULL) { + call calloc (PDM_XTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE) + call calloc (PDM_YTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE) + } else { + call realloc (PDM_XTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE) + call realloc (PDM_YTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE) + } + + # Calculate constant for incrementing the period. Give equally + # spaced frequencies. + + if (PDM_PMIN(pdmp) <= EPSILOND) + if (PDM_NTHPT(pdmp) >= 1) + PDM_PMIN(pdmp) = PDM_PMAX(PDMP)/PDM_NTHPT(pdmp) + else + call error (1, "alltheta: num thpts < 1") + + if (PDM_PMIN(pdmp) >= EPSILONR && PDM_NTHPT(pdmp) != 1) + factor = (PDM_PMAX(pdmp)/PDM_PMIN(pdmp))** + (1.0/(PDM_NTHPT(pdmp)-1)) + + # Calculate the ranges information from the sample string. + rg = rg_xrangesd (PDM_SAMPLE(pdmp), PDM_X(pdmp,1), PDM_NPT(pdmp)) + + # Call pdm_theta for each period and store the thetas and periods + # in the pdm data structure, calculate frequencies (1/p) as we go. + + period = PDM_PMIN(pdmp) + do i = 1, PDM_NTHPT(pdmp) { + PDM_YTH(pdmp,i) = pdm_theta (pdmp, rg, period) + if (porf == THETAPPLOT) + PDM_XTH(pdmp,i) = period + else { + if (period > EPSILOND) + PDM_XTH(pdmp,i) = 1.0d+0/period + else + call error (0, "alltheta: period very close to zero") + } + period = period * factor + } + + # Remove the scaling. + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = (PDM_DY(pdmp,i) * scale) + mmin + } + + # Put the mean back in. + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = PDM_DY(pdmp,i) + mmean + } +end diff --git a/noao/astutil/pdm/pdmampep.x b/noao/astutil/pdm/pdmampep.x new file mode 100644 index 00000000..08c2c4e1 --- /dev/null +++ b/noao/astutil/pdm/pdmampep.x @@ -0,0 +1,38 @@ +include +include +include +include "pdm.h" + +# PDM_AMPEP -- Calculate the amplitude and epoch for this data. + +procedure pdm_ampep (pdmp, period) + +pointer pdmp # PDM structure pointer +double period # period for which to calculate + +int i, isave +double npt, ymin, ymax +errchk pdm_phase + +begin + npt = PDM_NPT(pdmp) + + # Find the maximum and minimum values in the data. + # The difference is the amplitude. + + ymax = -MAX_DOUBLE + ymin = MAX_DOUBLE + do i = 1, npt { + if (PDM_INUSE(pdmp,i) == 0) + next + if (PDM_DY(pdmp,i) < ymin) + ymin = PDM_DY(pdmp,i) + if (PDM_DY(pdmp,i) > ymax) { + ymax = PDM_DY(pdmp,i) + isave = i + } + } + + PDM_AMPL(pdmp) = ymax - ymin + PDM_EPOCH(pdmp) = PDM_X(pdmp, isave) +end diff --git a/noao/astutil/pdm/pdmautorang.x b/noao/astutil/pdm/pdmautorang.x new file mode 100644 index 00000000..ba27d628 --- /dev/null +++ b/noao/astutil/pdm/pdmautorang.x @@ -0,0 +1,101 @@ +include +include +include +include +include "pdm.h" + +# PDM_AUTORANG -- Calculate the ranges division of the data automatically. + +int procedure pdm_autorang (pdmp) + +pointer pdmp # PDM structure pointer + +int npt, i, nrng +double sumsq, sum, var, sd, maxdif, meansep, rbegin, rend +int rngstrt +pointer sep, command, sp + +begin + call smark (sp) + call salloc (command, SZ_LINE, TY_CHAR) + + npt = PDM_NPT(pdmp) + + # Calculate the mean and standard deviation of the x-axis + # separation of the data points. (time intervals) + # Allocate an array of separations and fill it. + + call salloc (sep, npt-1, TY_DOUBLE) + do i = 1, npt-1 { + Memd[sep+i-1] = PDM_X(pdmp,i+1) - PDM_X(pdmp,i) + } + + sumsq = 0.0 + sum = 0.0 + for (i=1; i<=(npt-1); i=i+1) { + sumsq = sumsq + Memd[sep+i-1]**2 # Sum of squares. + sum = sum + Memd[sep+i-1] + } + if (npt != 1) { + var = (sumsq - sum**2/npt)/(npt - 1) # Variance. + sd = var**.5 # Standard Deviation. + } + + # Mean separation, maximum time diff. + if (npt != 1) { + meansep = (PDM_X(pdmp,npt) - PDM_X(pdmp,1)) / double(npt-1) + maxdif = meansep + sd * PDM_NSIGMA(pdmp) + } + + # Look through the separations and if we find one that is more + # than nsigma away from the mean on the plus side, divide the + # data at this point into another range. + + nrng = 0 + rngstrt = 1 + PDM_SAMPLE(pdmp) = EOS + do i = 1, npt - 1 { + if (Memd[sep+i-1] > maxdif) { + nrng = nrng + 1 + if (nrng > MAX_RANGES) { + call sfree (sp) + call error (0,"Max num ranges exceeded in autorange") + break + } + + # End of last range = x(i) + # If (i+1 != npts) beginning of next range = x(i+1) + # Remember where the next range starts. + + rbegin = PDM_X(pdmp,rngstrt) + rend = PDM_X(pdmp,i) + if ((i+1) < npt) + rngstrt = i+1 + + # Sprintf range info at end of string. + call sprintf (Memc[command], SZ_LINE, " %g:%g") + call pargd (rbegin) + call pargd (rend) + call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE) + } + } + + # Finish up last range if needed. + If (rngstrt <= npt && rngstrt > 1) { + rbegin = PDM_X(pdmp,rngstrt) + rend = PDM_X(pdmp,i) + + # Sprintf range info at end of string. + call sprintf (Memc[command], SZ_LINE, " %g:%g") + call pargd (rbegin) + call pargd (rend) + call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE) + } + + # If no ranges found, set the sample string to '*', all data. + if (nrng == 0) + call sprintf (PDM_SAMPLE(pdmp), SZ_LINE, "*") + + call sfree (sp) + return (nrng) +end diff --git a/noao/astutil/pdm/pdmbatch.x b/noao/astutil/pdm/pdmbatch.x new file mode 100644 index 00000000..6828b783 --- /dev/null +++ b/noao/astutil/pdm/pdmbatch.x @@ -0,0 +1,49 @@ +include +include +include +include "pdm.h" + +# PDM_BATCH -- Batch mode calculation. + +procedure pdm_batch (pdmp, file, infile, flip) + +pointer pdmp # PDM structure pointer +char file[SZ_FNAME] # file to store information +char infile[SZ_FNAME] # input data file name +bool flip # flip y-axis scale + +int fd +double pdm_signif(), signif +bool verbose +errchk pdm_alltheta, pdm_signif, pdm_ampep, pdm_phase + +begin + # Plot the data. + call pdm_dplot (pdmp, infile, flip) + + # Call pdm_alltheta to get the theta array; plot to metafile. + call pdm_alltheta (pdmp, THETAPPLOT) + call pdm_tplot (pdmp, THETAPPLOT, infile) + + # Call pdm_signif to calculate the significance of the theta minimum. + signif = pdm_signif (pdmp, PDM_MINR(pdmp)) + + # Write this significance out to the file. + fd = PDM_LFD(pdmp) + call fprintf (fd, "significance at this minimum = %g\n") + call pargd (signif) + call close (fd) + + # Call pdm_amplitudeepoch to get the amplitude and epoch. + call pdm_ampep (pdmp, PDM_MINR(pdmp)) + + # Call pdm_phase to get the phase curve. + # Plot to metafile. + + call pdm_phase(pdmp, PDM_MINR(pdmp), PDM_EPOCH(pdmp)) + call pdm_pplot (pdmp, PDM_MINR(pdmp), infile, flip) + + # Call pdm_show to print the parameter information to a file. + verbose = TRUE + call pdm_show (pdmp, file, verbose) +end diff --git a/noao/astutil/pdm/pdmclose.x b/noao/astutil/pdm/pdmclose.x new file mode 100644 index 00000000..b7e4541e --- /dev/null +++ b/noao/astutil/pdm/pdmclose.x @@ -0,0 +1,63 @@ +include +include "pdm.h" + +# PDM_CLOSE -- Close a PDM data structure. + +procedure pdm_close (pdmp, interactive) + +pointer pdmp # PDM structure pointer +bool interactive # interactive flag + +begin + # Close icfit pointers, curfit pointers, and graphics pointers, + # and the ranges pointer. + + if (PDM_ICD(pdmp) != NULL) + call ic_closed (PDM_ICD(pdmp)) + if (PDM_ICP(pdmp) != NULL) + call ic_closed (PDM_ICP(pdmp)) + if (PDM_CVD(pdmp) != NULL) + call dcvfree (PDM_CVD(pdmp)) + if (PDM_CVP(pdmp) != NULL) + call dcvfree (PDM_CVP(pdmp)) + if (PDM_GT(pdmp) != NULL) + call gt_free (PDM_GT(pdmp)) + if (PDM_GP(pdmp) != NULL) + call gclose (PDM_GP(pdmp)) + if (PDM_RG(pdmp) != NULL) + call rg_free (PDM_RG(pdmp)) + if (PDM_LFD(pdmp) != NULL) + call close (PDM_LFD(pdmp)) + if (!interactive) + if (PDM_PFD(pdmp) != NULL) + call close (PDM_PFD(pdmp)) + + # Free the data vectors. + if (PDM_XP(pdmp) != NULL) + call mfree (PDM_XP(pdmp), TY_DOUBLE) + if (PDM_ODYP(pdmp) != NULL) + call mfree (PDM_ODYP(pdmp), TY_DOUBLE) + if (PDM_DYP(pdmp) != NULL) + call mfree (PDM_DYP(pdmp), TY_DOUBLE) + if (PDM_ERRP(pdmp) != NULL) + call mfree (PDM_ERRP(pdmp), TY_REAL) + if (PDM_INUSEP(pdmp) != NULL) + call mfree (PDM_INUSEP(pdmp), TY_INT) + if (PDM_XTHP(pdmp) != NULL) + call mfree (PDM_XTHP(pdmp), TY_DOUBLE) + if (PDM_YTHP(pdmp) != NULL) + call mfree (PDM_YTHP(pdmp), TY_DOUBLE) + if (PDM_XPHP(pdmp) != NULL) + call mfree (PDM_XPHP(pdmp), TY_DOUBLE) + if (PDM_YPHP(pdmp) != NULL) + call mfree (PDM_YPHP(pdmp), TY_DOUBLE) + if (PDM_PHERRP(pdmp) != NULL) + call mfree (PDM_PHERRP(pdmp), TY_REAL) + if (PDM_SORTP(pdmp) != NULL) + call mfree (PDM_SORTP(pdmp), TY_INT) + if (PDM_SAMPLEP(pdmp) != NULL) + call mfree (PDM_SAMPLEP(pdmp), TY_CHAR) + + # Free the pdm data structure. + call mfree (pdmp, TY_STRUCT) +end diff --git a/noao/astutil/pdm/pdmcolon.x b/noao/astutil/pdm/pdmcolon.x new file mode 100644 index 00000000..6d438515 --- /dev/null +++ b/noao/astutil/pdm/pdmcolon.x @@ -0,0 +1,292 @@ +include +include +include +include +include "pdm.h" + +define PDM_KEYWORDS "|show|vshow|minp|maxp|minf|maxf|ntheta|sample|signif\ + |ampep|phase|unreject|alldata|origdata|" + +define SHOW 1 # Show parameter settings +define VSHOW 2 # Show verbose information +define PMIN 3 # Set min search period +define PMAX 4 # Set max search period +define FMIN 5 # Set min search frequency +define FMAX 6 # Set max search frequency +define NTHETA 7 # Set number of points for theta +define SAMPLE 8 # Set/show the sample ranges +define SIGNIF 9 # Find theta significance +define AMPEP 10 # Amplitude and Epoch +define PHASE 11 # Graph phase curve +define UNREJECT 12 # Unreject all the rejected data points +define ALLDATA 13 # Reset range to entire dataset +define ORIGDATA 14 # Reset data to origional dataset + +define SLOWTHRESH 500 # Threshold above which theta calc gets slow + +# PDM_COLON -- Decode colon commands. + +procedure pdm_colon (pdmp, cmdstr, ptype, infile, period, flip) + +pointer pdmp # PDM structure pointer +char cmdstr[ARB] # Command string +int ptype # plot type +char infile[SZ_FNAME] # input file name +double period # current working period +bool flip # flip the y-axis scale + +int nscan(), strdic() +int itemp, i +double temp, p1, signif, pdm_signif() +bool verbose +pointer cmd, sp +errchk pdm_signif, pdm_ampep, pdm_phase + +string keywords PDM_KEYWORDS + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + + # Unpack keyword from string, look up command in dictionary. + # Switch on command. Call the appropriate subroutines. + + switch (strdic (Memc[cmd], Memc[cmd], SZ_FNAME, keywords)) { + case SHOW: + # Show parameter settings. + call gargstr (Memc[cmd], SZ_LINE) + verbose = false + if (Memc[cmd] == EOS) { + call gdeactivate (PDM_GP(pdmp), AW_CLEAR) + iferr (call pdm_show (pdmp, "STDOUT", verbose)) + call erract (EA_WARN) + call greactivate (PDM_GP(pdmp), AW_PAUSE) + } else { + iferr (call pdm_show (pdmp, Memc[cmd], verbose)) + call erract (EA_WARN) + } + case VSHOW: + # Show verbose information. + call gargstr (Memc[cmd], SZ_LINE) + verbose = true + if (Memc[cmd] == EOS) { + call gdeactivate (PDM_GP(pdmp), AW_CLEAR) + iferr (call pdm_show (pdmp, "STDOUT", verbose)) + call erract (EA_WARN) + call greactivate (PDM_GP(pdmp), AW_PAUSE) + } else { + iferr (call pdm_show (pdmp, Memc[cmd], verbose)) + call erract (EA_WARN) + } + case PMIN: + # List or set minimum period. + call gargd (temp) + if (nscan() == 2) { + # Set the period minimum in structure. + PDM_PMIN(pdmp) = temp + if (temp >= EPSILOND) + PDM_FMAX(pdmp) = 1.0d+0/temp + # Save this minp out to the parameter. + call clputd ("minp", temp) + } else { + # Print out the period minimum from the structure. + call printf ("Current minimum period = %g\n") + call pargd (PDM_PMIN(pdmp)) + call flush (STDOUT) + } + case PMAX: + # List or set maximum period. + call gargd (temp) + if (nscan() == 2) { + # Set the period maximum in structure. + PDM_PMAX(pdmp) = temp + if (temp >= EPSILOND) + PDM_FMIN(pdmp) = 1.0d+0/temp + # Save this minp out to the parameter. + call clputd ("maxp", temp) + } else { + # Print out the period maximum from the structure. + call printf ("Current maximum period = %g\n") + call pargd (PDM_PMAX(pdmp)) + call flush (STDOUT) + } + case FMIN: + # List or set minimum frequency. + call gargd (temp) + if (nscan() == 2) { + # Set the frequency minimum in structure. + PDM_FMIN(pdmp) = temp + if (temp >= EPSILOND) + PDM_PMAX(pdmp) = 1.0d+0/temp + # Save this minp out to the parameter. + if (temp >= EPSILOND) + call clputd ("maxp", 1.0d+0/temp) + } else { + # Print out the frequency minimum from the structure. + call printf ("Current minimum frequency = %g\n") + call pargd (PDM_FMIN(pdmp)) + call flush (STDOUT) + } + case FMAX: + # List or set maximum frequency. + call gargd (temp) + if (nscan() == 2) { + # Set the frequency maximum in structure. + PDM_FMAX(pdmp) = temp + if (temp >= EPSILOND) + PDM_PMIN(pdmp) = 1.0d+0/temp + # Save this minp out to the parameter. + if (temp >= EPSILOND) + call clputd ("minp", 1.0d+0/temp) + } else { + # Print out the frequency maximum from the structure. + call printf ("Current maximum frequency = %g\n") + call pargd (PDM_FMAX(pdmp)) + call flush (STDOUT) + } + case NTHETA: + # Set/show number of theta points + call gargi (itemp) + if (nscan() == 2) { + # Set ntheta in structure. + PDM_NTHPT(pdmp) = itemp + if (itemp > SLOWTHRESH) + # Give message saying that with this ntheta, the + # theta calculation will take quite a while. + call printf ("Large ntheta => long calculation time \007\n") + } else { + # Print out the value of ntheta from the structure. + call printf ("Number of theta points = %g\n") + call pargi (PDM_NTHPT(pdmp)) + call flush (STDOUT) + } + case SAMPLE: + # List or set the sample points. + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call printf ("sample = %s\n") + call pargstr (PDM_SAMPLE(pdmp)) + } else { + if (ptype == DATAPLOT) { + call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp), + PDM_X(pdmp,1), PDM_NPT(pdmp), 0) + } + call strcpy (Memc[cmd], PDM_SAMPLE(pdmp), SZ_LINE) + if (ptype == DATAPLOT) { + call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp), + PDM_X(pdmp,1), PDM_NPT(pdmp), 1) + } + } + case SIGNIF: + # Calculate the significance of theta at period. + p1 = 0.0 + call gargd (temp) + if (nscan() == 2) # User entered a period. + p1 = temp + else { + # Use remembered period. + if (PDM_MINR(pdmp) >= EPSILOND) + p1 = PDM_MINR(pdmp) + else { + call printf ("No remembered minimum period. \007\n") + call flush (STDOUT) + } + } + # Calculate significance at cursor x position (per). + if (p1 >= EPSILOND) { + signif = pdm_signif (pdmp, p1) + # Print at bottom of screen. + call printf ("Significance at period %g = %g\n") + call pargd (p1) + call pargd (signif) + call flush (STDOUT) + } + case AMPEP: + # Calculate the amplitude and epoch for the data. + p1 = 0.0 + call gargd (temp) + if (nscan() == 2) # User entered a period. + p1 = temp + else { + # Use remembered period. + if (PDM_MINR(pdmp) >= EPSILOND) + p1 = PDM_MINR(pdmp) + else { + call printf ("No remembered minimum period. \007\n") + call flush (STDOUT) + } + } + # Calculate ampl & epoch at p1. + if (p1 >= EPSILOND) { + call pdm_ampep (pdmp, p1) + # Print at bottom of screen. + call printf("amplitude of data at period %g = %g, epoch = %g\n") + call pargd (p1) + call pargd (PDM_AMPL(pdmp)) + call pargd (PDM_EPOCH(pdmp)) + call flush (STDOUT) + } + case PHASE: + # Phase curve plot. + call gargd (temp) + if (nscan() == 2) { + # Calculate the phase curve, then make the plot. + call pdm_ampep (pdmp, temp) + call pdm_phase(pdmp, temp, PDM_EPOCH(pdmp)) + call pdm_pplot (pdmp, temp, infile, flip) + period = temp + } else { + # Use remembered period. + if (PDM_MINR(pdmp) >= EPSILOND) { + call pdm_ampep (pdmp, PDM_MINR(pdmp)) + call pdm_phase(pdmp, PDM_MINR(pdmp), PDM_EPOCH(pdmp)) + call pdm_pplot (pdmp, PDM_MINR(pdmp), infile, flip) + } else { + call printf ("No remembered minimum period. \007\n") + call flush (STDOUT) + } + } + ptype = PHASEPLOT + case UNREJECT: + # Copy original data vector into working data vector and + # set inuse array to all ones. + + do i = 1, PDM_NPT(pdmp) { + PDM_INUSE(pdmp,i) = 1 + PDM_DY(pdmp,i) = PDM_ODY(pdmp,i) + } + PDM_RESID(pdmp) = NO + + # Check type of plot and replot. + if (ptype == DATAPLOT) + call pdm_dplot (pdmp, infile, flip) + if (ptype == PHASEPLOT) + call pdm_pplot (pdmp, period, infile, flip) + case ALLDATA: + # Initialize the sample string and erase from the graph. + if (ptype == DATAPLOT) { + call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp), + PDM_X(pdmp,1), PDM_NPT(pdmp), 0) + } + call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE) + case ORIGDATA: + # Copy the original data vector into the working data vector. + do i = 1, PDM_NPT(pdmp) + PDM_DY(pdmp,i) = PDM_ODY(pdmp,i) + PDM_RESID(pdmp) = NO + + # Replot data. + call pdm_dplot (pdmp, infile, flip) + ptype = DATAPLOT + default: + # Error, unknown colon command; ring bell. + call printf ("\007") + call printf ("v/show minp/f maxp/f ntheta sample ") + call printf ("signif ampep phase unreject all/origdata\n") + } + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmcompare.x b/noao/astutil/pdm/pdmcompare.x new file mode 100644 index 00000000..b2d60978 --- /dev/null +++ b/noao/astutil/pdm/pdmcompare.x @@ -0,0 +1,25 @@ +include +include "pdm.h" + +# Compare procedure for qsort. + +int procedure pdm_compare (item1, item2) + +int item1 # index of first phase +int item2 # index of second phase + +pointer comarray +common /sortcom/ comarray +double p1, p2 + +begin + p1 = Memd[comarray+item1-1] + p2 = Memd[comarray+item2-1] + + if (p1 > p2) + return (1) + if (p1 == p2) + return (0) + if (p1 < p2) + return (-1) +end diff --git a/noao/astutil/pdm/pdmcursor.x b/noao/astutil/pdm/pdmcursor.x new file mode 100644 index 00000000..067ce966 --- /dev/null +++ b/noao/astutil/pdm/pdmcursor.x @@ -0,0 +1,383 @@ +include +include +include +include +include +include +include "pdm.h" + +define PROMPT "pdm options" + +# PDM_CURSOR -- Get the next command from the user in a graphics cursor loop. +# Perform the requested function. + +procedure pdm_cursor (pdmp, ptype, infile, flip) + +pointer pdmp # pointer to PDM data structure +int ptype # type of plot on the screen +char infile[SZ_FNAME] # input file name +bool flip # flip the y-axis scale + +real xx, yy +double x, y # cursor coordinates +double xmax, xmin +double period +int wcs # wcs to which coordinates belong +int key # keystroke value of cursor event +int stridxs(), gqverify(), ier +int clgcur(), i +int pdm_delete(), index +int pdm_undelete(), pdm_findmin() +double rx1, rx2 +double dcveval() +double pdm_signif(), signif +pointer weights # pointer to temporary weights array for icfit +pointer command # string value, if any +pointer sp, sptemp +errchk pdm_colon, icg_fit, pdm_fitphase, pdm_alltheta +errchk pdm_phase, pdm_signif, pdm_ampep, pdm_findmin + +begin + call smark (sp) + call salloc (command, SZ_LINE, TY_CHAR) + + while (clgcur ("cursor", xx, yy, wcs, key, Memc[command], + SZ_LINE) != EOF) { + + x = double(xx) + y = double(yy) + + # Switch on command, take appropriate action. + switch (key) { + case '?': + # List options. + call gpagefile (PDM_GP(pdmp), HELP, PROMPT) + case ':': + # Colon command. + call pdm_colon (pdmp, Memc[command], ptype, infile, + period, flip) + case 'h': + # Graph the data. + call pdm_dplot (pdmp, infile, flip) + ptype = DATAPLOT + case 'f': + # Call icfit on the data. + if (ptype == DATAPLOT) { + # Set the min/max ordinate values for icfit. + call alimd (PDM_X(pdmp,1), PDM_NPT(pdmp), xmin, xmax) + call ic_putr (PDM_ICD(pdmp), "xmin", real(xmin)) + call ic_putr (PDM_ICD(pdmp), "xmax", real(xmax)) + + # Allocate a temporary weights array and fill it with the + # in-use array values. + + call smark (sptemp) + call salloc (weights, PDM_NPT(pdmp), TY_DOUBLE) + do i = 1, PDM_NPT(pdmp) + Memd[weights+i-1] = double(PDM_INUSE(pdmp,i)) + + # Call icfit. + if (PDM_NPT(pdmp) >= 2) + call icg_fitd (PDM_ICD(pdmp), PDM_GP(pdmp), "cursor", + PDM_GT(pdmp), PDM_CVD(pdmp), PDM_X(pdmp,1), + PDM_DY(pdmp,1), Memd[weights], PDM_NPT(pdmp)) + + # Recover the weights array back into the in-use array. + do i = 1, PDM_NPT(pdmp) + PDM_INUSE(pdmp,i) = int(Memd[weights+i-1]) + + call sfree (sptemp) + + # Replot. + call pdm_dplot (pdmp, infile, flip) + call sfree (sptemp) + } else if (ptype == PHASEPLOT) { + # Call icfit on the phases. + call pdm_fitphase (pdmp) + + # Replot. + call pdm_pplot (pdmp, period, infile, flip) + } else + call printf ("Can't fit a THETA plot \007\n") + case 'i': + # Theta vs. frequency plot. + if (PDM_PMIN(pdmp) < EPSILOND && PDM_PMAX(pdmp) < EPSILOND) + call pdm_minmaxp(pdmp) + + # Calculate theta. + call pdm_alltheta (pdmp, THETAFPLOT) + + # Graph theta vs frequency. + call pdm_tplot (pdmp, THETAFPLOT, infile) + ptype = THETAFPLOT + case 'k': + # Theta vs. period plot. + if (PDM_PMIN(pdmp) < EPSILOND && PDM_PMAX(pdmp) < EPSILOND) + call pdm_minmaxp(pdmp) + + # Calculate theta. + call pdm_alltheta (pdmp, THETAPPLOT) + + # Graph theta vs frequency. + call pdm_tplot (pdmp, THETAPPLOT, infile) + ptype = THETAPPLOT + case 'p': + # Phase curve plot. + if (ptype == THETAPPLOT) { + # Find the epoch, calculate the phases w.r.t. this epoch. + call pdm_ampep (pdmp, x) + call pdm_phase(pdmp, x, PDM_EPOCH(pdmp)) + call pdm_pplot (pdmp, x, infile, flip) + ptype = PHASEPLOT + period = x + } else if (ptype == THETAFPLOT) { + # Find the epoch, calculate the phases w.r.t. this epoch. + call pdm_ampep (pdmp, 1.0d+0/x) + call pdm_phase(pdmp, 1.0d+0/x, PDM_EPOCH(pdmp)) + call pdm_pplot (pdmp, 1.0d+0/x, infile, flip) + ptype = PHASEPLOT + period = 1.0d+0/x + } else + call printf ("Wrong type of plot on screen for p key\007\n") + case 'd': + # Delete the point nearest the cursor. + index = pdm_delete (pdmp, x, y, ptype) + case 'u': + # Undelete the point nearest the cursor. + index = pdm_undelete (pdmp, x, y, ptype) + case 'j': + # Subtract the fit from the data and use residuals. + if (PDM_CVD(pdmp) == NULL) + call printf ("Fit has not been done. \007\n") + else if (PDM_RESID(pdmp) == YES) + call printf ("Already using residuals. \007\n") + else { + # For each point, calculate the fit function and subtract + # it from the data. + + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = PDM_DY(pdmp,i) - + dcveval (PDM_CVD(pdmp), PDM_X(pdmp,i)) + } + PDM_RESID(pdmp) = YES + if (ptype == DATAPLOT) + call pdm_dplot (pdmp, infile, flip) + } + case 's': + # Set sample regions with the cursor. + if (ptype == DATAPLOT) { + if (stridxs ("*", PDM_SAMPLE(pdmp)) > 0) + PDM_SAMPLE(pdmp) = EOS + + rx1 = x + call printf ("again:\n") + if (clgcur ("cursor", xx, yy, wcs, key, Memc[command], + SZ_LINE) == EOF) + break + rx2 = double(xx) + + call sprintf (Memc[command], SZ_LINE, " %g:%g") + call pargd (rx1) + call pargd (rx2) + + call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE) + call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp), + PDM_X(pdmp,1), PDM_NPT(pdmp), 1) + call printf (" \n") + call gflush (PDM_GP(pdmp)) + } else + call printf ("Wrong type of plot for s key \007\n") + case 't': + # Initialize the sample string and erase from the graph. + if (ptype == DATAPLOT) { + call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp), + PDM_X(pdmp,1), PDM_NPT(pdmp), 0) + } + call gflush (PDM_GP(pdmp)) + call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE) + call gflush (PDM_GP(pdmp)) + case 'g': + # Significance of theta at cursor x position. + if (ptype == THETAPPLOT) { + # Calculate significance at cursor x position (per). + signif = pdm_signif (pdmp, x) + # Print at bottom of screen. + call printf ("Significance at cursor = %g\n") + call pargd (signif) + } else if (ptype == THETAFPLOT) { + # Calculate significance at cursor x position (per). + signif = pdm_signif (pdmp, 1.0d+0/x) + # Print at bottom of screen. + call printf ("Significance at cursor = %g\n") + call pargd (signif) + } else if (ptype == PHASEPLOT) { + # Calculate significance at current period/frequency + signif = pdm_signif (pdmp, period) + # Print at bottom of screen + call printf ("Significance of this period = %g\n") + call pargd (signif) + } else { + # Data plot. + call printf ("Wrong type of plot for g key \007\n") + } + case 'a': + # Amplitude and epoch at cursor x position. + if (ptype == THETAPPLOT) { + # Calculate ampl & epoch at cursor x position. + call pdm_ampep (pdmp, x) + # Print at bottom of screen. + call printf ("amplitude of data = %g, epoch = %g\n") + call pargd (PDM_AMPL(pdmp)) + call pargd (PDM_EPOCH(pdmp)) + } else if (ptype == THETAFPLOT) { + # Calculate ampl & epoch at cursor x position. + call pdm_ampep (pdmp, 1.0d+0/x) + # Print at bottom of screen. + call printf ("amplitude of data = %g, epoch = %g\n") + call pargd (PDM_AMPL(pdmp)) + call pargd (PDM_EPOCH(pdmp)) + } else if (ptype == PHASEPLOT) { + # Calculate ampl & epoch at current period/frequency + call pdm_ampep (pdmp, period) + # Print at bottom of screen. + call printf ("amplitude of data = %g, epoch = %g\n") + call pargd (PDM_AMPL(pdmp)) + call pargd (PDM_EPOCH(pdmp)) + } else { + # Data plot. + call printf ("Wrong type of plot for a key \007\n") + } + case ',': + # Set minp or minf to cursor x position. + if (ptype == THETAPPLOT) { + PDM_PMIN(pdmp) = x + PDM_FMIN(pdmp) = 1.0d+0/x + # Print at bottom of screen. + call printf ("minp now %g\n") + call pargd (PDM_PMIN(pdmp)) + } else if (ptype == THETAFPLOT) { + PDM_FMAX(pdmp) = x + PDM_PMAX(pdmp) = 1.0d+0/x + # Print at bottom of screen. + call printf ("minf now %g\n") + call pargd (PDM_FMAX(pdmp)) + } else { + # Data plot or phase plot. + call printf ("Wrong type of plot for , key \007\n") + } + case '.': + # Set maxp or maxf to cursor x position. + if (ptype == THETAPPLOT) { + PDM_PMAX(pdmp) = x + PDM_FMAX(pdmp) = 1.0d+0/x + # Print at bottom of screen. + call printf ("maxp now %g\n") + call pargd (PDM_PMAX(pdmp)) + } else if (ptype == THETAFPLOT) { + PDM_FMIN(pdmp) = x + PDM_PMIN(pdmp) = 1.0d+0/x + # Print at bottom of screen. + call printf ("maxf now %g\n") + call pargd (PDM_FMIN(pdmp)) + } else { + # Data plot or phase plot. + call printf ("Wrong type of plot for . key \007\n") + } + case 'm': + # Mark range and find minimum in this range. + if (ptype == THETAFPLOT || ptype == THETAPPLOT) { + rx1 = x + call printf ("again:\n") + if (clgcur ("cursor", xx, yy, wcs, key, Memc[command], + SZ_LINE) == EOF) + break + rx2 = double(xx) + index = pdm_findmin(pdmp, ptype, rx1, rx2, 1, + PDM_NTHPT(pdmp)) + PDM_MINR(pdmp) = PDM_XTH(pdmp,index) + call printf ("period at minimum = %g, frequency = %g\n") + call pargd (PDM_XTH(pdmp,index)) + call pargd (1.0d+0/PDM_XTH(pdmp,index)) + + } else + call printf ("Wrong type of plot for m key. \007\n") + case 'r': + # Check type of plot and replot. + if (ptype == DATAPLOT) + call pdm_dplot (pdmp, infile, flip) + if (ptype == THETAPPLOT) + call pdm_tplot (pdmp, THETAPPLOT, infile) + if (ptype == THETAFPLOT) + call pdm_tplot (pdmp, THETAFPLOT, infile) + if (ptype == PHASEPLOT) + call pdm_pplot (pdmp, period, infile, flip) + case 'e': + # Toggle error bars. + if (PDM_EB(pdmp) == NO) + PDM_EB(pdmp) = YES + else + PDM_EB(pdmp) = NO + + # Check type of plot and replot. + if (ptype == PHASEPLOT) + call pdm_pplot (pdmp, period, infile, flip) + case 'x': + # Remove a trend from the data by fitting a straight line to + # the data and removing it. + + # Set the min/max ordinate values for icfit. + call alimd (PDM_X(pdmp,1), PDM_NPT(pdmp), xmin, xmax) + call dcvinit (PDM_CVD(pdmp), SPLINE1, 1, xmin, xmax) + + # Allocate a weights array and fill it. + call smark (sptemp) + call salloc (weights, PDM_NPT(pdmp), TY_DOUBLE) + do i = 1, PDM_NPT(pdmp) + Memd[weights+i-1] = PDM_INUSE(pdmp,i) + + call dcvfit (PDM_CVD(pdmp), PDM_X(pdmp,1), PDM_DY(pdmp,1), + Memd[weights], PDM_NPT(pdmp), WTS_USER, ier) + if (ier != 0) { + call eprintf ("error in dcvfit\n") + call erract (EA_WARN) + } + + # Subtract the fit from the data. + if (PDM_RESID(pdmp) == YES) { + call printf ("Already using residuals. \007\n") + next + } else { + # For each point, calculate the fit function and subtract + # it from the data. + + do i = 1, PDM_NPT(pdmp) { + PDM_DY(pdmp,i) = PDM_DY(pdmp,i) - + dcveval (PDM_CVD(pdmp), PDM_X(pdmp,i)) + } + PDM_RESID(pdmp) = YES + } + + call dcvfree (PDM_CVD(pdmp)) + call sfree (sptemp) + + # Replot. + call pdm_dplot (pdmp, infile, flip) + ptype = DATAPLOT + case 'z': + # Flip the y-axis scale. + flip = !flip + case 'q': + # Quit. Exit PDM. + if (gqverify() == YES) + break + default: + # Error: unknown command: ring bell. + call printf ("\007\n") + call printf ("? for help or (h,f,i,k,p") + call printf (",d,u,j,s,t,g,a,m,r,x,q,:)\n") + } + } + + call flush (STDOUT) + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmdelete.x b/noao/astutil/pdm/pdmdelete.x new file mode 100644 index 00000000..5052a2e7 --- /dev/null +++ b/noao/astutil/pdm/pdmdelete.x @@ -0,0 +1,103 @@ +include +include +include +include +include "pdm.h" + +define MSIZE 2.0 # Mark size + +# PDM_DELETE -- Delete the point from the plot and set it's inuse array +# entry to zero. + +int procedure pdm_delete (pdmp, cx, cy, ptype) + +pointer pdmp # pointer to PDM data structure +double cx, cy # device cursor coordinates +int ptype # plot type + +pointer gp +real x, y +int npts, i, j, index +real x0, y0, r2, r2min + +begin + gp = PDM_GP(pdmp) + npts = PDM_NPT(pdmp) + + if (ptype == DATAPLOT) { + # Transform world cursor coordinates to NDC. + call gctran (gp, real(cx), real(cy), x, y, 1, 0) + + # Search for nearest point in-use. + j = 0 + r2min = MAX_REAL + do i = 1, npts { + if (PDM_INUSE(pdmp,i) == 0) + next + + call gctran (gp, real(PDM_X(pdmp,i)), real(PDM_DY(pdmp,i)), + x0, y0, 1, 0) + + r2 = (x0 - x) ** 2 + (y0 - y) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Mark the deleted point with a cross and set the weight to zero. + if (j != 0) { + call gscur (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j))) + call gmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)), + GM_CROSS, MSIZE, MSIZE) + PDM_INUSE(pdmp,j) = 0 + call gflush(gp) + } + + return (j) + + } else if (ptype == PHASEPLOT) { + # Transform world cursor coordinates to NDC. + call gctran (gp, real(cx), real(cy), x, y, 1, 0) + + # Search for nearest point in-use. + j = 0 + r2min = MAX_REAL + do i = 1, npts { + index = PDM_SORT(pdmp,i) + if (PDM_INUSE(pdmp,index) == 0) + next + + call gctran (gp, real(PDM_XPH(pdmp,i)), real(PDM_YPH(pdmp,i)), + x0, y0, 1, 0) + + r2 = (x0 - x) ** 2 + (y0 - y) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Mark the deleted point with a cross and set the weight to zero. + # We mark two points since all of the points are displayed twice + # on the phase plot. + + if (j != 0) { + call gscur (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j))) + call gmark (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j)), GM_CROSS, MSIZE, MSIZE) + call gscur (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j))) + call gmark (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)), + GM_CROSS, MSIZE, MSIZE) + + # Calculate which point this corresponds to. + index = PDM_SORT(pdmp,j) + PDM_INUSE(pdmp,index) = 0 + call gflush (gp) + } + + return (index) + } else + return (0) +end diff --git a/noao/astutil/pdm/pdmdplot.x b/noao/astutil/pdm/pdmdplot.x new file mode 100644 index 00000000..c0b53d19 --- /dev/null +++ b/noao/astutil/pdm/pdmdplot.x @@ -0,0 +1,101 @@ +include +include +include +include +include "pdm.h" + +define MSIZE 2.0 # Mark size +define EDGESCL 5.0 # Percent extra space around plot data + +# PDM_DPLOT -- Plot the data on the screen. + +procedure pdm_dplot (pdmp, filename, flip) + +pointer pdmp # pointer to PDM data structure +char filename[SZ_FNAME] # name of the input data file +bool flip # flip the y-axis scale + +int npt, i +real x1, x2, y1, y2, scldif, sclspc +pointer gp, xtemp, ytemp +pointer title +pointer system_id, sp + +begin + call smark (sp) + call salloc (system_id, SZ_LINE, TY_CHAR) + call salloc (title, 4*SZ_LINE, TY_CHAR) + call salloc (xtemp, PDM_NPT(pdmp), TY_REAL) + call salloc (ytemp, PDM_NPT(pdmp), TY_REAL) + + npt = PDM_NPT(pdmp) + gp = PDM_GP(pdmp) + call gclear (gp) + + # Scale the wcs. + do i = 1, PDM_NPT(pdmp) { + Memr[xtemp+i-1] = real(PDM_X(pdmp,i)) + Memr[ytemp+i-1] = real(PDM_DY(pdmp,i)) + } + call gascale (gp, Memr[xtemp], npt, 1) + call gascale (gp, Memr[ytemp], npt, 2) + + # Get the X and Y boundaries. + call ggwind (gp, x1, x2, y1, y2) + + # Add boundry space. + scldif = x2 - x1 + sclspc = scldif * (EDGESCL / 100.) + x1 = x1 - sclspc + x2 = x2 + sclspc + scldif = y2 - y1 + sclspc = scldif * (EDGESCL / 100.) + y1 = y1 - sclspc + y2 = y2 + sclspc + + # Flip the y-axis scale if flip = TRUE. + if (flip) + call gswind (gp, x1, x2, y2, y1) + else + call gswind (gp, x1, x2, y1, y2) + + # Multiline title, save in an array and sprintf to it. + # Get the system identification. + + call sysid (Memc[system_id], SZ_LINE) + call sprintf (Memc[title], 4*SZ_LINE, + "%s\nFile = %s\n%s\nnumpts = %d") + call pargstr (Memc[system_id]) + call pargstr (filename) + if (PDM_RESID(pdmp) == YES) + call pargstr ("Data with fit removed") + else + call pargstr ("Data") + call pargi (npt) + + # Draw the axes. + call glabax (gp, Memc[title], "obs time", "magnitude") + + # Make the plot. + if (npt <= PDM_PLUSPOINT(pdmp)) { + call gpmark (gp, Memr[xtemp], Memr[ytemp], npt, + GM_PLUS, MSIZE, MSIZE) + } else { + call gpmark (gp, Memr[xtemp], Memr[ytemp], npt, + GM_POINT, 1.0, 1.0) + } + + # Call the routine to mark the ranges if they are in effect. + call rg_gxmarkd (gp, PDM_SAMPLE(pdmp), PDM_X(pdmp,1), npt, 1) + + # Draw an x over any deleted points. + do i = 1, npt { + if (PDM_INUSE(pdmp,i) == 0) { + call gscur (gp, Memr[xtemp+i-1], Memr[ytemp+i-1]) + call gmark (gp, Memr[xtemp+i-1], Memr[ytemp+i-1], GM_CROSS, + MSIZE, MSIZE) + } + } + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmfindmin.x b/noao/astutil/pdm/pdmfindmin.x new file mode 100644 index 00000000..4053bf43 --- /dev/null +++ b/noao/astutil/pdm/pdmfindmin.x @@ -0,0 +1,57 @@ +include +include +include +include "pdm.h" + +# PDM_FINDMIN -- Find the minimum value of the abscissa. + +int procedure pdm_findmin (pdmp, ptype, startint, endint, is, ie) + +pointer pdmp # pointer to PDM data structure +int ptype # plot type +double startint, endint # start and end ordinates +int is, ie # start and end indexs + +double miny, dy, dx +int i, isave, npt +pointer xpt, ypt + +begin + # Dereference npt. + npt = PDM_NPT(pdmp) + + # Dereference the appropriate abcissa and ordinate. + switch (ptype) { + case DATAPLOT: + xpt = PDM_XP(pdmp) + ypt = PDM_DYP(pdmp) + case THETAPPLOT: + xpt = PDM_XTHP(pdmp) + ypt = PDM_YTHP(pdmp) + case THETAFPLOT: + ypt = PDM_YTHP(pdmp) + xpt = PDM_XTHP(pdmp) + case PHASEPLOT: + ypt = PDM_YPHP(pdmp) + xpt = PDM_XPHP(pdmp) + } + + # Search the abscissas between startint and endint + # for the minimum value. + + isave = 1 + miny = MAX_DOUBLE + do i = is, ie { + dx = Memd[xpt+i-1] + dy = Memd[ypt+i-1] + if (dx > startint && dx < endint) { + if (dy < miny) { + miny = dy + isave = i + } + } + } + + # Return the corresponding index value. + return (isave) +end diff --git a/noao/astutil/pdm/pdmfitphase.x b/noao/astutil/pdm/pdmfitphase.x new file mode 100644 index 00000000..6d9aa901 --- /dev/null +++ b/noao/astutil/pdm/pdmfitphase.x @@ -0,0 +1,43 @@ +include +include +include +include +include "pdm.h" + +# PDM_FITPHASE -- Call ICFIT on the Phase Curve. + +procedure pdm_fitphase (pdmp) + +pointer pdmp # pointer to PDM data structure + +double xmin, xmax +pointer weights +int i, npt +errchk calloc, icg_fit + +begin + # Dereference some pointers. + npt = PDM_NPT(pdmp) + + # Calloc a phase in-use array. + call alimd (PDM_XPH(pdmp,1), npt, xmin, xmax) + call ic_putr (PDM_ICP(pdmp), "xmin", real(xmin)) + call ic_putr (PDM_ICP(pdmp), "xmax", real(xmax)) + call calloc (weights, npt, TY_DOUBLE) + + # Permute the data in-use array and save it in the phase weights array. + do i = 1, npt + Memd[weights+i-1] = double(PDM_INUSE(pdmp,PDM_SORT(pdmp,i))) + + # Call icfit on the phase curve (pass the phase x, y, wts) + if (npt >= 2) + call icg_fitd (PDM_ICP(pdmp), PDM_GP(pdmp), "cursor", + PDM_GT(pdmp), PDM_CVP(pdmp), PDM_XPH(pdmp,1), + PDM_YPH(pdmp,1), Memd[weights], npt) + + # Update the data in-use array appropriately. + do i = 1, npt + PDM_INUSE(pdmp,PDM_SORT(pdmp,i)) = int(Memd[weights+i-1]) + + call mfree (weights, TY_DOUBLE) +end diff --git a/noao/astutil/pdm/pdmgdata.x b/noao/astutil/pdm/pdmgdata.x new file mode 100644 index 00000000..79e8d49f --- /dev/null +++ b/noao/astutil/pdm/pdmgdata.x @@ -0,0 +1,136 @@ +include +include +include +include "pdm.h" + +define SZ_BUF 100 + +# PDM_GDATA -- Get Data from the input files. + +int procedure pdm_gdata (pdmp, infile) + +pointer pdmp # pointer to PDM data structure +char infile[SZ_LINE] # input data file name + +int fntopnb(), list, clgfil() +int n, ncols, lineno, buflen +int open(), getline(), nscan() +int fd +pointer nextfile +pointer lbuf, ip, sp +errchk realloc, fntopnb, open + +begin + # Get a line buffer. + call smark (sp) + call salloc (lbuf, SZ_LINE, TY_CHAR) + call salloc (nextfile, SZ_LINE, TY_CHAR) + + # Open the input file as a list of files. + list = fntopnb (infile, 0) + + # Initialize some variables. + n = 0 + ncols = 0 + lineno = 0 + + # For each input file in the list, read the data. + while (clgfil (list, Memc[nextfile], SZ_FNAME) != EOF) { + + # Open this input file. + fd = open (Memc[nextfile], READ_ONLY, TEXT_FILE) + + # Read in the data from this file. + while (getline (fd, Memc[lbuf]) != EOF) { + # Skip white space and blank lines. + lineno = lineno + 1 + for (ip = lbuf; IS_WHITE(Memc[ip]); ip = ip + 1) + ; + if (Memc[ip] == '\n' || Memc[ip] == EOS) + next + + if (n == 0) { + buflen = SZ_BUF + iferr { + call calloc (PDM_XP(pdmp), buflen, TY_DOUBLE) + call calloc (PDM_DYP(pdmp), buflen, TY_DOUBLE) + call calloc (PDM_ODYP(pdmp), buflen, TY_DOUBLE) + call calloc (PDM_INUSEP(pdmp), buflen, TY_INT) + call calloc (PDM_ERRP(pdmp), buflen, TY_REAL) + } then + call erract (EA_FATAL) + } else if (n + 1 > buflen) { + buflen = buflen + SZ_BUF + call realloc (PDM_XP(pdmp), buflen, TY_DOUBLE) + call realloc (PDM_DYP(pdmp), buflen, TY_DOUBLE) + call realloc (PDM_ODYP(pdmp), buflen, TY_DOUBLE) + call realloc (PDM_INUSEP(pdmp), buflen, TY_INT) + call realloc (PDM_ERRP(pdmp), buflen, TY_REAL) + } + + # Read data from the file, put it in the data structure. + call sscan (Memc[ip]) + call gargd (PDM_X(pdmp,n+1)) + call gargd (PDM_ODY(pdmp,n+1)) + call gargr (PDM_ERR(pdmp,n+1)) + PDM_INUSE(pdmp,n+1) = 1 + PDM_DY(pdmp,n+1) = PDM_ODY(pdmp,n+1) + + # If this is line one, then determine the number of columns. + if (ncols == 0 && nscan() > 0) + ncols = nscan() + + # Check this line against the number of columns and do the + # appropriate thing. + + switch (nscan()) { + case 0: + call printf ("no args; %s, line %d: %s\n") + call pargstr (Memc[nextfile]) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + case 1: + if (ncols >= 2) { + call eprintf ("only one arg; %s, line %d: %s\n") + call pargstr (Memc[nextfile]) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + } else { + PDM_ODY(pdmp,n+1) = PDM_X(pdmp,n+1) + PDM_DY(pdmp,n+1) = PDM_X(pdmp,n+1) + PDM_X(pdmp,n+1) = n + 1.0d+0 + PDM_ERR(pdmp,n+1) = 0.0 + } + case 2: + if (ncols == 3) { + call eprintf ("only two args; %s, line %d: %s\n") + call pargstr (Memc[nextfile]) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + } else { + PDM_ODY(pdmp,n+1) = PDM_ODY(pdmp,n+1) + PDM_DY(pdmp,n+1) = PDM_ODY(pdmp,n+1) + PDM_X(pdmp,n+1) = PDM_X(pdmp,n+1) + PDM_ERR(pdmp,n+1) = 0.0 + } + + } + + n = n + 1 + } + call close (fd) + } + + call realloc (PDM_XP(pdmp), n, TY_DOUBLE) + call realloc (PDM_DYP(pdmp), n, TY_DOUBLE) + call realloc (PDM_ODYP(pdmp), n, TY_DOUBLE) + call realloc (PDM_INUSEP(pdmp), n, TY_INT) + call realloc (PDM_ERRP(pdmp), n, TY_REAL) + + call fntclsb (list) + call sfree (sp) + return (n) +end diff --git a/noao/astutil/pdm/pdmminmaxp.x b/noao/astutil/pdm/pdmminmaxp.x new file mode 100644 index 00000000..01cdcd2c --- /dev/null +++ b/noao/astutil/pdm/pdmminmaxp.x @@ -0,0 +1,43 @@ +include +include +include +include "pdm.h" + +# PDM_MINMAXP -- Calculate the minimum and maximum periods automatically. + +procedure pdm_minmaxp (pdmp) + +pointer pdmp # pointer to PDM data structure + +int npt, i +double minp, maxp +pointer sep, sp + +begin + call smark (sp) + npt = PDM_NPT(pdmp) + + # Allocate an array of separations and fill it. Find the minimum + # separation as we go. + + call salloc (sep, npt-1, TY_DOUBLE) + maxp = PDM_X(pdmp,npt) - PDM_X(pdmp,1) + minp = maxp + do i = 1, npt-1 { + Memd[sep+i-1] = PDM_X(pdmp,i+1) - PDM_X(pdmp,i) + if (Memd[sep+i-1] < minp) + minp = Memd[sep+i-1] + } + + # Set minp equal to twice this minimum (Nyquist criterion). Set fmax. + PDM_PMIN(pdmp) = 2.0d+0 * minp + if (PDM_PMIN(pdmp) != 0.0d+0) + PDM_FMAX(pdmp) = 1.0d+0/PDM_PMIN(pdmp) + + # Set maxp equal to 4 times maxp. Set fmin. + PDM_PMAX(pdmp) = 4.0d+0 * maxp + if (PDM_PMAX(pdmp) != 0.0d+0) + PDM_FMIN(pdmp) = 1.0d+0/PDM_PMAX(pdmp) + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmopen.x b/noao/astutil/pdm/pdmopen.x new file mode 100644 index 00000000..aa58ad5f --- /dev/null +++ b/noao/astutil/pdm/pdmopen.x @@ -0,0 +1,47 @@ +include +include "pdm.h" + +# PDM_OPEN -- Open a new PDM data structure. + +int procedure pdm_open (device, batchfile, metafile, interactive) + +char device[SZ_FNAME] # graphics device +char batchfile[SZ_FNAME] # file to store batch information +char metafile[SZ_FNAME] # file to store plots +bool interactive # interactive flag + +pointer pdmp +pointer gt_init(), gopen() +int open() +errchk ic_open, gopen, open, calloc, malloc + +begin + # Allocate a pdm data structure. + call calloc (pdmp, PDM_LENSTRUCT, TY_STRUCT) + + # Set up icfit structure pointers for data and phase curve fits. + call ic_open (PDM_ICD(pdmp)) + call ic_open (PDM_ICP(pdmp)) + + # Set up gtools and the gio structure pointers. + PDM_GT(pdmp) = gt_init () + if (interactive) { + PDM_GP(pdmp) = gopen (device, NEW_FILE, STDGRAPH) + PDM_LFD(pdmp) = open ("STDOUT", APPEND, TEXT_FILE) + } else { + PDM_LFD(pdmp) = open (batchfile, APPEND, TEXT_FILE) + PDM_PFD(pdmp) = open (metafile, APPEND, BINARY_FILE) + PDM_GP(pdmp) = gopen ("stdvdm", NEW_FILE, PDM_PFD(pdmp)) + } + + # Allocate space for the sample string and put a '*' in it. + call malloc (PDM_SAMPLEP(pdmp), SZ_LINE, TY_CHAR) + call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE) + + # Booleans. + PDM_RESID(pdmp) = NO + PDM_RANGE(pdmp) = YES + PDM_EB(pdmp) = NO + + return (pdmp) +end diff --git a/noao/astutil/pdm/pdmphase.x b/noao/astutil/pdm/pdmphase.x new file mode 100644 index 00000000..76b3e0d2 --- /dev/null +++ b/noao/astutil/pdm/pdmphase.x @@ -0,0 +1,72 @@ +include +include +include +include +include "pdm.h" + +# PDM_PHASE -- Calculate Phase Curve for period (compressed light curve). + +procedure pdm_phase (pdmp, period, epoch) + +pointer pdmp # pointer to PDM data structure +double period # period to calculate the phase for +double epoch # epoch of this data + +int j, offset, temp +double p +pointer npt, phaseint, sp +errchk calloc, realloc + +begin + call smark (sp) + npt = PDM_NPT(pdmp) + call salloc (phaseint, npt, TY_DOUBLE) + + + # Allocate space for the output phase data (ordinate and abscissa) + # in the pdm data structure. + + if (PDM_XPHP(pdmp) == NULL) { + call calloc (PDM_XPHP(pdmp), 2*npt, TY_DOUBLE) + call calloc (PDM_YPHP(pdmp), 2*npt, TY_DOUBLE) + call calloc (PDM_PHERRP(pdmp), 2*npt, TY_REAL) + } else { + call realloc (PDM_XPHP(pdmp), 2*npt, TY_DOUBLE) + call realloc (PDM_YPHP(pdmp), 2*npt, TY_DOUBLE) + call realloc (PDM_PHERRP(pdmp), 2*npt, TY_REAL) + } + + # Set up the sort array and a temporary array for the phases. + if (PDM_SORTP(pdmp) == NULL) + call calloc (PDM_SORTP(pdmp), npt, TY_INT) + else + call realloc (PDM_SORTP(pdmp), npt, TY_INT) + + # Calculate the phases for all the points. + for (j=1; j<=npt; j=j+1) { + PDM_SORT(pdmp,j) = j + if (period > EPSILONR) { + temp = (int(epoch/period)+1)*period + p = (PDM_X(pdmp,j) - epoch + temp)/period + } + Memd[phaseint+j-1] = double(p - int(p)) + } + + # Sort the phase array into ascending order and permute + # the index array (sort). + + call pdm_sort (phaseint, PDM_SORTP(pdmp), npt) + + # Store the data in the pdm data structure. + do j = 1, npt { + offset = PDM_SORT(pdmp,j) + PDM_YPH(pdmp,j) = PDM_DY(pdmp,offset) + PDM_XPH(pdmp,j) = Memd[phaseint+offset-1] + PDM_PHERR(pdmp,j) = PDM_ERR(pdmp,offset) + PDM_YPH(pdmp,j+npt) = PDM_DY(pdmp,offset) + PDM_XPH(pdmp,j+npt) = Memd[phaseint+offset-1] + 1.0 + PDM_PHERR(pdmp,j+npt) = PDM_ERR(pdmp,offset) + } + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmpplot.x b/noao/astutil/pdm/pdmpplot.x new file mode 100644 index 00000000..78d0ee04 --- /dev/null +++ b/noao/astutil/pdm/pdmpplot.x @@ -0,0 +1,120 @@ +include +include +include +include +include "pdm.h" + +define MSIZE 2.0 # Mark size +define EDGESCL 5.0 # Percent extra space around plot data + +# PDM_PPLOT -- Plot the Phase curve. + +procedure pdm_pplot (pdmp, period, filename, flip) + +pointer pdmp # pointer to PDM data structure +double period # period for which to plot the phase curve +char filename[SZ_FNAME] # name of the input data file +bool flip # flip the y-axis scale + +int npt, i, index +double frequency +real x1, x2, y1, y2, scldif, sclspc +pointer gp, xtemp, ytemp, etemp +pointer title, system_id, sp + +begin + npt = PDM_NPT(pdmp) + gp = PDM_GP(pdmp) + call gclear (PDM_GP(pdmp)) + + call smark (sp) + call salloc (system_id, SZ_LINE, TY_CHAR) + call salloc (title, PDM_SZ_TITLE, TY_CHAR) + call salloc (xtemp, 2*npt, TY_REAL) + call salloc (ytemp, 2*npt, TY_REAL) + call salloc (etemp, 2*npt, TY_REAL) + + do i = 1, 2*npt { + Memr[xtemp+i-1] = PDM_XPH(pdmp,i) + Memr[ytemp+i-1] = PDM_YPH(pdmp,i) + Memr[etemp+i-1] = PDM_PHERR(pdmp,i) + } + + # Scale the wcs. + call gascale (gp, Memr[xtemp], 2*npt, 1) + call gascale (gp, Memr[ytemp], 2*npt, 2) + + # Get the boundaries in X and Y. + call ggwind (gp, x1, x2, y1, y2) + + # Add boundry space. + scldif = x2 - x1 + sclspc = scldif * (EDGESCL / 100.) + x1 = x1 - sclspc + x2 = x2 + sclspc + scldif = y2 - y1 + sclspc = scldif * (EDGESCL / 100.) + y1 = y1 - sclspc + y2 = y2 + sclspc + + # Flip the y-axis scale if flip = TRUE. + if (flip) + call gswind (gp, x1, x2, y2, y1) + else + call gswind (gp, x1, x2, y1, y2) + + # Multiline title, save in an array and sprintf to it. + # Get the system identification. + + call sysid (Memc[system_id], SZ_LINE) + + # Calculate the frequency. + if (period != 0.0) + frequency = 1.0d+0/period + + call sprintf (Memc[title], PDM_SZ_TITLE, +"%s\nPhase curve at period %12.12g, frequency %12.12g\nFile = %s, numpts = %d") + call pargstr (Memc[system_id]) + call pargd (period) + call pargd (frequency) + call pargstr (filename) + call pargi (npt) + + # Draw the axes. + call glabax (gp, Memc[title], "phase", "magnitude") + + # Make the plot. If error bars are turned on, draw them. + if (PDM_EB(pdmp) == YES && Memr[etemp] > EPSILONR) { + call gpmark (gp, Memr[xtemp], Memr[ytemp], + 2*npt, GM_CIRCLE, 1.0, 1.0) + call gpmark (gp, Memr[xtemp], Memr[ytemp], + 2*npt, GM_CIRCLE+GM_FILL, 1.0, 1.0) + do i = 1, 2*npt + call gpmark (gp, Memr[xtemp+i-1], Memr[ytemp+i-1], + 1, GM_VEBAR, MSIZE, -(2.0*Memr[etemp+i-1])) + } else { + + if (npt <= PDM_PLUSPOINT(pdmp)) { + call gpmark (gp, Memr[xtemp], Memr[ytemp], + 2*npt, GM_PLUS, MSIZE, MSIZE) + } else { + call gpmark (gp, Memr[xtemp], Memr[ytemp], + 2*npt, GM_POINT, 1.0, 1.0) + } + } + + # Draw an x over any deleted points. + do i = 1, npt { + index = PDM_SORT(pdmp,i) + if (PDM_INUSE(pdmp,index) == 0) { + call gscur (gp, Memr[xtemp+i-1]+1, Memr[ytemp+i-1]) + call gmark (gp, Memr[xtemp+i-1]+1, Memr[ytemp+i-1], GM_CROSS, + MSIZE, MSIZE) + call gscur (gp, Memr[xtemp+i-1], Memr[ytemp+i-1]) + call gmark (gp, Memr[xtemp+i-1], Memr[ytemp+i-1], GM_CROSS, + MSIZE, MSIZE) + } + } + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmranperm.x b/noao/astutil/pdm/pdmranperm.x new file mode 100644 index 00000000..1d80264e --- /dev/null +++ b/noao/astutil/pdm/pdmranperm.x @@ -0,0 +1,56 @@ +include +include +include +include +include +include "pdm.h" + +# PDM_RANPERM -- Make a random permutation of the data vector. +# This is the algorithm: +# Starting at the beginning of the input array, +# Do this numdatapoints times { +# Find the next random position. +# Call urand for a random number between one and numdatapoints(hash). +# Move forward in the input array this number of places (mod numpts). +# While (marker array for this position has a zero) { +# Move to the next position (linear rehash). +# } +# Put the input array value found at this position into the +# output array, set the marker array value corresponding +# to this position to zero. +# } + + +procedure pdm_ranperm (inarray, inuse, outarray, outinuse, numpts, seed) + +int numpts # number of points in the data +double inarray[numpts] # data to be permuted +double outarray[numpts] # output permuted data to this array +int inuse[numpts] # the PDM in-use array +int outinuse[numpts] # output scrambled in-use array +long seed # a seed for the random number generator + +int count, pos +real urand() +pointer p, sp # array to keep track of which have been used + +begin + # Allocate the output array and the marker array. + # Fill the marker array with ones (amovki) + + call smark (sp) + call salloc (p, numpts, TY_INT) + call amovki (1, Memi[p], numpts) # Set this array to all ones. + + pos = 0 + do count = 1, numpts { + pos = mod(pos+int(urand(seed)*numpts)+1, numpts) # Hash. + while (Memi[p+pos] == 0) # Linear rehash. + pos = mod(pos+1, numpts) + outarray[count] = inarray[pos+1] + outinuse[count] = inuse[pos+1] + Memi[p+pos] = 0 + } + + call sfree (sp) +end diff --git a/noao/astutil/pdm/pdmshow.x b/noao/astutil/pdm/pdmshow.x new file mode 100644 index 00000000..6f8af37a --- /dev/null +++ b/noao/astutil/pdm/pdmshow.x @@ -0,0 +1,56 @@ +include +include +include +include "pdm.h" + +# PDM_SHOW -- Print information to file. + +procedure pdm_show (pdmp, file, verbose) + +pointer pdmp # pointer to PDM data structure +char file[ARB] # file to put the show information in +bool verbose # verbose output flag + +int fd, open(), i +errchk open() + +begin + # Open the output file. + fd = open (file, APPEND, TEXT_FILE) + + # Print information from the data structure. + call fprintf (fd, + "minimum period searched = %12.12g, maximum = %g12.12\n") + call pargd (PDM_PMIN(pdmp)) + call pargd (PDM_PMAX(pdmp)) + call fprintf (fd, + "period = %12.12g, amplitude = %12.12g, epoch = %12.12g\n") + call pargd (PDM_MINR(pdmp)) + call pargd (PDM_AMPL(pdmp)) + call pargd (PDM_EPOCH(pdmp)) + + if (verbose) { + # Print the working data set out as x,y,in-use triplets. + call fprintf (fd, "The working data vector is as follows: \n") + do i = 1, PDM_NPT(pdmp) { + call fprintf ( fd, "index = %d, x = %12.12g, y = %12.12g\n") + call pargi (i) + call pargd (PDM_X(pdmp,i)) + call pargd (PDM_DY(pdmp,i)) + } + + if (PDM_XPHP(pdmp) != NULL) { + # Print the phasecurve out as x,y pairs + call fprintf (fd, "\nThe phase curve vector is as follows: \n") + do i = 1, PDM_NPT(pdmp) { + call fprintf ( fd, "index = %d, x = %12.12g, y = %12.12g\n") + call pargi (i) + call pargd (PDM_XPH(pdmp,i)) + call pargd (PDM_YPH(pdmp,i)) + } + } + } + + # Close the output file. + call close (fd) +end diff --git a/noao/astutil/pdm/pdmsignif.x b/noao/astutil/pdm/pdmsignif.x new file mode 100644 index 00000000..13a17acc --- /dev/null +++ b/noao/astutil/pdm/pdmsignif.x @@ -0,0 +1,61 @@ +include +include +include +include +include +include "pdm.h" + +define NUMTRIES 100 + +# PDM_SIGNIF -- Calculate the significance of the theta statistic for +# a certain period. Use the "Method of Randomization". + +double procedure pdm_signif (pdmp, period) + +pointer pdmp # pointer to PDM data structure +double period # period at which to find significance + +int lesscount, i, npt +double otheta, theta, pdm_theta(), pdm_thetaran() +long seed +pointer rg, pt, inuse, oinuse, rg_xrangesd(), sp +errchk pdm_statistics, pdm_theta(), pdm_ranperm, pdm_thetaran() + +begin + # Do NUMTRIES random permutations on the data, calculate the theta + # statistic on the scrambled data for this period. Return the + # fraction of these permutaions which yield thetas less than + # the unmixed data. + + call smark (sp) + npt = PDM_NPT(pdmp) + + # Make sure the statistics are up to date. + call pdm_statistics (pdmp) + + # Allocate a temporary array for the scrambled data, allocate a + # temporary copy of the inuse array and copy the real inuse array + # into it. + + call salloc (pt, npt, TY_DOUBLE) + call salloc (inuse, npt, TY_INT) + call salloc (oinuse, npt, TY_INT) + call amovi (PDM_INUSE(pdmp,1), Memi[inuse], npt) + lesscount = 0 + + # Calculate the ranges information from the sample string. + rg = rg_xrangesd (PDM_SAMPLE(pdmp), PDM_X(pdmp,1), npt) + otheta = pdm_theta (pdmp, rg, period) + seed = 1.0 + + do i = 1, NUMTRIES { + call pdm_ranperm (PDM_DY(pdmp,1), Memi[inuse], Memd[pt], + Memi[oinuse], npt, seed) + theta = pdm_thetaran (pdmp, pt, oinuse, rg, period) + if (theta < otheta) + lesscount = lesscount + 1 + } + + call sfree (sp) + return (1.0d+0 - (double(lesscount)/double(NUMTRIES))) +end diff --git a/noao/astutil/pdm/pdmsort.x b/noao/astutil/pdm/pdmsort.x new file mode 100644 index 00000000..93c59aec --- /dev/null +++ b/noao/astutil/pdm/pdmsort.x @@ -0,0 +1,20 @@ +include +include "pdm.h" + +# PDM_SORT -- Sort the phases into ascending order. + +procedure pdm_sort (array, sort, numpts) + +pointer array # array to sort +pointer sort # array to contain the sort indexes +int numpts # number of points to sort + +pointer comarray +common /sortcom/ comarray +extern pdm_compare() + +begin + comarray = array + if (numpts > 1) + call qsort (Memi[sort], numpts, pdm_compare) +end diff --git a/noao/astutil/pdm/pdmstats.x b/noao/astutil/pdm/pdmstats.x new file mode 100644 index 00000000..aeabacaa --- /dev/null +++ b/noao/astutil/pdm/pdmstats.x @@ -0,0 +1,37 @@ +include +include +include +include "pdm.h" + +# PDM_STATISTICS -- Calculate the sum of squares and the variance of the data. + +procedure pdm_statistics (pdmp) + +pointer pdmp # pointer to PDM data structure + +int npt, i, j +double var, sumx2, sum + +begin + npt = PDM_NPT(pdmp) + + # Calculate the sum of squares and the variance of the data. + sumx2 = 0.0 + sum = 0.0 + j = 0 + + do i = 1, npt { + if (PDM_INUSE(pdmp,i) == 1) { + sumx2 = sumx2 + PDM_DY(pdmp,i)**2 # Sum of squares. + sum = sum + PDM_DY(pdmp,i) + j = j + 1 + } + } + + if (j != 1) + var = (sumx2 - sum**2/double(j))/double(j - 1) # Variance. + + # Put these two values in the data structure. + PDM_SUMSQ(pdmp) = sumx2 + PDM_DVAR(pdmp) = var +end diff --git a/noao/astutil/pdm/pdmtheta.x b/noao/astutil/pdm/pdmtheta.x new file mode 100644 index 00000000..b2b295bc --- /dev/null +++ b/noao/astutil/pdm/pdmtheta.x @@ -0,0 +1,120 @@ +include +include +include +include +include "pdm.h" + +# PDM_THETA -- Calculate the theta statistic for this period. +# Theta is a measure of the dispersion of data phases about a mean +# light curve. + +double procedure pdm_theta (pdmp, rg, period) + +pointer pdmp # pointer to PDM data structure +pointer rg # segments structure pointer +double period # period at which to find theta + +int i, j, k, l, m +double s2 +int ndof, segst, segend +bool bin10 +double theta +pointer sumbin, sum2bin, numbin, sp +errchk binem + +begin + # Allocate bin storage. + call smark (sp) + call salloc (sumbin, 10, TY_DOUBLE) + call salloc (sum2bin, 10, TY_DOUBLE) + call salloc (numbin, 10, TY_INT) + + ndof = 0 + s2 = 0 + + # Do loop on the segments. + do i = 1, RG_NRGS(rg) { + + # Calculate segst, segend, bin10. + segst = min(RG_X2(rg,i),RG_X1(rg,i)) + segend = max(RG_X2(rg,i),RG_X1(rg,i)) + bin10 = ((segend - segst) >= BIN10) + + # Calculate the number of points in each bin and the sum of + # the bins. + + call binem (period, bin10, PDM_XP(pdmp), PDM_DYP(pdmp), + segst, segend, PDM_INUSEP(pdmp), sumbin, sum2bin, numbin) + + # Calculate sigma**2 for this period. + for (j=0; j<=9; j=j+1) { + k = numbin+j + l = sumbin+j + m = sum2bin+j + if (Memi[k] > 1) { + ndof = ndof + Memi[k] - 1 + s2 = s2 + (Memd[m] - Memd[l]**2 / Memi[k]) + } + } + } + + # Calculate theta. + theta = (s2 / double(ndof)) / PDM_DVAR(pdmp) + + call sfree (sp) + + return (theta) + +end + + +# BINEM -- Put the data points into the appropriate bins. + +procedure binem(incper, bin10, x, y, segst, segend, inuse, sumbin, sum2bin, + numbin) + +double incper # period increment +bool bin10 # use 5 or 10 bins flag +pointer x # ordinates +pointer y # abcissas +pointer inuse # PDM in-use array +int segst, segend # segment start and segment end +pointer sumbin, sum2bin, numbin # pointers to bins of sum, sum2, and number + +int bin1, bin2, j, k, l, m +double p, phase, p0 + +begin + do j = 1, 10 { + Memi[numbin+j-1] = 0 + Memd[sumbin+j-1] = 0.0 + Memd[sum2bin+j-1] = 0.0 + } + + #p0 = Memd[x] + call alimd (Memd[x+segst-1], segend-segst+1, p0, p) + do j = segst, segend { + if (Memi[inuse+j-1] == 0) + next + p = (Memd[x+j-1] - p0)/incper + phase = double(p - int(p)) + if (bin10) { + bin1 = mod(int(10.*phase+0.5), 10) + } else { + bin1 = 2 * int(5. * phase) + 1 + bin2 = 2 * (mod(int(5. * phase + 0.5), 5)) + k = numbin+bin2 + l = sumbin+bin2 + m = sum2bin+bin2 + Memi[k] = Memi[k] + 1 + Memd[l] = Memd[l] + Memd[y+j-1] + Memd[m] = Memd[m] + Memd[y+j-1] ** 2 + } + k = numbin+bin1 + l = sumbin+bin1 + m = sum2bin+bin1 + Memi[k] = Memi[k] + 1 + Memd[l] = Memd[l] + Memd[y+j-1] + Memd[m] = Memd[m] + Memd[y+j-1] ** 2 + } +end diff --git a/noao/astutil/pdm/pdmthetaran.x b/noao/astutil/pdm/pdmthetaran.x new file mode 100644 index 00000000..28e63b34 --- /dev/null +++ b/noao/astutil/pdm/pdmthetaran.x @@ -0,0 +1,118 @@ +include +include +include +include +include "pdm.h" + +# PDMTHETARAN -- This program is a copy of pdmtheta but can be used on +# scrambled data. + +double procedure pdm_thetaran (pdmp, y, inuse, rg, period) + +pointer pdmp # pointer to PDM data structure +pointer y # pointer to abcissas +pointer inuse # pointer to PDM in-use array +pointer rg # pointer to ranges structure +double period # period to calculate theta for + +int i, j, k, l, m +double s2 +int ndof, segst, segend +bool bin10 +double theta +pointer sumbin, sum2bin, numbin, sp +errchk binemran + +begin + # Allocate bin storage. + call smark (sp) + call salloc (sumbin, 10, TY_DOUBLE) + call salloc (sum2bin, 10, TY_DOUBLE) + call salloc (numbin, 10, TY_INT) + + s2 = 0 + ndof = 0 + + # Do loop on the segments. + do i = 1, RG_NRGS(rg) { + + # Calculate segst, segend, bin10. + segst = min(RG_X2(rg,i),RG_X1(rg,i)) + segend = max(RG_X2(rg,i),RG_X1(rg,i)) + bin10 = ((segend - segst) >= BIN10) + + # Calculate the number of points in each bin and the sum of + # the bins. + + call binemran (period, bin10, PDM_XP(pdmp), y, segst, + segend, inuse, sumbin, sum2bin, numbin) + + # Calculate sigma**2 for this period. + for (j=0; j<=9; j=j+1) { + k = numbin+j + l = sumbin+j + m = sum2bin+j + if (Memi[k] > 1) { + ndof = ndof + Memi[k] - 1 + s2 = s2 + (Memd[m] - Memd[l]**2 / Memi[k]) + } + } + } + + # Calculate theta. + theta = (s2 / double(ndof)) / PDM_DVAR(pdmp) + + call sfree (sp) + return (theta) +end + + +# BINEMRAN -- Put the data points into the appropriate bins (scrambled data). + +procedure binemran (incper, bin10, x, y, segst, segend, inuse, sumbin, sum2bin, + numbin) + +double incper +bool bin10 +pointer x +pointer y +pointer inuse +int segst, segend +pointer sumbin, sum2bin, numbin + +int bin1, bin2, j, k, l, m +double p, phase, p0 + +begin + do j = 1, 10 { + Memi[numbin+j-1] = 0 + Memd[sumbin+j-1] = 0.0 + Memd[sum2bin+j-1] = 0.0 + } + + p0 = Memd[x] + do j = segst, segend { + if (Memi[inuse+j-1] == 0) + next + p = (Memd[x+j-1] - p0)/incper + phase = double(p - int(p)) + if (bin10) { + bin1 = mod(int(10.*phase+0.5d+0), 10) + } else { + bin1 = 2 * int(5.0d+0 * phase) + 1 + bin2 = 2 * (mod(int(5.0d+0 * phase + 0.5d+0), 5)) + k = numbin+bin2 + l = sumbin+bin2 + m = sum2bin+bin2 + Memi[k] = Memi[k] + 1 + Memd[l] = Memd[l] + Memd[y+j-1] + Memd[m] = Memd[m] + Memd[y+j-1] ** 2 + } + k = numbin+bin1 + l = sumbin+bin1 + m = sum2bin+bin1 + Memi[k] = Memi[k] + 1 + Memd[l] = Memd[l] + Memd[y+j-1] + Memd[m] = Memd[m] + Memd[y+j-1] ** 2 + } +end diff --git a/noao/astutil/pdm/pdmtplot.x b/noao/astutil/pdm/pdmtplot.x new file mode 100644 index 00000000..e8c3ce8a --- /dev/null +++ b/noao/astutil/pdm/pdmtplot.x @@ -0,0 +1,139 @@ +include +include +include +include +include "pdm.h" + +define EDGESCL 5.0 # Percent extra space around plot data + +# PDM_TPLOT -- Plot the data on the screen. + +procedure pdm_tplot (pdmp, porf, filename) + +pointer pdmp # pointer to PDM data structure +int porf # period or frequency flag +char filename[SZ_FNAME] # name of the input data file + +int nthpt +real x1, x2, y1, y2, scldif, sclspc +pointer gp, xtemp, ytemp +pointer title +char system_id[SZ_LINE] +int indx, pdm_findmin() +errchk malloc, pdm_findmin() + +begin + # Dereference some structure stuff. + nthpt = PDM_NTHPT(pdmp) + gp = PDM_GP(pdmp) + + call malloc (title, PDM_SZ_TITLE, TY_CHAR) + call malloc (xtemp, nthpt, TY_REAL) + call malloc (ytemp, nthpt, TY_REAL) + call gclear (gp) + + do indx = 1, nthpt { + Memr[xtemp+indx-1] = real(PDM_XTH(pdmp,indx)) + Memr[ytemp+indx-1] = real(PDM_YTH(pdmp,indx)) + } + + if (porf == THETAPPLOT) { + # Scale the wcs. + call gascale (gp, Memr[xtemp], nthpt, 1) + call gascale (gp, Memr[ytemp], nthpt, 2) + + # Get the X and Y boundary values. + call ggwind (gp, x1, x2, y1, y2) + + # Add boundry space. + scldif = x2 - x1 + sclspc = scldif * (EDGESCL / 100.) + x1 = x1 - sclspc + x2 = x2 + sclspc + scldif = y2 - y1 + sclspc = scldif * (EDGESCL / 100.) + y1 = y1 - sclspc + y2 = y2 + sclspc + + call gswind (gp, x1, x2, y1, y2) + + # Find the minimum value and save it in the remembered minimum. + indx = pdm_findmin (pdmp, THETAPPLOT, PDM_PMIN(pdmp), + PDM_PMAX(pdmp), 1, nthpt) + PDM_MINR(pdmp) = PDM_XTH(pdmp,indx) + + # Multiline title, save in an array and sprintf to it. + # Get the system identification. + + call sysid (system_id, SZ_LINE) + call sprintf (Memc[title], PDM_SZ_TITLE, + "%s\nFile = %s, minimum = %12.12g\n%s\nnumpts = %d") + call pargstr (system_id) + call pargstr (filename) + call pargd (PDM_MINR(pdmp)) + call pargstr ("Theta vs Period") + call pargi (nthpt) + + # Draw the axes. + call glabax (gp, Memc[title], "period", "theta") + + # Make the plot. + call gpline (gp, Memr[xtemp], Memr[ytemp], nthpt) + + # Put the cursor at the minimum. + call gscur (gp, Memr[xtemp+indx-1], Memr[ytemp+indx-1]) + } else { + # Scale the wcs. + call gascale (gp, Memr[xtemp], nthpt, 1) + call gascale (gp, Memr[ytemp], nthpt, 2) + + # Get the X and Y boundary values. + call ggwind (gp, x1, x2, y1, y2) + + # Add boundry space. + scldif = x2 - x1 + sclspc = scldif * (EDGESCL / 100.) + x1 = x1 - sclspc + x2 = x2 + sclspc + scldif = y2 - y1 + sclspc = scldif * (EDGESCL / 100.) + y1 = y1 - sclspc + y2 = y2 + sclspc + + call gswind (gp, x1, x2, y1, y2) + + # Find the minimum value and save it in the remembered minimum. + indx = pdm_findmin (pdmp, THETAFPLOT, PDM_FMIN(pdmp), + PDM_FMAX(pdmp), 1, nthpt) + if (PDM_XTH(pdmp,indx) > EPSILOND) + PDM_MINR(pdmp) = 1./PDM_XTH(pdmp,indx) + + # Multiline title, save in an array and sprintf to it. + # Get the system identification. + + call sysid (system_id, SZ_LINE) + call sprintf (Memc[title], PDM_SZ_TITLE, + "%s\nFile = %s, minimum = %12.12g\n%s\nnumpts = %d") + call pargstr (system_id) + call pargstr (filename) + if (PDM_MINR(pdmp) > EPSILOND) + call pargd (1.0d+0/PDM_MINR(pdmp)) + else + call pargd (0.0d+0) + call pargstr ("Theta vs Frequency") + call pargi (nthpt) + + # Draw the axes. + call glabax (gp, Memc[title], "frequency", "theta") + + # Make the plot. + call gpline (gp, Memr[xtemp], Memr[ytemp], nthpt) + + # Put the cursor at the minimum. + call gscur (gp, Memr[xtemp+indx-1], Memr[ytemp+indx-1]) + } + + call mfree (title, TY_CHAR) + call mfree (xtemp, TY_REAL) + call mfree (ytemp, TY_REAL) +end diff --git a/noao/astutil/pdm/pdmundelete.x b/noao/astutil/pdm/pdmundelete.x new file mode 100644 index 00000000..2f5f0e31 --- /dev/null +++ b/noao/astutil/pdm/pdmundelete.x @@ -0,0 +1,124 @@ +include +include +include +include +include "pdm.h" + +define MSIZE 2.0 # Mark size + +# PDM_UNDELETE -- Undelete the nearest deleted point from the plot +# and set it's inuse array entry to one (in-use). + +int procedure pdm_undelete (pdmp, cx, cy, ptype) + +pointer pdmp # pointer to PDM data structure +double cx, cy # device coordinates of point to undelete +int ptype # plot type flag + +pointer gp +real x, y +int npt, i, j, index +real x0, y0, r2, r2min + +begin + gp = PDM_GP(pdmp) + npt = PDM_NPT(pdmp) + + if (ptype == DATAPLOT) { + # Transform world cursor coordinates to NDC. + call gctran (gp, real(cx), real(cy), x, y, 1, 0) + + # Search for nearest point not in-use. + j = 0 + r2min = MAX_REAL + do i = 1, npt { + if (PDM_INUSE(pdmp,i) == 1) + next + + call gctran (gp, real(PDM_X(pdmp,i)), real(PDM_DY(pdmp,i)), + x0, y0, 1, 0) + + r2 = (x0 - x) ** 2 + (y0 - y) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Unmark the deleted point and reset the weight. + if (j != 0) { + call gscur (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j))) + call gseti (gp, G_PMLTYPE, GL_CLEAR) + call gmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)), + GM_CROSS, MSIZE, MSIZE) + call gseti (gp, G_PMLTYPE, GL_SOLID) + if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp)) + call gpmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)), + 1, GM_PLUS, MSIZE, MSIZE) + else + call gpmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)), + 1, GM_POINT, 1.0, 1.0) + PDM_INUSE(pdmp,j) = 1 + call gflush (gp) + } + + return (j) + } else if (ptype == PHASEPLOT) { + # Transform world cursor coordinates to NDC. + call gctran (gp, real(cx), real(cy), x, y, 1, 0) + + # Search for nearest point not in-use. + j = 0 + r2min = MAX_REAL + do i = 1, npt { + index = PDM_SORT(pdmp,i) + if (PDM_INUSE(pdmp,index) == 1) + next + + call gctran (gp, real(PDM_XPH(pdmp,i)), real(PDM_YPH(pdmp,i)), + x0, y0, 1, 0) + + r2 = (x0 - x) ** 2 + (y0 - y) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Unmark the deleted point and reset the weight. + if (j != 0) { + call gscur (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j))) + call gseti (gp, G_PMLTYPE, GL_CLEAR) + call gmark (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j)), GM_CROSS, MSIZE, MSIZE) + call gseti (gp, G_PMLTYPE, GL_SOLID) + if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp)) + call gpmark (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j)), 1, GM_PLUS, MSIZE, MSIZE) + else + call gpmark (gp, real(PDM_XPH(pdmp,j))+1.0, + real(PDM_YPH(pdmp,j)), 1, GM_POINT, 1.0, 1.0) + + call gscur (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j))) + call gseti (gp, G_PMLTYPE, GL_CLEAR) + call gmark (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)), + GM_CROSS, MSIZE, MSIZE) + call gseti (gp, G_PMLTYPE, GL_SOLID) + if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp)) + call gpmark (gp, real(PDM_XPH(pdmp,j)), + real(PDM_YPH(pdmp,j)), 1, GM_PLUS, MSIZE, MSIZE) + else + call gpmark (gp, real(PDM_XPH(pdmp,j)), + real(PDM_YPH(pdmp,j)), 1, GM_POINT, 1.0, 1.0) + + # Calculate which point this corresponds to. + index = PDM_SORT(pdmp,j) + PDM_INUSE(pdmp,index) = 1 + call gflush (gp) + } + + return (index) + } else + return (0) +end diff --git a/noao/astutil/pdm/t_pdm.x b/noao/astutil/pdm/t_pdm.x new file mode 100644 index 00000000..7498804a --- /dev/null +++ b/noao/astutil/pdm/t_pdm.x @@ -0,0 +1,72 @@ + +include +include +include +include "pdm.h" + +# PDM -- Phase Dispersion Minimization. Find periodicities in light +# curve data. Root program. + +procedure t_pdm() + +char infile[SZ_FNAME] # input data file +char metafile[SZ_FNAME] # batch metafile +char batchfile[SZ_FNAME] # batch log file +char device[SZ_FNAME] # plot device +bool interactive # interactive or batch flag +bool autoranges # autoranges flag + +int ptype # plot type +pointer pdmp, pdm_open() # structure stuff +bool flip, clgetb() +real clgetr() +int clgeti(), pdm_gdata(), pdm_autorang() +errchk pdm_open, pdm_autorang, pdm_batch, pdm_cursor + +begin + # Get some of the CL parameters and open the pdm data structure. + call clgstr ("infiles", infile, SZ_FNAME) + call clgstr ("metafile", metafile, SZ_FNAME) + call clgstr ("batchfile", batchfile, SZ_FNAME) + call clgstr ("device", device, SZ_FNAME) + interactive = clgetb ("interactive") + flip = clgetb ("flip") + pdmp = pdm_open (device, batchfile, metafile, interactive) + + # Clio to get more parameters. + PDM_PMIN(pdmp) = double(clgetr ("minp")) + call clputr ("minp", real(PDM_PMIN(pdmp))) + if (PDM_PMIN(pdmp) > EPSILOND) + PDM_FMAX(pdmp) = 1.0d+0/PDM_PMIN(pdmp) + else + PDM_FMAX(pdmp) = 0.0d+0 + PDM_PMAX(pdmp) = double(clgetr ("maxp")) + call clputr ("maxp", real(PDM_PMAX(pdmp))) + if (PDM_PMAX(pdmp) > EPSILOND) + PDM_FMIN(pdmp) = 1.0d+0/PDM_PMAX(pdmp) + else + PDM_FMIN(pdmp) = 0.0d+0 + PDM_NTHPT(pdmp) = clgeti ("ntheta") + autoranges = clgetb ("autoranges") + PDM_NSIGMA(pdmp) = double(clgetr ("nsigma")) + PDM_PLUSPOINT(pdmp) = clgeti ("pluspoint") + + # Read in the data. + PDM_NPT(pdmp) = pdm_gdata (pdmp, infile) + + # If the autoranges flag is set, call the autorange subroutine. + if (autoranges) + PDM_NRANGE(pdmp) = pdm_autorang (pdmp) + + if (!interactive) + call pdm_batch (pdmp, batchfile, infile, flip) + else { + # Plot the data on the screen and call the cursor loop. + call pdm_dplot (pdmp, infile, flip) + ptype = DATAPLOT + call pdm_cursor(pdmp, ptype, infile, flip) + } + + # Close the pdm data structure. + call pdm_close (pdmp, interactive) +end -- cgit