[69] | 1 | ;+ |
---|
| 2 | ; |
---|
[125] | 3 | ; @file_comments |
---|
[238] | 4 | ; returns the values of the first derivative of |
---|
[186] | 5 | ; the interpolating function at the points X2i. It is a double |
---|
[69] | 6 | ; precision array. |
---|
| 7 | ; |
---|
| 8 | ; Given the arrays X and Y, which tabulate a function (with the X[i] |
---|
[242] | 9 | ; and Y[i] in ascending order), and given an input value X2, the |
---|
[238] | 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. |
---|
[69] | 13 | ; |
---|
[125] | 14 | ; @examples |
---|
[69] | 15 | ; |
---|
[371] | 16 | ; IDL> y2 = spl_fstdrv(x, y, yscd, x2) |
---|
| 17 | ; |
---|
[125] | 18 | ; @param x {in}{required} |
---|
[242] | 19 | ; An n-elements (at least 2) input vector that specifies the |
---|
[125] | 20 | ; tabulate points in ascending order. |
---|
[69] | 21 | ; |
---|
[125] | 22 | ; @param y {in}{required} |
---|
[242] | 23 | ; f(x) = y. An n-elements input vector that specifies the values |
---|
[125] | 24 | ; of the tabulated function F(Xi) corresponding to Xi. |
---|
[69] | 25 | ; |
---|
[125] | 26 | ; @param yscd {in}{required} |
---|
[262] | 27 | ; The output from <proidl>SPL_INIT</proidl> for the specified X and Y. |
---|
[69] | 28 | ; |
---|
[242] | 29 | ; @param x2 {in}{required} {type= scalar or array} |
---|
[125] | 30 | ; The input values for which the first derivative values are desired. |
---|
[186] | 31 | ; |
---|
[125] | 32 | ; @returns |
---|
[372] | 33 | ; y2: f'(x2) = y2. |
---|
[69] | 34 | ; |
---|
[101] | 35 | ; @history |
---|
[372] | 36 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005 |
---|
[118] | 37 | ; |
---|
[231] | 38 | ; @version |
---|
| 39 | ; $Id$ |
---|
[118] | 40 | ; |
---|
[69] | 41 | ;- |
---|
| 42 | FUNCTION spl_fstdrv, x, y, yscd, x2 |
---|
| 43 | ; |
---|
[238] | 44 | compile_opt idl2, strictarrsubs |
---|
| 45 | ; |
---|
[114] | 46 | ; compute the first derivative of the spline function |
---|
[69] | 47 | ; |
---|
| 48 | nx = n_elements(x) |
---|
| 49 | ny = n_elements(y) |
---|
| 50 | ; x must have at least 2 elements |
---|
[125] | 51 | IF nx LT 2 THEN stop |
---|
[69] | 52 | ; y must have the same number of elements than x |
---|
| 53 | IF nx NE ny THEN stop |
---|
[125] | 54 | ; define loc in a way that |
---|
[69] | 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) |
---|
[125] | 59 | ; change loc in order to |
---|
[69] | 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 |
---|
[125] | 62 | loc = 0 > temporary(loc) < (nx-2) |
---|
[69] | 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] |
---|
[125] | 75 | ; beware of the computation precision... |
---|
[69] | 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 |
---|