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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
define ORDER 4 # The order of the spline
# EXTREMA -- Find the extrema in an array of x and y points.
# The input data points are fitted with a cubic interpolation spline. The
# spline is then searched for points where the first derivative changes sign.
# The minimum step size of this search is controlled by the parameter dx.
# The positions of these extrema are returned in the x array, the value of the
# spline at the extrema are returned in the y array, and the curvature or
# second derivative of the spline at the extrema are returned in the
# curvature array. The function returns the number of extrema found.
int procedure extrema (x, y, curvature, npts, dx)
real x[npts], y[npts] # Input data points and output extrema
real curvature[npts] # 2nd deriv. of cubic spline at extrema
int npts # Number of input data points
real dx # Precision of extrema positions
int i, ier, nextrema
real xeval, left_deriv, right_deriv
pointer sp, bspln, q
real seval()
errchk salloc, seval
begin
# Allocate working arrays for spline routines
call smark (sp)
call salloc (bspln, 2 * npts + 30, TY_REAL)
call salloc (q, (2 * ORDER - 1) * npts, TY_REAL)
# Calculate the spline coefficients
call spline (x, y, npts, Memr[bspln], Memr[q], ORDER, ier)
if (ier != 0) {
call sfree (sp)
return (0)
}
# Initialize the curvature array
call aclrr (curvature, npts)
# Find the extrema defined by a change in sign in the first derivative.
nextrema = 0
left_deriv = seval (x[1], 1, Memr[bspln])
do i = 2, npts {
xeval = x[i]
right_deriv = seval (xeval, 1, Memr[bspln])
if (left_deriv * right_deriv <= 0.) {
for (xeval = x[i - 1] + dx; xeval <= x[i]; xeval = xeval + dx) {
right_deriv = seval (xeval, 1, Memr[bspln])
if (left_deriv * right_deriv <= 0.)
break
left_deriv = right_deriv
}
nextrema = nextrema + 1
x[nextrema] = xeval
y[nextrema] = seval (xeval, 0, Memr[bspln])
curvature[nextrema] = seval (xeval, 2, Memr[bspln])
if (curvature[nextrema] == 0.)
nextrema = nextrema - 1
if (nextrema == npts)
break
}
left_deriv = right_deriv
}
call sfree (sp)
return (nextrema)
end
|
