source: trunk/SRC/Interpolation/spl_incr.pro @ 228

Last change on this file since 228 was 163, checked in by navarro, 18 years ago

header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)

  • Property svn:keywords set to Id
File size: 19.8 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5;
6; @file_comments
7;
8; Given the arrays X and Y, which tabulate a function (with the X[i]
9; AND Y[i] in ascending order), and given an input value X2, the
10; SPL_INCR function returns an interpolated value for the given values
11; of X2. The interpolation method is based on cubic spline, corrected
12; in a way that interpolated values are also monotonically increasing.
13;
14; @param x1 {in}{required}
15; An n-element (at least 2) input vector that specifies the tabulate points in
16; a strict ascending order.
17;
18; @param y1 {in}{required}
19; f(x) = y. An n-element input vector that specifies the values
20; of the tabulated function F(Xi) corresponding to Xi. As f is
21; supposed to be monotonically increasing, y values must be
22; monotonically increasing. y can have equal consecutive values.
23;
24; @param x2 {in}{required}
25; The input values for which the interpolated values are
26; desired. Its values must be strictly monotonically increasing.
27;
28; @param der2
29;
30; @param x
31;
32; @returns
33;
34;    y2: f(x2) = y2. Double precision array
35;
36; @restrictions
37; It might be possible that y2[i+1]-y2[i] has very small negative
38; values (amplitude smaller than 1.e-6)...
39;
40; @examples
41;
42; IDL> n = 100L
43; IDL> x = (dindgen(n))^2
44; IDL> y = abs(randomn(0, n))
45; IDL> y[n/2:n/2+1] = 0.
46; IDL> y[n-n/3] = 0.
47; IDL> y[n-n/6:n-n/6+5] = 0.
48; IDL> y = total(y, /cumulative, /double)
49; IDL> x2 = dindgen((n-1)^2)
50; IDL> n2 = n_elements(x2)
51; IDL> print, min(y[1:n-1]-y[0:n-2]) LT 0
52; IDL> y2 = spl_incr( x, y, x2)
53; IDL> splot, x, y, xstyle = 1, ystyle = 1, ysurx=.25, petit = [1, 2, 1], /land
54; IDL> oplot, x2, y2, color = 100
55; IDL> c = y2[1:n2-1] - y2[0:n2-2]
56; IDL> print, min(c) LT 0
57; IDL> print, min(c, max = ma), ma
58; IDL> splot,c,xstyle=1,ystyle=1, yrange=[-.01,.05], ysurx=.25, petit = [1, 2, 2], /noerase
59; IDL> oplot,[0, n_elements(c)], [0, 0], linestyle = 1
60;
61; @history
62;  Sebastien Masson (smasson\@lodyc.jussieu.fr): May-Dec 2005
63;
64; @version $Id$
65;
66;-
67;------------------------------------------------------------
68;------------------------------------------------------------
69;------------------------------------------------------------
70FUNCTION pure_concave, x1, x2, y1, y2, der2, x
71; X^n type
72;
73  compile_opt idl2, strictarrsubs
74;
75  xx = (double(x)-double(x1))/(double(x2)-double(x1))
76  f = (double(x2)-double(x1))/(double(y2)-double(y1))
77  n = der2*temporary(f)
78  res = xx^(n)
79;   IF check_math() GT 0 THEN BEGIN
80;       zero = where(abs(res) LT 1.e-10)
81;       IF zero[0] NE -1 THEN res[zero] = 0.0d
82;   END
83  res = temporary(res)*(double(y2)-double(y1))+y1
84;
85;  IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop
86  RETURN, res
87END
88
89;+
90; @param x1 {in}{required}
91; An n-element (at least 2) input vector that specifies the tabulate points in
92; a strict ascending order.
93;
94; @param y1 {in}{required}
95; f(x) = y. An n-element input vector that specifies the values
96;    of the tabulated function F(Xi) corresponding to Xi. As f is
97;    supposed to be monotonically increasing, y values must be
98;    monotonically increasing. y can have equal consecutive values.
99;
100; @param x2 {in}{required}
101; The input values for which the interpolated values are
102; desired. Its values must be strictly monotonically increasing.
103;
104; @param der2
105; @param x
106;
107;-
108FUNCTION pure_convex, x1, x2, y1, y2, der2, x
109; 1-(1-X)^n type
110;
111  compile_opt idl2, strictarrsubs
112;
113  xx = 1.0d - (double(x)-double(x1))/(double(x2)-double(x1))
114  f = (double(x2)-double(x1))/(double(y2)-double(y1))
115  n = der2*temporary(f)
116  res = xx^(n)
117;   IF check_math() GT 0 THEN BEGIN
118;       zero = where(abs(res) LT 1.e-10)
119;       IF zero[0] NE -1 THEN res[zero] = 0.0d
120;   END
121  res = 1.0d - temporary(res)
122  res = temporary(res)*(y2-y1)+y1
123;
124;  IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop
125  RETURN, res
126END
127
128;+
129; @param x
130; @param y
131; @param x2
132; @keyword YP0 The first derivative of the interpolating function at the
133;    point X0. If YP0 is omitted, the second derivative at the
134;    boundary is set to zero, resulting in a "natural spline."
135; @keyword YPN_1 The first derivative of the interpolating function at the
136;    point Xn-1. If YPN_1 is omitted, the second derivative at the
137;    boundary is set to zero, resulting in a "natural spline."
138;-
139FUNCTION spl_incr, x, y, x2, YP0 = yp0, YPN_1 = ypn_1
140;
141  compile_opt idl2, strictarrsubs
142;
143;---------------------------------
144; check and initialization ...
145;---------------------------------
146;
147  nx = n_elements(x)
148  ny = n_elements(y)
149  nx2 = n_elements(x2)
150; x must have at least 2 elements
151  IF nx LT 2 THEN stop
152; y must have the same number of elements than x
153  IF nx NE ny THEN stop
154; x be monotonically increasing
155  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop
156; x2 be monotonically increasing
157  IF N_ELEMENTS(X2) GE 2 THEN $
158  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop
159; y be monotonically increasing
160  IF min(y[1:ny-1]-y[0:ny-2]) LT 0 THEN stop
161;---------------------------------
162; first check: check if two consecutive values are equal
163;---------------------------------
164  bad = where(y[1:ny-1]-y[0:ny-2] EQ 0, cntbad)
165  IF cntbad NE 0 THEN BEGIN
166; define the results: y2
167      y2 = dblarr(nx2)
168; define xinx2: see help of value_locate
169;  if xinx2[i] eq -1   :                 x[bad[i]] <  x2[0]
170;  if xinx2[i] eq nx2-1:                 x[bad[i]] >= x2[nx2-1]
171;  else                : x2[xinx2[i]] <= x[bad[i]] <  x2[xinx2[i]+1]
172    xinx2 = value_locate(x2, x[bad])
173    xinx2_1 = value_locate(x2, x[bad+1])
174;
175; left side ... if there is x2 values smaller that x[bad[0]].
176; we force ypn_1 = 0.0d
177    IF xinx2[0] NE -1 THEN BEGIN
178      IF bad[0] EQ 0 THEN BEGIN
179        IF xinx2[0] NE 0 THEN stop
180        y2[0] = y[0]
181      ENDIF ELSE BEGIN
182        y2[0:xinx2[0]] $
183          = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $
184                     , yp0 = yp0, ypn_1 = 0.0d)
185      ENDELSE
186    ENDIF
187; flat section
188    IF xinx2_1[0] NE -1 THEN $
189      y2[(xinx2[0]+1) < xinx2_1[0] : xinx2_1[0]] = y[bad[0]]
190; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in
191; more than 2 pieces...
192      IF cntbad GT 1 THEN BEGIN
193; we take care of the piece located wetween bad[ib-1]+1 and bad[ib]
194        FOR ib = 1, cntbad-1 DO BEGIN
195; if there is x2 values smaller that x[bad[ib]], then the x2 values
196; located between bad[ib-1]+1 and bad[ib] are (xinx2_1[ib-1]+1:xinx2[ib]
197; and if we don't have two consecutive flat sections
198          IF xinx2[ib] NE -1 AND (bad[ib-1] NE bad[ib]-1) THEN begin
199            y2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
200              = spl_incr(x[bad[ib-1]+1:bad[ib]], y[bad[ib-1]+1:bad[ib]] $
201                         , x2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
202                         , yp0 = 0.0d, ypn_1 = 0.0d)
203          ENDIF
204; flat section
205          IF xinx2_1[ib] NE -1 THEN $
206            y2[(xinx2[ib]+1) < xinx2_1[ib] : xinx2_1[ib]] = y[bad[ib]]
207        ENDFOR
208      ENDIF
209; right side ... if there is x2 values larger that x[bad[cntbad-1]+1].
210; we force yp0 = 0.0d
211      IF xinx2_1[cntbad-1] NE nx2-1 THEN $
212        y2[xinx2_1[cntbad-1]+1:nx2-1] $
213        = spl_incr(x[bad[cntbad-1]+1:nx-1], y[bad[cntbad-1]+1:nx-1] $
214                        , x2[xinx2_1[cntbad-1]+1:nx2-1] $
215                        , yp0 = 0.0d, ypn_1 = ypn_1new)
216
217    RETURN, y2
218
219  ENDIF
220;-----------
221; compute the second derivative of the cubic spline on each x.
222;-----------
223  yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double)
224;---------------------------------
225; second check: none of the first derivative on x values must be negative.
226;---------------------------------
227;
228; compute the first derivative on x
229;
230  yifrst = spl_fstdrv(x, y, yscd, x)
231;
232; we force the negative first derivative to 0 by calling again
233; spl_incr with the keywords yp0 and ypn_1 to specify the
234; first derivative equal to 0
235;
236  bad = where(yifrst LT 0.0d, cntbad)
237  IF cntbad NE 0 THEN BEGIN
238;
239; we define the new values of the keyword ypn_1:
240; if the first derivative of the last value of x is negative
241; we define the new values of the keyword ypn_1 to 0.0d0
242    IF bad[cntbad-1] EQ nx-1 THEN BEGIN
243      ypn_1new = 0.0d
244; we remove this case from the list
245      IF cntbad GE 2 THEN bad = bad[0:cntbad-2]
246      cntbad = cntbad-1
247; else we take the value of ypn_1 if it was already defined
248    ENDIF ELSE IF n_elements(ypn_1) NE 0 THEN ypn_1new = ypn_1
249;
250; we define the new values of the keyword yp0:
251; if the first derivative of the first value of x is negative
252; we define the new values of the keyword yp0 to 0.0
253    IF bad[0] EQ 0 THEN BEGIN
254      yp0new = 0.0d
255; we remove this case from the list
256      IF cntbad GE 2 THEN bad = bad[1:cntbad-1]
257      cntbad = cntbad-1
258; else we take the value of yp0 if it was already defined
259    ENDIF ELSE IF n_elements(yp0) NE 0 THEN yp0new = yp0
260;
261; if all the negative derivative corresponded to one of the cases above,
262; then we can directly call spl_incr with the new yp0new and ypn_1new
263    IF cntbad LE 0 THEN BEGIN
264      y2 = spl_incr(x, y, x2, yp0 = yp0new, ypn_1 = ypn_1new)
265;
266; else: there is still cases with negative derivative ...
267; we will cut spl_incr in n spl_incr and specify yp0, ypn_1
268; for each of this n spl_incr
269    ENDIF ELSE BEGIN
270; define xinx2: see help of value_locate
271;  if xinx2[i] eq -1   :                 x[bad[i]] <  x2[0]
272;  if xinx2[i] eq nx2-1:                 x[bad[i]] >= x2[nx2-1]
273;  else                : x2[xinx2[i]] <= x[bad[i]] <  x2[xinx2[i]+1]
274      xinx2 = value_locate(x2, x[bad])
275      y2 = dblarr(nx2)
276; left side ... if there is x2 values smaller that x[bad[0]].
277; we force ypn_1 = 0.0d
278      IF xinx2[0] NE -1 THEN $
279        y2[0:xinx2[0]] $
280        = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $
281                        , yp0 = yp0new, ypn_1 = 0.0d)
282; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in
283; more than 2 pieces -> we have middle pieces for which
284; we force yp0 = 0.0d and ypn_1 = 0.0d
285      IF cntbad GT 1 THEN BEGIN
286; we take care of the piece located wetween bad[ib-1] and bad[ib]
287        FOR ib = 1, cntbad-1 DO BEGIN
288; if there is x2 values smaller that x[bad[ib]], then the x2 values
289; located between bad[ib-1] and bad[ib] are (xinx2[ib-1]+1:xinx2[ib]
290          IF xinx2[ib] NE -1 THEN begin
291            y2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
292              = spl_incr(x[bad[ib-1]:bad[ib]], y[bad[ib-1]:bad[ib]] $
293                              , x2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
294                              , yp0 = 0.0d, ypn_1 = 0.0d)
295          endif
296        ENDFOR
297      ENDIF
298; right side ... if there is x2 values larger that x[bad[cntbad-1]].
299; we force yp0 = 0.0d
300      IF xinx2[cntbad-1] NE nx2-1 THEN $
301        y2[xinx2[cntbad-1]+1:nx2-1] $
302        = spl_incr(x[bad[cntbad-1]:nx-1], y[bad[cntbad-1]:nx-1] $
303                        , x2[xinx2[cntbad-1]+1:nx2-1] $
304                        , yp0 = 0.0d, ypn_1 = ypn_1new)
305    ENDELSE
306; we return the checked and corrected value of yfrst
307;       FOR i = 0, nx-1 DO BEGIN
308;         same = where(abs(x2- x[i]) LT 1.e-10, cnt)
309; ;        IF cnt NE 0 THEN y2[same] = y[i]
310;       ENDFOR
311    RETURN, y2
312  ENDIF
313;
314; we can be in this part of the code only if:
315;  (1) spl_incr is called by itself
316;  (2) none are the first derivative in x are negative (because they have been
317;      checked and corrected by the previous call to spl_incr, see above)
318;---------------------------------
319; third check: we have to make sure that the first derivative cannot
320;               have negative values between on x[0] and x[nx-1]
321;---------------------------------
322;
323; first we compute the first derivative, next we correct the values
324; where we know that the first derivative can be negative.
325;
326  y2 = spl_interp(x, y, yscd, x2, /double)
327;
328; between x[i] and x[i+1], the cubic spline is a cubic function:
329; y  =  a*X^3 +  b*X^2 + c*X + d
330; y' = 3a*X^2 + 2b*X   + c
331; y''= 6a*X   + 2b
332; if we take X = x[i+1]-x[i] then
333;    d = y[i]; c = y'[i]; b = 0.5 * y''[i],
334;    a = 1/6 * (y''[i+1]-y''[i])/(x[i+1]-x[i])
335;
336; y'[i] and y'[i+1] are positive so y' can be negative
337; between x[i] and x[i+1] only if
338;   1) a > 0
339;            ==> y''[i+1] > y''[i]
340;   2) y' reach its minimum value between  x[i] and x[i+1]
341;      -> 0 < - b/(3a) < x[i+1]-x[i]
342;            ==> y''[i+1] > 0 > y''[i]
343;
344; we do a first selection by looking for those points...
345;
346  loc = lindgen(nx-1)
347  maybebad = where(yscd[loc] LE 0.0d AND yscd[loc+1] GE 0.0d, cntbad)
348;
349  IF cntbad NE 0 THEN BEGIN
350
351    mbbloc = loc[maybebad]
352
353    aaa = (yscd[mbbloc+1]-yscd[mbbloc])/(6.0d*(x[mbbloc+1]-x[mbbloc]))
354    bbb = 0.5d * yscd[mbbloc]
355    ccc = yifrst[mbbloc]
356    ddd = y[mbbloc]
357;
358; definitive selection:
359; y' can become negative if and only if (2b)^2 - 4(3a)c > 0
360; y' can become negative if and only if    b^2  - (3a)c > 0
361;
362    delta = bbb*bbb - 3.0d*aaa*ccc
363;
364    bad = where(delta GT 0, cntbad)
365;
366    IF cntbad NE 0 THEN BEGIN
367      delta = delta[bad]
368      aaa = aaa[bad]
369      bbb = bbb[bad]
370      ccc = ccc[bad]
371      ddd = ddd[bad]
372      bad = maybebad[bad]
373; define xinx2_1: see help of value_locate
374;  if xinx2_1[i] eq -1   :                   x[bad[i]] <  x2[0]
375;  if xinx2_1[i] eq nx2-1:                   x[bad[i]] >= x2[nx2-1]
376;  else                  : x2[xinx2_1[i]] <= x[bad[i]] <  x2[xinx2_1[i]+1]
377      xinx2_1 = value_locate(x2, x[bad])
378; define xinx2_2: see help of value_locate
379;  if xinx2_2[i] eq -1   :                   x[bad[i]+1] <  x2[0]
380;  if xinx2_2[i] eq nx2-1:                   x[bad[i]+1] >= x2[nx2-1]
381;  else                  : x2[xinx2_2[i]] <= x[bad[i]+1] <  x2[xinx2_2[i]+1]
382      xinx2_2 = value_locate(x2, x[bad+1])
383; to avoid the particular case when x2 = x[bad[i]]
384; and there is no other x2 point until x[bad[i]+1]
385      xinx2_1 = xinx2_1 < (xinx2_2-1)
386;
387      FOR ib = 0, cntbad-1 DO BEGIN
388;
389; at least one of the x2 points must be located between
390; x[bad[ib]] and x[bad[ib]+1]
391        IF x2[0] LE x[bad[ib]+1] AND x2[nx2-1] GE x[bad[ib]] THEN BEGIN
392;
393          CASE 1 OF
394            yifrst[bad[ib]+1] EQ 0.0d:BEGIN
395; case pur convex: we use the first derivative of 1-(1-x)^n
396; and ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]]
397              y2[xinx2_1[ib]+1:xinx2_2[ib]] $
398                = pure_convex(x[bad[ib]], x[bad[ib]+1] $
399                              , y[bad[ib]], y[bad[ib]+1]  $
400                              , yifrst[bad[ib]] $
401                              , x2[xinx2_1[ib]+1:xinx2_2[ib]])
402            END
403            yifrst[bad[ib]] EQ 0.0d:BEGIN
404; case pur concave: we use the first derivative of x^n
405; and ajust n to get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1]
406              y2[xinx2_1[ib]+1:xinx2_2[ib]] $
407                = pure_concave(x[bad[ib]], x[bad[ib]+1] $
408                               , y[bad[ib]], y[bad[ib]+1] $
409                               , yifrst[bad[ib]+1] $
410                               , x2[xinx2_1[ib]+1:xinx2_2[ib]])
411            END
412            ELSE:BEGIN
413; in those cases, the first derivative has 2 zero between
414; x[bad[ib]] and x[bad[ib]+1]. We look for the minimum value of the
415; first derivative that correspond to the inflection point of y
416              xinfl = -bbb[ib]/(3.0d*aaa[ib])
417; we compute the y value for xinfl
418              yinfl = aaa[ib]*xinfl*xinfl*xinfl + bbb[ib]*xinfl*xinfl $
419                + ccc[ib]*xinfl + ddd[ib]
420;
421              CASE 1 OF
422; if y[xinfl] smaller than y[bad[ib]] then we conserve y2 until
423; the first zero of y2 and from this point we use x^n and ajust n to
424; get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1]
425                yinfl LT y[bad[ib]]:BEGIN
426; value of the first zero (y'[xzero]=0)
427                  xzero = (-bbb[ib]-sqrt(delta[ib]))/(3.0d*aaa[ib])
428; value of y[xzero]...
429                  yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $
430                    + ccc[ib]*xzero + ddd[ib]
431; if yzero > y[bad[ib]+1] then we cannot applay the method we want to
432; apply => we use then convex-concave case by changing by hand the
433; value of yinfl and xinfl
434                  IF yzero GT y[bad[ib]+1] THEN BEGIN
435                    yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]])
436                    xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]])
437                    GOTO, convexconcave
438                  ENDIF
439; define xinx2_3: see help of value_locate
440;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
441;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
442;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
443                  xinx2_3 = value_locate(x2, x[bad[ib]]+xzero)
444; to avoid the particular case when x2 = x[bad[ib]]+xzero
445; and there is no other x2 point until x[bad[ib]+1]
446                  xinx2_3 = xinx2_3 < (xinx2_2[ib]-1)
447                  IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN
448                    y2[xinx2_3+1:xinx2_2[ib]] $
449                      = pure_concave(x[bad[ib]]+xzero, x[bad[ib]+1] $
450                                     , yzero, y[bad[ib]+1] $
451                                     , yifrst[bad[ib]+1] $
452                                     , x2[xinx2_3+1:xinx2_2[ib]])
453                  ENDIF
454                END
455; if y[xinfl] bigger than y[bad[ib]+1] then we conserve y2 from
456; the second zero of y2 and before this point we use 1-(1-x)^n and
457; ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]]
458                yinfl GT y[bad[ib]+1]:BEGIN
459; value of the second zero (y'[xzero]=0)
460                  xzero = (-bbb[ib]+sqrt(delta[ib]))/(3.0d*aaa[ib])
461; value of y[xzero]...
462                  yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $
463                    + ccc[ib]*xzero + ddd[ib]
464; if yzero < y[bad[ib]] then we cannot applay the method we want to
465; apply => we use then convex-concave case by changing by hand the
466; value of yinfl and xinfl
467                  IF yzero lt y[bad[ib]] THEN BEGIN
468                    yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]])
469                    xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]])
470                    GOTO, convexconcave
471                  ENDIF
472; define xinx2_3: see help of value_locate
473;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
474;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
475;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
476                  xinx2_3 = value_locate(x2, x[bad[ib]]+xzero)
477                  IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN
478                    y2[xinx2_1[ib]+1:xinx2_3] $
479                      = pure_convex(x[bad[ib]], x[bad[ib]]+xzero  $
480                                    , y[bad[ib]], yzero   $
481                                    , yifrst[bad[ib]] $
482                                    , x2[xinx2_1[ib]+1:xinx2_3])
483                  ENDIF
484                END
485                ELSE:BEGIN
486convexconcave:
487; define xinx2_3: see help of value_locate
488;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
489;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
490;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
491                  xinx2_3 = value_locate(x2, x[bad[ib]]+xinfl)
492
493                  IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN
494                    y2[xinx2_1[ib]+1:xinx2_3] $
495                      = pure_convex(x[bad[ib]], x[bad[ib]]+xinfl  $
496                                    , y[bad[ib]], yinfl  $
497                                    , yifrst[bad[ib]] $
498                                    , x2[xinx2_1[ib]+1:xinx2_3])
499
500                  ENDIF
501                  IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN
502                    y2[xinx2_3+1:xinx2_2[ib]] $
503                      = pure_concave(x[bad[ib]]+xinfl, x[bad[ib]+1] $
504                                     , yinfl, y[bad[ib]+1] $
505                                     , yifrst[bad[ib]+1] $
506                                     , x2[xinx2_3+1:xinx2_2[ib]])
507                  ENDIF
508                END
509              ENDCASE
510
511            END
512          ENDCASE
513        ENDIF
514      ENDFOR
515
516    ENDIF
517  ENDIF
518;
519  RETURN, y2
520;
521;------------------------------------------------------------------
522;------------------------------------------------------------------
523;
524END
Note: See TracBrowser for help on using the repository browser.