Context Navigation

spl_incr.pro @ 136

Last change on this file since 136 was 136, checked in by pinsard, 18 years ago
some improvements and corrections in some .pro file according to aspell and idldoc log file
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	;------------------------------------------------------------
70	FUNCTION 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
87	END
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	;-
108	FUNCTION 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
126	END
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	;-
139	FUNCTION spl_incr, x, y, x2, YP0 = yp0, YPN_1 = ypn_1
140	;
141	compile_opt idl2, strictarrsubs
142	;
143	;---------------------------------
144	; check and initialisation ...
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 = aX^3 + bX^2 + c*X + d
330	; y' = 3aX^2 + 2bX + 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 = bbbbbb - 3.0daaa*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]xinflxinflxinfl + 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]xzeroxzeroxzero + 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]xzeroxzeroxzero + 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
486	convexconcave:
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	;
524	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format