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
|
include <math/curfit.h>
# MS_SMOOTH -- Smooth MULTISPEC parameters with the CURFIT package.
# MS_SET_SMOOTH -- Initialize and define function for smoothing.
# MS_FREE_SMOOTH -- Free allocated memory from smoothing.
# This procedure is numerical and does not depend on the MULTISPEC
# package.
procedure ms_smooth (x, y)
real x[ARB] # Array of x values
real y[ARB] # Array of y values
int curve_type # Curfit function
int order # Order of function
real xmin # Minimum x value
real xmax # Maximum x value
int npoints # Number of points in fits
int i, npts, ier
real xmn, xmx
pointer cv, w
real cveval()
data cv/NULL/, w/NULL/
begin
# Check for a valid curfit pointer.
if (cv == NULL)
call error (0, "param_smooth: Undefined smoothing function")
# Zero and fit the data with uniform weights.
call cvzero (cv)
# call cvfit (cv, x, y, Memr[w], npts, WTS_UNIFORM, ier)
# Accumulate points and check for out of bounds points.
do i = 1, npts
if ((x[i] >= xmn) && (x[i] <= xmx))
call cvaccum (cv, x[i], y[i], Memr[w+i-1], WTS_UNIFORM)
call cvsolve (cv, ier)
if (ier != OK)
call error (0, "param_smooth: Error in function fit")
# Evaluate fit placing fit values back in y array.
# call cvvector (cv, x, y, npts)
do i = 1, npts
if ((x[i] >= xmn) && (x[i] <= xmx))
y[i] = cveval (cv, x[i])
return
entry ms_set_smooth (xmin, xmax, npoints)
# Set or reset curfit data structure and allocate memory for weights.
if (cv != NULL)
call cvfree (cv)
if (w == NULL)
call malloc (w, npoints, TY_REAL)
# Determine curve_type and order.
call clgcurfit ("function", "order", curve_type, order)
# Initialize curfit data structure and record number of points.
xmn = xmin
xmx = xmax
call cvinit (cv, curve_type, order, xmn, xmx)
npts = npoints
return
entry ms_free_smooth ()
# Free allocated memory.
if (cv != NULL)
call cvfree (cv)
if (w != NULL)
call mfree (w, TY_REAL)
end
|
