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