aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/splot/splotfun.x
blob: 4c94350f273527c0c38948cd63f5fec56ed2397e (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
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