1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ; returns the values of the first derivative of |
---|
5 | ; the interpolating function at the points X2i. It is a double |
---|
6 | ; precision array. |
---|
7 | ; |
---|
8 | ; Given the arrays X and Y, which tabulate a function (with the X[i] |
---|
9 | ; and Y[i] in ascending order), and given an input value X2, the |
---|
10 | ; <pro>spl_incr</pro> function returns an interpolated value for the given |
---|
11 | ; values of X2. The interpolation method is based on cubic spline, corrected |
---|
12 | ; in a way that interpolated value are also in ascending order. |
---|
13 | ; |
---|
14 | ; @examples |
---|
15 | ; |
---|
16 | ; IDL> y2 = spl_fstdrv(x, y, yscd, x2) |
---|
17 | ; |
---|
18 | ; @param x {in}{required} |
---|
19 | ; An n-elements (at least 2) input vector that specifies the |
---|
20 | ; tabulate points in ascending order. |
---|
21 | ; |
---|
22 | ; @param y {in}{required} |
---|
23 | ; f(x) = y. An n-elements input vector that specifies the values |
---|
24 | ; of the tabulated function F(Xi) corresponding to Xi. |
---|
25 | ; |
---|
26 | ; @param yscd {in}{required} |
---|
27 | ; The output from <proidl>SPL_INIT</proidl> for the specified X and Y. |
---|
28 | ; |
---|
29 | ; @param x2 {in}{required} {type= scalar or array} |
---|
30 | ; The input values for which the first derivative values are desired. |
---|
31 | ; |
---|
32 | ; @returns |
---|
33 | ; y2: f'(x2) = y2. |
---|
34 | ; |
---|
35 | ; @history |
---|
36 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005 |
---|
37 | ; |
---|
38 | ; @version |
---|
39 | ; $Id$ |
---|
40 | ; |
---|
41 | ;- |
---|
42 | FUNCTION spl_fstdrv, x, y, yscd, x2 |
---|
43 | ; |
---|
44 | compile_opt idl2, strictarrsubs |
---|
45 | ; |
---|
46 | ; compute the first derivative of the spline function |
---|
47 | ; |
---|
48 | nx = n_elements(x) |
---|
49 | ny = n_elements(y) |
---|
50 | ; x must have at least 2 elements |
---|
51 | IF nx LT 2 THEN stop |
---|
52 | ; y must have the same number of elements than x |
---|
53 | IF nx NE ny THEN stop |
---|
54 | ; define loc in a way that |
---|
55 | ; if loc[i] eq -1 : x2[i] < x[0] |
---|
56 | ; if loc[i] eq nx2-1: x2[i] >= x[nx-1] |
---|
57 | ; else : x[loc[i]] <= x2[i] < x[loc[i]+1] |
---|
58 | loc = value_locate(x, x2) |
---|
59 | ; change loc in order to |
---|
60 | ; use x[0] and x[1] even if x2[i] < x[0] -> extrapolation |
---|
61 | ; use x[nx-2] and x[nx-1] even if x2[i] >= x[nx-1] -> extrapolation |
---|
62 | loc = 0 > temporary(loc) < (nx-2) |
---|
63 | ; distance between to consecutive x |
---|
64 | deltax = x[loc+1]-x[loc] |
---|
65 | ; distance between to consecutive y |
---|
66 | deltay = y[loc+1]-y[loc] |
---|
67 | ; relative distance between x2[i] and x[loc[i]+1] |
---|
68 | a = (x[loc+1]-x2)/deltax |
---|
69 | ; relative distance between x2[i] and x[loc[i]] |
---|
70 | b = 1.0d - a |
---|
71 | ; compute the first derivative on x (see numerical recipes Chap 3.3) |
---|
72 | yfrst = temporary(deltay)/deltax $ |
---|
73 | - 1.0d/6.0d * (3.0d*a*a - 1.0d) * deltax * yscd[loc] $ |
---|
74 | + 1.0d/6.0d * (3.0d*b*b - 1.0d) * deltax * yscd[loc+1] |
---|
75 | ; beware of the computation precision... |
---|
76 | ; force near zero values to be exactly 0.0 |
---|
77 | zero = where(abs(yfrst) LT 1.e-10) |
---|
78 | IF zero[0] NE -1 THEN yfrst[zero] = 0.0d |
---|
79 | |
---|
80 | RETURN, yfrst |
---|
81 | END |
---|