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