;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ;+ ; ; @file_comments ; ; Given the arrays X and Y, which tabulate a function (with the X[i] ; AND Y[i] in ascending order), and given an input value X2, the ; SPL_INCR function returns an interpolated value for the given values ; of X2. The interpolation method is based on cubic spline, corrected ; in a way that interpolated values are also monotonically increasing. ; ; @examples y2 = spl_incr(x, y, x2) ; ; @param x1 {in}{required} An n-element (at least 2) input vector that specifies the ; tabulate points in a strict ascending order. ; ; @param y1 {in}{required} f(x) = y. An n-element input vector that specifies the values ; of the tabulated function F(Xi) corresponding to Xi. As f is ; supposed to be monotonically increasing, y values must be ; monotonically increasing. y can have equal consecutive values. ; ; @param x2 {in}{required} The input values for which the interpolated values are ; desired. Its values must be strictly monotonically increasing. ; ; ; ; ; @returns ; ; y2: f(x2) = y2. Double precision array ; ; @restrictions ; It might be possible that y2[i+1]-y2[i] has very small negative ; values (amplitude smaller than 1.e-6)... ; ; @examples ; ; n = 100L ; x = (dindgen(n))^2 ; y = abs(randomn(0, n)) ; y[n/2:n/2+1] = 0. ; y[n-n/3] = 0. ; y[n-n/6:n-n/6+5] = 0. ; y = total(y, /cumulative, /double) ; x2 = dindgen((n-1)^2) ; n2 = n_elements(x2) ; print, min(y[1:n-1]-y[0:n-2]) LT 0 ; y2 = spl_incr( x, y, x2) ; splot, x, y, xstyle = 1, ystyle = 1, ysurx=.25, petit = [1, 2, 1], /land ; oplot, x2, y2, color = 100 ; c = y2[1:n2-1] - y2[0:n2-2] ; print, min(c) LT 0 ; print, min(c, max = ma), ma ; splot,c,xstyle=1,ystyle=1, yrange=[-.01,.05], ysurx=.25, petit = [1, 2, 2], /noerase ; oplot,[0, n_elements(c)], [0, 0], linestyle = 1 ; ; @history ; Sebastien Masson (smasson\@lodyc.jussieu.fr): May-Dec 2005 ;- ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ FUNCTION pure_concave, x1, x2, y1, y2, der2, x ; X^n type xx = (double(x)-double(x1))/(double(x2)-double(x1)) f = (double(x2)-double(x1))/(double(y2)-double(y1)) n = der2*temporary(f) res = xx^(n) ; IF check_math() GT 0 THEN BEGIN ; zero = where(abs(res) LT 1.e-10) ; IF zero[0] NE -1 THEN res[zero] = 0.0d ; END res = temporary(res)*(double(y2)-double(y1))+y1 ; ; IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop RETURN, res END FUNCTION pure_convex, x1, x2, y1, y2, der2, x ; 1-(1-X)^n type xx = 1.0d - (double(x)-double(x1))/(double(x2)-double(x1)) f = (double(x2)-double(x1))/(double(y2)-double(y1)) n = der2*temporary(f) res = xx^(n) ; IF check_math() GT 0 THEN BEGIN ; zero = where(abs(res) LT 1.e-10) ; IF zero[0] NE -1 THEN res[zero] = 0.0d ; END res = 1.0d - temporary(res) res = temporary(res)*(y2-y1)+y1 ; ; IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop RETURN, res END ;+ ; @keyword YP0 The first derivative of the interpolating function at the ; point X0. If YP0 is omitted, the second derivative at the ; boundary is set to zero, resulting in a "natural spline." ; @keyword YPN_1 The first derivative of the interpolating function at the ; point Xn-1. If YPN_1 is omitted, the second derivative at the ; boundary is set to zero, resulting in a "natural spline." ;- FUNCTION spl_incr, x, y, x2, YP0 = yp0, YPN_1 = ypn_1 ; ;--------------------------------- ; check and initialisation ... ;--------------------------------- nx = n_elements(x) ny = n_elements(y) nx2 = n_elements(x2) ; x must have at least 2 elements IF nx LT 2 THEN stop ; y must have the same number of elements than x IF nx NE ny THEN stop ; x be monotonically increasing IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop ; x2 be monotonically increasing IF N_ELEMENTS(X2) GE 2 THEN $ IF min(x2[1:nx2-1]-x2[0:nx2-2]) LE 0 THEN stop ; y be monotonically increasing IF min(y[1:ny-1]-y[0:ny-2]) LT 0 THEN stop ;--------------------------------- ; first check: check if two consecutive values are equal ;--------------------------------- bad = where(y[1:ny-1]-y[0:ny-2] EQ 0, cntbad) IF cntbad NE 0 THEN BEGIN ; define the results: y2 y2 = dblarr(nx2) ; define xinx2: see help of value_locate ; if xinx2[i] eq -1 : x[bad[i]] < x2[0] ; if xinx2[i] eq nx2-1: x[bad[i]] >= x2[nx2-1] ; else : x2[xinx2[i]] <= x[bad[i]] < x2[xinx2[i]+1] xinx2 = value_locate(x2, x[bad]) xinx2_1 = value_locate(x2, x[bad+1]) ; ; left side ... if there is x2 values smaller that x[bad[0]]. ; we force ypn_1 = 0.0d IF xinx2[0] NE -1 THEN BEGIN IF bad[0] EQ 0 THEN BEGIN IF xinx2[0] NE 0 THEN stop y2[0] = y[0] ENDIF ELSE BEGIN y2[0:xinx2[0]] $ = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $ , yp0 = yp0, ypn_1 = 0.0d) ENDELSE ENDIF ; flat section IF xinx2_1[0] NE -1 THEN $ y2[(xinx2[0]+1) < xinx2_1[0] : xinx2_1[0]] = y[bad[0]] ; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in ; more than 2 pieces... IF cntbad GT 1 THEN BEGIN ; we take care of the piece located wetween bad[ib-1]+1 and bad[ib] FOR ib = 1, cntbad-1 DO BEGIN ; if there is x2 values smaller that x[bad[ib]], then the x2 values ; located between bad[ib-1]+1 and bad[ib] are (xinx2_1[ib-1]+1:xinx2[ib] ; and if we don't have two consecutive flat sections IF xinx2[ib] NE -1 AND (bad[ib-1] NE bad[ib]-1) THEN begin y2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $ = spl_incr(x[bad[ib-1]+1:bad[ib]], y[bad[ib-1]+1:bad[ib]] $ , x2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $ , yp0 = 0.0d, ypn_1 = 0.0d) ENDIF ; flat section IF xinx2_1[ib] NE -1 THEN $ y2[(xinx2[ib]+1) < xinx2_1[ib] : xinx2_1[ib]] = y[bad[ib]] ENDFOR ENDIF ; right side ... if there is x2 values larger that x[bad[cntbad-1]+1]. ; we force yp0 = 0.0d IF xinx2_1[cntbad-1] NE nx2-1 THEN $ y2[xinx2_1[cntbad-1]+1:nx2-1] $ = spl_incr(x[bad[cntbad-1]+1:nx-1], y[bad[cntbad-1]+1:nx-1] $ , x2[xinx2_1[cntbad-1]+1:nx2-1] $ , yp0 = 0.0d, ypn_1 = ypn_1new) RETURN, y2 ENDIF ;----------- ; compute the second derivative of the cubic spline on each x. ;----------- yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double) ;--------------------------------- ; second check: none of the first derivative on x values must be negative. ;--------------------------------- ; ; compute the first derivative on x ; yifrst = spl_fstdrv(x, y, yscd, x) ; ; we force the negative first derivative to 0 by calling again ; spl_incr with the keywords yp0 and ypn_1 to specify the ; first derivative equal to 0 ; bad = where(yifrst LT 0.0d, cntbad) IF cntbad NE 0 THEN BEGIN ; ; we define the new values of the keyword ypn_1: ; if the first derivative of the last value of x is negative ; we define the new values of the keyword ypn_1 to 0.0d0 IF bad[cntbad-1] EQ nx-1 THEN BEGIN ypn_1new = 0.0d ; we remove this case from the list IF cntbad GE 2 THEN bad = bad[0:cntbad-2] cntbad = cntbad-1 ; else we take the value of ypn_1 if it was already defined ENDIF ELSE IF n_elements(ypn_1) NE 0 THEN ypn_1new = ypn_1 ; ; we define the new values of the keyword yp0: ; if the first derivative of the first value of x is negative ; we define the new values of the keyword yp0 to 0.0 IF bad[0] EQ 0 THEN BEGIN yp0new = 0.0d ; we remove this case from the list IF cntbad GE 2 THEN bad = bad[1:cntbad-1] cntbad = cntbad-1 ; else we take the value of yp0 if it was already defined ENDIF ELSE IF n_elements(yp0) NE 0 THEN yp0new = yp0 ; ; if all the negative derivative corresponded to one of the cases above, ; then we can directly call spl_incr with the new yp0new and ypn_1new IF cntbad LE 0 THEN BEGIN y2 = spl_incr(x, y, x2, yp0 = yp0new, ypn_1 = ypn_1new) ; ; else: there is still cases with negative derivative ... ; we will cut spl_incr in n spl_incr and specify yp0, ypn_1 ; for each of this n spl_incr ENDIF ELSE BEGIN ; define xinx2: see help of value_locate ; if xinx2[i] eq -1 : x[bad[i]] < x2[0] ; if xinx2[i] eq nx2-1: x[bad[i]] >= x2[nx2-1] ; else : x2[xinx2[i]] <= x[bad[i]] < x2[xinx2[i]+1] xinx2 = value_locate(x2, x[bad]) y2 = dblarr(nx2) ; left side ... if there is x2 values smaller that x[bad[0]]. ; we force ypn_1 = 0.0d IF xinx2[0] NE -1 THEN $ y2[0:xinx2[0]] $ = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $ , yp0 = yp0new, ypn_1 = 0.0d) ; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in ; more than 2 pieces -> we have middle pieces for which ; we force yp0 = 0.0d and ypn_1 = 0.0d IF cntbad GT 1 THEN BEGIN ; we take care of the piece located wetween bad[ib-1] and bad[ib] FOR ib = 1, cntbad-1 DO BEGIN ; if there is x2 values smaller that x[bad[ib]], then the x2 values ; located between bad[ib-1] and bad[ib] are (xinx2[ib-1]+1:xinx2[ib] IF xinx2[ib] NE -1 THEN begin y2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $ = spl_incr(x[bad[ib-1]:bad[ib]], y[bad[ib-1]:bad[ib]] $ , x2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $ , yp0 = 0.0d, ypn_1 = 0.0d) endif ENDFOR ENDIF ; right side ... if there is x2 values larger that x[bad[cntbad-1]]. ; we force yp0 = 0.0d IF xinx2[cntbad-1] NE nx2-1 THEN $ y2[xinx2[cntbad-1]+1:nx2-1] $ = spl_incr(x[bad[cntbad-1]:nx-1], y[bad[cntbad-1]:nx-1] $ , x2[xinx2[cntbad-1]+1:nx2-1] $ , yp0 = 0.0d, ypn_1 = ypn_1new) ENDELSE ; we return the checked and corrected value of yfrst ; FOR i = 0, nx-1 DO BEGIN ; same = where(abs(x2- x[i]) LT 1.e-10, cnt) ; ; IF cnt NE 0 THEN y2[same] = y[i] ; ENDFOR RETURN, y2 ENDIF ; ; we can be in this part of the code only if: ; (1) spl_incr is called by itself ; (2) none are the first derivative in x are negative (because they have been ; checked and corrected by the previous call to spl_incr, see above) ;--------------------------------- ; third check: we have to make sure that the first derivative cannot ; have negative values between on x[0] and x[nx-1] ;--------------------------------- ; ; first we compute the first derivative, next we correct the values ; where we know that the first derivative can be negative. ; y2 = spl_interp(x, y, yscd, x2, /double) ; ; between x[i] and x[i+1], the cubic spline is a cubic function: ; y = a*X^3 + b*X^2 + c*X + d ; y' = 3a*X^2 + 2b*X + c ; y''= 6a*X + 2b ; if we take X = x[i+1]-x[i] then ; d = y[i]; c = y'[i]; b = 0.5 * y''[i], ; a = 1/6 * (y''[i+1]-y''[i])/(x[i+1]-x[i]) ; ; y'[i] and y'[i+1] are positive so y' can be negative ; between x[i] and x[i+1] only if ; 1) a > 0 ; ==> y''[i+1] > y''[i] ; 2) y' reach its minimum value between x[i] and x[i+1] ; -> 0 < - b/(3a) < x[i+1]-x[i] ; ==> y''[i+1] > 0 > y''[i] ; ; we do a first selection by looking for those points... ; loc = lindgen(nx-1) maybebad = where(yscd[loc] LE 0.0d AND yscd[loc+1] GE 0.0d, cntbad) ; IF cntbad NE 0 THEN BEGIN mbbloc = loc[maybebad] aaa = (yscd[mbbloc+1]-yscd[mbbloc])/(6.0d*(x[mbbloc+1]-x[mbbloc])) bbb = 0.5d * yscd[mbbloc] ccc = yifrst[mbbloc] ddd = y[mbbloc] ; ; definitive selection: ; y' can become negative if and only if (2b)^2 - 4(3a)c > 0 ; y' can become negative if and only if b^2 - (3a)c > 0 ; delta = bbb*bbb - 3.0d*aaa*ccc ; bad = where(delta GT 0, cntbad) ; IF cntbad NE 0 THEN BEGIN delta = delta[bad] aaa = aaa[bad] bbb = bbb[bad] ccc = ccc[bad] ddd = ddd[bad] bad = maybebad[bad] ; define xinx2_1: see help of value_locate ; if xinx2_1[i] eq -1 : x[bad[i]] < x2[0] ; if xinx2_1[i] eq nx2-1: x[bad[i]] >= x2[nx2-1] ; else : x2[xinx2_1[i]] <= x[bad[i]] < x2[xinx2_1[i]+1] xinx2_1 = value_locate(x2, x[bad]) ; define xinx2_2: see help of value_locate ; if xinx2_2[i] eq -1 : x[bad[i]+1] < x2[0] ; if xinx2_2[i] eq nx2-1: x[bad[i]+1] >= x2[nx2-1] ; else : x2[xinx2_2[i]] <= x[bad[i]+1] < x2[xinx2_2[i]+1] xinx2_2 = value_locate(x2, x[bad+1]) ; to avoid the particular case when x2 = x[bad[i]] ; and there is no other x2 point until x[bad[i]+1] xinx2_1 = xinx2_1 < (xinx2_2-1) ; FOR ib = 0, cntbad-1 DO BEGIN ; ; at least one of the x2 points must be located between ; x[bad[ib]] and x[bad[ib]+1] IF x2[0] LE x[bad[ib]+1] AND x2[nx2-1] GE x[bad[ib]] THEN BEGIN ; CASE 1 OF yifrst[bad[ib]+1] EQ 0.0d:BEGIN ; case pur convex: we use the first derivative of 1-(1-x)^n ; and ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]] y2[xinx2_1[ib]+1:xinx2_2[ib]] $ = pure_convex(x[bad[ib]], x[bad[ib]+1] $ , y[bad[ib]], y[bad[ib]+1] $ , yifrst[bad[ib]] $ , x2[xinx2_1[ib]+1:xinx2_2[ib]]) END yifrst[bad[ib]] EQ 0.0d:BEGIN ; case pur concave: we use the first derivative of x^n ; and ajust n to get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1] y2[xinx2_1[ib]+1:xinx2_2[ib]] $ = pure_concave(x[bad[ib]], x[bad[ib]+1] $ , y[bad[ib]], y[bad[ib]+1] $ , yifrst[bad[ib]+1] $ , x2[xinx2_1[ib]+1:xinx2_2[ib]]) END ELSE:BEGIN ; in those cases, the first derivative has 2 zero between ; x[bad[ib]] and x[bad[ib]+1]. We look for the minimum value of the ; first derivative that correspond to the inflection point of y xinfl = -bbb[ib]/(3.0d*aaa[ib]) ; we compute the y value for xinfl yinfl = aaa[ib]*xinfl*xinfl*xinfl + bbb[ib]*xinfl*xinfl $ + ccc[ib]*xinfl + ddd[ib] ; CASE 1 OF ; if y[xinfl] smaller than y[bad[ib]] then we conserve y2 until ; the first zero of y2 and from this point we use x^n and ajust n to ; get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1] yinfl LT y[bad[ib]]:BEGIN ; value of the first zero (y'[xzero]=0) xzero = (-bbb[ib]-sqrt(delta[ib]))/(3.0d*aaa[ib]) ; value of y[xzero]... yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $ + ccc[ib]*xzero + ddd[ib] ; if yzero > y[bad[ib]+1] then we cannot applay the method we want to ; apply => we use then convex-concave case by changing by hand the ; value of yinfl and xinfl IF yzero GT y[bad[ib]+1] THEN BEGIN yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]]) xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]]) GOTO, convexconcave ENDIF ; define xinx2_3: see help of value_locate ; if xinx2_3[ib] eq -1 : x[bad[ib]]+xzero < x2[0] ; if xinx2_3[ib] eq nx2-1: x[bad[ib]]+xzero >= x2[nx2-1] ; else : x2[xinx2_3] <= x[bad[ib]]+xzero < x2[xinx3_2+1] xinx2_3 = value_locate(x2, x[bad[ib]]+xzero) ; to avoid the particular case when x2 = x[bad[ib]]+xzero ; and there is no other x2 point until x[bad[ib]+1] xinx2_3 = xinx2_3 < (xinx2_2[ib]-1) IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN y2[xinx2_3+1:xinx2_2[ib]] $ = pure_concave(x[bad[ib]]+xzero, x[bad[ib]+1] $ , yzero, y[bad[ib]+1] $ , yifrst[bad[ib]+1] $ , x2[xinx2_3+1:xinx2_2[ib]]) ENDIF END ; if y[xinfl] bigger than y[bad[ib]+1] then we conserve y2 from ; the second zero of y2 and before this point we use 1-(1-x)^n and ; ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]] yinfl GT y[bad[ib]+1]:BEGIN ; value of the second zero (y'[xzero]=0) xzero = (-bbb[ib]+sqrt(delta[ib]))/(3.0d*aaa[ib]) ; value of y[xzero]... yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $ + ccc[ib]*xzero + ddd[ib] ; if yzero < y[bad[ib]] then we cannot applay the method we want to ; apply => we use then convex-concave case by changing by hand the ; value of yinfl and xinfl IF yzero lt y[bad[ib]] THEN BEGIN yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]]) xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]]) GOTO, convexconcave ENDIF ; define xinx2_3: see help of value_locate ; if xinx2_3[ib] eq -1 : x[bad[ib]]+xzero < x2[0] ; if xinx2_3[ib] eq nx2-1: x[bad[ib]]+xzero >= x2[nx2-1] ; else : x2[xinx2_3] <= x[bad[ib]]+xzero < x2[xinx3_2+1] xinx2_3 = value_locate(x2, x[bad[ib]]+xzero) IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN y2[xinx2_1[ib]+1:xinx2_3] $ = pure_convex(x[bad[ib]], x[bad[ib]]+xzero $ , y[bad[ib]], yzero $ , yifrst[bad[ib]] $ , x2[xinx2_1[ib]+1:xinx2_3]) ENDIF END ELSE:BEGIN convexconcave: ; define xinx2_3: see help of value_locate ; if xinx2_3[ib] eq -1 : x[bad[ib]]+xzero < x2[0] ; if xinx2_3[ib] eq nx2-1: x[bad[ib]]+xzero >= x2[nx2-1] ; else : x2[xinx2_3] <= x[bad[ib]]+xzero < x2[xinx3_2+1] xinx2_3 = value_locate(x2, x[bad[ib]]+xinfl) IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN y2[xinx2_1[ib]+1:xinx2_3] $ = pure_convex(x[bad[ib]], x[bad[ib]]+xinfl $ , y[bad[ib]], yinfl $ , yifrst[bad[ib]] $ , x2[xinx2_1[ib]+1:xinx2_3]) ENDIF IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN y2[xinx2_3+1:xinx2_2[ib]] $ = pure_concave(x[bad[ib]]+xinfl, x[bad[ib]+1] $ , yinfl, y[bad[ib]+1] $ , yifrst[bad[ib]+1] $ , x2[xinx2_3+1:xinx2_2[ib]]) ENDIF END ENDCASE END ENDCASE ENDIF ENDFOR ENDIF ENDIF ; RETURN, y2 ; ;------------------------------------------------------------------ ;------------------------------------------------------------------ ; END