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
|
include <error.h>
include <mach.h>
include <smw.h>
# Function Mode for STEK
# FUN_DO -- Branch and execute proper function
procedure fun_do (key, sh1, y, n, w0, wpc)
int key
pointer sh1
real y[n]
int n
double w0, wpc
char spec2[SZ_FNAME]
int i, nline, nband, nap, strlen()
real const, clgetr()
pointer im, mw, sh2
bool wave_scl
errchk getimage, shdr_rebin
begin
switch (key) {
case 'a': # Absolute value
do i = 1, n
y[i] = abs (y[i])
case 'd': # Dexp (base 10)
const = log10 (MAX_REAL)
do i = 1, n
if (abs (y[i]) < const)
y[i] = 10.0 ** y[i]
else if (y[i] >= const)
y[i] = MAX_REAL
else
y[i] = 0.0
case 'e': # Exp base e
const = log (MAX_REAL)
do i = 1, n
if (abs (y[i]) < const)
y[i] = exp (y[i])
else if (y[i] >= const)
y[i] = MAX_REAL
else
y[i] = 0.0
case 'i': # Inverse
do i = 1, n
if (y[i] != 0.0)
y[i] = 1.0/y[i]
else
y[i] = 0.0
case 'l': # Log10
do i = 1, n
if (y[i] > 0.0)
y[i] = log10 (y[i])
else
y[i] = -0.5
case 'm': # Multiply by constant
const = clgetr ("constant")
call amulkr (y, const, y, n)
case 'n': # Log base e
do i = 1, n
if (y[i] > 0.0)
y[i] = log (y[i])
else
y[i] = -0.5
case 'p': # Add constant
const = clgetr ("constant")
call aaddkr (y, const, y, n)
case 's': # Square root
do i = 1, n
if (y[i] >= 0.0)
y[i] = sqrt (y[i])
else
y[i] = 0.0
case '+', '-', '*', '/': # Binary operations
call printf ("Second spectrum ")
call clgstr ("spec2", spec2, SZ_FNAME)
if (strlen (spec2) == 0)
return
wave_scl = true
nline = 0
nband = 0
nap = 0
im = NULL
mw = NULL
sh2 = NULL
call getimage (spec2, nline, nband, nap, wave_scl, w0, wpc,
"angstroms", im, mw, sh2, NULL)
call shdr_rebin (sh2, sh1)
switch (key) {
case '+':
call aaddr (y, Memr[SY(sh2)], y, n)
case '-':
call asubr (y, Memr[SY(sh2)], y, n)
case '*':
call amulr (y, Memr[SY(sh2)], y, n)
case '/':
do i = 1, n
if (Memr[SY(sh2)+i-1] == 0.0)
y[i] = 0.0
else
y[i] = y[i] / Memr[SY(sh2)+i-1]
}
call shdr_close (sh2)
call smw_close (mw)
call imunmap (im)
# Redraw
case 'r':
;
default:
call error (0, "Unknown function")
}
end
# FUN_HELP
procedure fun_help ()
begin
call printf ("q=quit l,n=log10,e d,e=d,exp s=sqrt a=abs i=1/s")
call printf (" p=+k m=*k +,-,*,/=2spec ops\n")
end
|