aboutsummaryrefslogtreecommitdiff
path: root/noao/astutil/pdm/pdmcolon.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/astutil/pdm/pdmcolon.x')
-rw-r--r--noao/astutil/pdm/pdmcolon.x292
1 files changed, 292 insertions, 0 deletions
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 <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+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
generated by cgit (git 2.34.1) at 2025-09-20 10:41:35 -0400