aboutsummaryrefslogtreecommitdiff
path: root/noao/astutil/pdm/pdmcolon.x
blob: 6d438515c871f97d4a13357007a68d522c4e6120 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
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