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