Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

spl_incr.pro @ 98

Last change on this file since 98 was 69, checked in by smasson, 18 years ago
debug + new xxx
Property svn:executable set to ``*
File size: 19.3 KB

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

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

Download in other formats:

Original Format