[69] | 1 | ;+ |
---|
[232] | 2 | ; |
---|
[101] | 3 | ; @file_comments |
---|
[69] | 4 | ; Given the arrays X and Y, which tabulate a function (with the X[i] |
---|
| 5 | ; AND Y[i] in ascending order), and given an input value X2, the |
---|
[242] | 6 | ; <pro>spl_incr</pro> function returns an interpolated value for the given values |
---|
[69] | 7 | ; of X2. The interpolation method is based on cubic spline, corrected |
---|
| 8 | ; in a way that integral of the interpolated values is the same as the |
---|
[297] | 9 | ; integral of the input values. (-> for example to build daily data |
---|
[69] | 10 | ; from monthly mean and keep the monthly mean of the computed daily |
---|
[297] | 11 | ; data equal to the original values) |
---|
[69] | 12 | ; |
---|
[125] | 13 | ; @param x {in}{required} |
---|
[242] | 14 | ; An n-elements (at least 2) input vector that specifies the tabulate points in |
---|
[125] | 15 | ; a strict ascending order. |
---|
[69] | 16 | ; |
---|
[163] | 17 | ; @param yin {in}{required}{type=array} |
---|
[125] | 18 | ; an array with one element less than x. y[i] represents the |
---|
| 19 | ; mean value between x[i] and x[i+1]. if /GE0 is activated, y must |
---|
| 20 | ; have positive values. |
---|
[69] | 21 | ; |
---|
[125] | 22 | ; @param x2 {in}{required} |
---|
| 23 | ; The input values for which the interpolated values are desired. |
---|
| 24 | ; Its values must be strictly monotonically increasing. |
---|
[69] | 25 | ; |
---|
[125] | 26 | ; @keyword GE0 |
---|
| 27 | ; to force that y2 is always GE than 0. In that case, y must also be GE than 0. |
---|
[69] | 28 | ; |
---|
[125] | 29 | ; @keyword YP0 |
---|
| 30 | ; The first derivative of the interpolating function at the |
---|
| 31 | ; point X0. If YP0 is omitted, the second derivative at the |
---|
| 32 | ; boundary is set to zero, resulting in a "natural spline." |
---|
[69] | 33 | ; |
---|
[125] | 34 | ; @keyword YPN_1 |
---|
| 35 | ; The first derivative of the interpolating function at the |
---|
| 36 | ; point Xn-1. If YPN_1 is omitted, the second derivative at the |
---|
| 37 | ; boundary is set to zero, resulting in a "natural spline." |
---|
[69] | 38 | ; |
---|
[125] | 39 | ; @returns |
---|
[136] | 40 | ; y2: the mean value between two consecutive values of x2. This |
---|
| 41 | ; array has one element less than y2. y2 has double precision. |
---|
[69] | 42 | ; |
---|
[101] | 43 | ; @restrictions |
---|
[136] | 44 | ; It might be possible that y2 has very small negative values |
---|
| 45 | ; (amplitude smaller than 1.e-6)... |
---|
[69] | 46 | ; |
---|
[125] | 47 | ; @examples |
---|
[69] | 48 | ; |
---|
| 49 | ; 12 monthly values of precipitations into daily values: |
---|
| 50 | ; |
---|
[121] | 51 | ; IDL> yr1 = 1990 |
---|
| 52 | ; IDL> yr2 = 1992 |
---|
| 53 | ; IDL> nyr = yr2-yr1+1 |
---|
| 54 | ; IDL> n1 = 12*nyr+1 |
---|
| 55 | ; IDL> x = julday(1+findgen(n1), replicate(1, n1) $ |
---|
| 56 | ; IDL> , replicate(yr1, n1), fltarr(n1)) |
---|
| 57 | ; IDL> n2 = 365*nyr + total(leapyr(yr1+indgen(nyr))) + 1 |
---|
| 58 | ; IDL> x2 = julday(replicate(1, n2), 1+findgen(n2) $ |
---|
| 59 | ; IDL> , replicate(yr1, n2), fltarr(n2)) |
---|
| 60 | ; IDL> y = abs(randomn(0, n1-1)) |
---|
| 61 | ; IDL> y2 = spl_keep_mean(x, y, x2, /ge0) |
---|
[69] | 62 | |
---|
[121] | 63 | ; IDL> print, min(x, max = ma), ma |
---|
| 64 | ; IDL> print, min(x2, max = ma), ma |
---|
| 65 | ; IDL> print, vairdate([min(x, max = ma), ma]) |
---|
| 66 | ; IDL> print, total(y*(x[1:n1-1]-x[0:n1-2])) |
---|
| 67 | ; IDL> print, total(y2*(x2[1:n2-1]-x2[0:n2-2])) |
---|
[69] | 68 | ; |
---|
[101] | 69 | ; @history |
---|
| 70 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005 |
---|
[118] | 71 | ; |
---|
[231] | 72 | ; @version |
---|
| 73 | ; $Id$ |
---|
[118] | 74 | ; |
---|
[69] | 75 | ;- |
---|
[231] | 76 | ; |
---|
[69] | 77 | FUNCTION spl_keep_mean, x, yin, x2, YP0 = yp0, YPN_1 = ypn_1, GE0 = ge0 |
---|
| 78 | ; |
---|
[114] | 79 | compile_opt idl2, strictarrsubs |
---|
| 80 | ; |
---|
[69] | 81 | ;--------------------------------- |
---|
[163] | 82 | ; check and initialization ... |
---|
[69] | 83 | ;--------------------------------- |
---|
| 84 | ; |
---|
| 85 | nx = n_elements(x) |
---|
| 86 | ny = n_elements(yin) |
---|
| 87 | nx2 = n_elements(x2) |
---|
| 88 | ; x must have at least 2 elements |
---|
[125] | 89 | IF nx LT 2 THEN stop |
---|
[69] | 90 | ; x2 must have at least 2 elements |
---|
[125] | 91 | IF nx2 LT 2 THEN stop |
---|
[69] | 92 | ; x be monotonically increasing |
---|
[125] | 93 | IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop |
---|
[69] | 94 | ; x2 be monotonically increasing |
---|
[125] | 95 | IF min(x2[1:nx2-1]-x2[0:nx2-2]) LE 0 THEN stop |
---|
[69] | 96 | ; |
---|
| 97 | ;--------------------------------- |
---|
| 98 | ; compute the integral of y |
---|
| 99 | ;--------------------------------- |
---|
| 100 | ; if spl_keep_mean is called by the user (and not by itself) we must compute |
---|
| 101 | ; the integral of y. yin must have one element less than x |
---|
| 102 | IF nx NE ny+1 THEN stop |
---|
| 103 | y = double(yin)*double(x[1:nx-1]-x[0:nx-2]) |
---|
| 104 | y = [0.0d, temporary(y)] |
---|
| 105 | y = total(temporary(y), /cumulative, /double) |
---|
| 106 | ; |
---|
| 107 | ;--------------------------------- |
---|
| 108 | ; compute the "spline" interpolation |
---|
| 109 | ;--------------------------------- |
---|
| 110 | ; |
---|
| 111 | IF keyword_set(ge0) THEN BEGIN |
---|
| 112 | ; if the want that the interpolated values are always >= 0, we must |
---|
| 113 | ; have yin >= 0.0d |
---|
| 114 | IF min(yin) LT 0 THEN stop |
---|
| 115 | ; call spl_incr |
---|
| 116 | y2 = spl_incr(x, temporary(y), x2, yp0 = yp0, ypn_1 = ypn_1) |
---|
| 117 | ENDIF ELSE BEGIN |
---|
| 118 | yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double) |
---|
| 119 | y2 = spl_interp(x, y, temporary(yscd), x2, /double) |
---|
| 120 | ENDELSE |
---|
| 121 | ; ; |
---|
| 122 | ;--------------------------------- |
---|
| 123 | ; Compute the derivative of y |
---|
| 124 | ;--------------------------------- |
---|
| 125 | ; |
---|
| 126 | yfrst = (y2[1:nx2-1]-y2[0:nx2-2])/(x2[1:nx2-1]-x2[0:nx2-2]) |
---|
| 127 | ; it can happen that we have very small negative values (-1.e-6 for ex) |
---|
| 128 | ; yfrst = 0.0d > temporary(yfrst) |
---|
| 129 | RETURN, yfrst |
---|
| 130 | ; |
---|
| 131 | END |
---|