Context Navigation

spl_incr.pro @ 114

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

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

Download in other formats:

Original Format