SUBROUTINE MIN_QUADR_POLIN(AL0,AL1,F0,F1,DF0,A1,A2,X_MIN) C C---Given f(a),f'(a),f(b) a1,a2 finds minimum of quadratic interpolating C---polynom in the interval (min(a1,a2),max(a1,a2) IMPLICIT NONE REAL AL0,Al1,F0,F1,DF0,A1,A2,X_MIN C REAL A,B,C,F_X0,F_X1,X0,X1,AL1M0 REAL DF10,A11,A22,F_A1,F_A2,F_MIN,F_C,X_MIN1 REAL PUT_IN_INTERVAL EXTERNAL PUT_IN_INTERVAL C C---If AL0 = AL1 (unlikely case. Calling routine should take care of that) C---then take neares point in interval (A1,A2) as a minimiser. It does not C---matter which pont is used IF(AL0.EQ.AL1) THEN X_MIN = PUT_IN_INTERVAL(A1,A2,AL0) RETURN ENDIF C C---Convert derivative to X = (AL-Al0)/(AL1-AL0), AL0 --> 0.0 AL1 --> 1.0 AL1M0 = AL1 - AL0 DF10 = DF0*AL1M0 A11 = (A1-AL0)/AL1M0 A22 = (A2-AL0)/AL1M0 CALL QUAD_INTER(F0,F1,DF10,A,B,C) C C---Calculate values of interpolating polynom for all possible points X0 = PUT_IN_INTERVAL(A11,A22,0.0) X1 = PUT_IN_INTERVAL(A11,A22,1.0) F_A1 = C+A11*(B+A*A22) F_A2 = C+A22*(B+A*A22) F_X0 = C+X0*(B+A*X0) F_X1 = C+X1*(B+A*X1) X_MIN = X0 F_MIN = F_X0 C CALL MIN_OF_2(X1,X_MIN,F_X1,F_MIN) CALL MIN_OF_2(X0,X_MIN,F_X0,F_MIN) CALL MIN_OF_2(X1,X_MIN,F_X1,F_MIN) C C--Find minimum of interpolating polynom in the given interval IF(A.NE.0.0) THEN X_MIN1= -B/(2.0*A) X_MIN1 = PUT_IN_INTERVAL(A11,A22,X_MIN) F_C = A*X_MIN1**2 + B*X_MIN1 + C CALL MIN_OF_2(X_MIN1,X_MIN,F_C,F_MIN) ENDIF X_MIN = AL0 + X_MIN*(AL1-AL0) RETURN END C SUBROUTINE QUAD_INTER(F0,F1,DF0,A,B,C) C c---Given f(0),f'(0) and f(1) calculates coefficients for quadratic C---interpolation for z = (X-X0)/(X1-X0) REAL F0,F1,DF0,A,B,C C A = F1-F0-DF0 B = DF0 C = F0 RETURN END C SUBROUTINE MIN_CUBE_INTER(X0,X1,F0,F1,DF0,DF1,A1,A2,X_MIN) C C---given a, b, f(a), f'(a), f(b), f'(b) and a1, a2 finds minimum of C---cubic interpolating polynom in the interval (min(a1,a2),max(a1,a2)) C---inclusive IMPLICIT NONE REAL X0,X1,F0,F1,DF0,DF1,A1,A2,X_MIN C REAL A,B,C,D,X_MIN1,X_MIN2,X1MX0,DF10,DF11,A11,A22,X_0,X_1, & F_X0,F_X1,F_A1,F_A2,F_MIN,F_M1,F_M2 DOUBLE PRECISION D_LOC REAL EPS_LOC REAL PUT_IN_INTERVAL EXTERNAL PUT_IN_INTERVAL DATA EPS_LOC /1.0E-8/ C IF(X1.EQ.X0) THEN X_MIN = X1 RETURN ENDIF C C---Convert argument to switch X0 --> 0.0, X1 --> 1.0 X1MX0 = X1-X0 DF10 = DF0*X1MX0 DF11 = DF1*X1MX0 A11 = (A1-X0)/X1MX0 A22 = (A2-X0)/X1MX0 C CALL CUBE_INTER(F0,F1,DF10,DF11,A,B,C,D) X_0 = PUT_IN_INTERVAL(A11,A22,0.0) X_1 = PUT_IN_INTERVAL(A11,A22,1.0) F_X0 = D + X_0*(C+X_0*(B+X_0*A)) F_X1 = D + X_1*(C+X_1*(B+X_1*A)) F_A1 = D + A11*(C+A11*(B+A11*A)) F_A2 = D + A22*(C+A22*(B+A22*A)) C C---First find minimum of function in four points X_0,X_1,A11,A22 C---These points migth coincide X_MIN = X_0 F_MIN = F_X0 CALL MIN_OF_2(X_1,X_MIN,F_X1,F_MIN) CALL MIN_OF_2(A11,X_MIN,F_A1,F_MIN) CALL MIN_OF_2(A22,X_MIN,F_A2,F_MIN) C C----Now find minimum of cubic polynom and compare it with previous C----values IF(ABS(A).GT.EPS_LOC) THEN D_LOC = B**2-3.0*A*C IF(D_LOC.GE.0.0) THEN X_MIN1 = (-B+SQRT(D_LOC))/(3*A) X_MIN1 = PUT_IN_INTERVAL(A11,A22,X_MIN1) F_M1 = D + X_MIN1*(C + X_MIN1*(B + X_MIN1*A)) CALL MIN_OF_2(X_MIN1,X_MIN,F_M1,F_MIN) X_MIN2 = (-B-SQRT(D_LOC))/(3*A) X_MIN2 = PUT_IN_INTERVAL(A11,A22,X_MIN2) F_M2 = D + X_MIN2*(C + X_MIN2*(B+ X_MIN2*A)) CALL MIN_OF_2(X_MIN2,X_MIN,F_M2,F_MIN) ENDIF ELSE IF(ABS(B).GT.EPS_LOC) THEN X_MIN1 = -C/(2.0*B) X_MIN1 = PUT_IN_INTERVAL(A11,A22,X_MIN1) F_M1 = D + X_MIN1*(C + X_MIN1*B) CALL MIN_OF_2(X_MIN1,X_MIN,F_M1,F_MIN) ENDIF ENDIF C C---Convert to original scale X_MIN = X0 + X_MIN*X1MX0 cd STOP RETURN END C SUBROUTINE CUBE_INTER(F0,F1,DF0,DF1,A,B,C,D) C C---Given f(0), f'(0), f(1) and f'(1) calculates coefficient of C---cubic interpolation for z = (X-X0)/(X1-X0) IMPLICIT NONE REAL F0,F1,DF0,DF1,A,B,C,D C D = F0 C = DF0 B = 3*(F1-F0) - 2.0*DF0-DF1 A = DF0 + DF1 - 2.0*(F1-F0) RETURN END C REAL FUNCTION PUT_IN_INTERVAL(A1,A2,X) C C--Puts X into interval (A1,A2). if it is not in (A1,A2) then C--takes nearest point. A1 could be less or greater than A2 IMPLICIT NONE REAL A1,A2,X C PUT_IN_INTERVAL = AMIN1(AMAX1(A1,A2),AMAX1(AMIN1(A1,A2),X)) END C SUBROUTINE MIN_OF_2(X_1,X_2,F_1,F_2) C C---It changes values of F_2 and X_2 to F_1 and X_1 if F_2.Gt.F_1 IMPLICIT NONE REAL X_1,X_2,F_1,F_2 C IF(F_2.GT.F_1) THEN X_2 = X_1 F_2 = F_1 ENDIF RETURN END C SUBROUTINE CUB_SPLINE_COEFS(N_POINTS,X_SAMPLE,Y_SAMPLE,COEFS_S) IMPLICIT NONE C C----Coefficients for cubic spline interpolation C----This routine uses not-aknot conditions on the boundary. C----see de Boor, A practical guide to splines (1978) C---- REAL X_SAMPLE(*),Y_SAMPLE(*),COEFS_S(4,*) INTEGER N_POINTS C C----Local variables INTEGER I,N_POINTS1 REAL TEMP,DELTA,DFDX1,DFDX2 C C---Consider N = 1 C---It will have to be considered in interpolating subroutine IF(N_POINTS.LE.1) THEN COEFS_S(1,1) = Y_SAMPLE(1) COEFS_S(2,1) = 0.0 COEFS_S(3,1) = 0.0 COEFS_S(4,1) = 0.0 RETURN ENDIF C C---Points should be distinct and increasing N_POINTS1 = N_POINTS-1 DO I=1,N_POINTS1 IF(X_SAMPLE(I).GE.X_SAMPLE(I+1)) THEN CALL ERRWRT(-1,'Points must be increasing and distinct') CALL ERRWRT(1,'Error in CUB_SPLINE_COEFS') ENDIF ENDDO DO I=1,N_POINTS COEFS_S(1,I) = Y_SAMPLE(I) ENDDO C C---N=2. Linear interpolation IF(N_POINTS.EQ.2) THEN COEFS_S(2,1) = & (Y_SAMPLE(2)-Y_SAMPLE(1))/(X_SAMPLE(2)-X_SAMPLE(1)) COEFS_S(3,1) = 0.0 COEFS_S(4,1) = 0.0 RETURN ENDIF C C--N>=3. general case DO I=2,N_POINTS COEFS_S(3,I) = X_SAMPLE(I) - X_SAMPLE(I-1) COEFS_S(4,I) = (COEFS_S(1,I)-COEFS_S(1,I-1))/COEFS_S(3,I) ENDDO COEFS_S(4,1) = COEFS_S(3,3) COEFS_S(3,1) = COEFS_S(3,2) + COEFS_S(3,3) COEFS_S(2,1) = ((COEFS_S(3,2)+2.0*COEFS_S(3,1))*COEFS_S(4,2)* & COEFS_S(3,3) + COEFS_S(3,2)*COEFS_S(3,2)* & COEFS_S(4,3))/COEFS_S(3,1) DO I=2,N_POINTS1 TEMP = -COEFS_S(3,I+1)/COEFS_S(4,I-1) COEFS_S(2,I) = TEMP*COEFS_S(2,I-1)+3.0*(COEFS_S(3,I)* & COEFS_S(4,I+1)+COEFS_S(3,I+1)*COEFS_S(4,I)) COEFS_S(4,I) = TEMP*COEFS_S(3,I-1)+2.0*(COEFS_S(3,I)+ & COEFS_S(3,I+1)) ENDDO C C--not-a-knot at the last point IF(N_POINTS.GT.3) THEN TEMP = COEFS_S(3,N_POINTS1) + COEFS_S(3,N_POINTS) COEFS_S(2,N_POINTS) = ((COEFS_S(3,N_POINTS) + 2.0*TEMP)* & COEFS_S(4,N_POINTS)*COEFS_S(3,N_POINTS1)+ & COEFS_S(3,N_POINTS)**2*(COEFS_S(1,N_POINTS1)- & COEFS_S(1,N_POINTS-2))/COEFS_S(3,N_POINTS1))/ & TEMP TEMP = -TEMP/COEFS_S(4,N_POINTS1) COEFS_S(4,N_POINTS) = COEFS_S(3,N_POINTS1) ELSE COEFS_S(2,N_POINTS) = 2.0*COEFS_S(4,N_POINTS) COEFS_S(4,N_POINTS) = 1.0 TEMP = -1.0/COEFS_S(4,N_POINTS1) ENDIF COEFS_S(4,N_POINTS) = TEMP*COEFS_S(3,N_POINTS1) + & COEFS_S(4,N_POINTS) COEFS_S(2,N_POINTS) = (COEFS_S(2,N_POINTS1)*TEMP + & COEFS_S(2,N_POINTS))/COEFS_S(4,N_POINTS) C C---Back substitution DO I=N_POINTS1,1,-1 COEFS_S(2,I) = (COEFS_S(2,I)-COEFS_S(3,I)*COEFS_S(2,I+1))/ & COEFS_S(4,I) ENDDO C C---Now generate coefficients for cubic spline (picewise cubic functions) DO I=2,N_POINTS DELTA = COEFS_S(3,I) DFDX1 = (COEFS_S(1,I) - COEFS_S(1,I-1))/DELTA DFDX2 = COEFS_S(2,I-1) + COEFS_S(2,I) - 2.0*DFDX1 COEFS_S(3,I-1) = (DFDX1-COEFS_S(2,I-1)-DFDX2)/DELTA COEFS_S(4,I-1) = (DFDX2/DELTA)/DELTA ENDDO C C---Normal termination RETURN END C SUBROUTINE CUBIC_SPLINE_VALUE1(N_POINTS,X_POINTS,X_CURRENT, & I_INTERVAL,Y_CURRENT,COEFS_S) IMPLICIT NONE C c---Calculates value of function at point X_CURRENT using cubic interpolation c---Interpolation coefficient are assumed to be known. Interval number is known C REAL X_POINTS(*),COEFS_S(4,*),X_CURRENT,Y_CURRENT INTEGER I_INTERVAL,N_POINTS,IERROR C C----Local variables INTEGER I_INT1 REAL DX C c----Calculate value of function at the given point I_INT1 = MAX(1,MIN(N_POINTS-1,I_INTERVAL)) DX = X_CURRENT-X_POINTS(I_INT1) Y_CURRENT = COEFS_S(1,I_INT1) + & DX*(COEFS_S(2,I_INT1)+ & DX*(COEFS_S(3,I_INT1)+ & DX* COEFS_S(4,I_INT1))) RETURN END C SUBROUTINE CUBIC_SPLINE_VALUE2(N_POINTS,X_POINTS,X_CURRENT, & Y_CURRENT,COEFS_S) IMPLICIT NONE C c---Calculates value of function at point X_CURRENT using cubic interpolation c---Interpolation coefficient are assumed to be known. Interval number is not C---known C REAL X_POINTS(*),COEFS_S(4,*),X_CURRENT,Y_CURRENT INTEGER I_INTERVAL,N_POINTS,IERROR C C----Local variables REAL DX INTEGER K1,K2,K C C----Find interval where X_CURRENT belongs. First two conditions C----means that extrapolation will be used. Not a very good idea. IF(X_CURRENT.LE.X_POINTS(1)) THEN K1 = 1 ELSEIF(X_CURRENT.GE.X_POINTS(N_POINTS)) THEN K1= N_POINTS - 1 ELSE K1 = 1 K2 = N_POINTS-1 1 CONTINUE IF(K2.GT.K1+1) THEN K = (K1+K2)/2 IF(X_POINTS(K).GT.X_CURRENT) THEN K2 = K ELSE K1 = K ENDIF GOTO 1 ENDIF ENDIF C c----Calculate value of function at the given point DX = X_CURRENT-X_POINTS(K1) Y_CURRENT = COEFS_S(1,K1) + & DX*(COEFS_S(2,K1)+ & DX*(COEFS_S(3,K1)+ & DX* COEFS_S(4,K1))) RETURN END C SUBROUTINE LINTER_VALUE2(N_POINTS,X_POINTS,Y_POINTS,X_CURRENT, & K1,Y_CURRENT) IMPLICIT NONE C c---Calculates value of function at point X_CURRENT using cubic interpolation c---Interpolation coefficient are assumed to be known. Interval number is not C---known C REAL X_POINTS(*),Y_POINTS(*),X_CURRENT,Y_CURRENT INTEGER I_INTERVAL,N_POINTS,IERROR C C----Local variables REAL DX,A,B INTEGER K1,K2,K C C----Find interval where X_CURRENT belongs. First two conditions C----means that extrapolation will be used. Not a very good idea. cd WRITE(*,*)N_POINTS,K1 IF(N_POINTS.EQ.1) THEN Y_CURRENT = Y_POINTS(1) X_CURRENT = X_POINTS(1) K1 =1 RETURN ENDIF IF(K1.GE.1.AND.K1.LT.N_POINTS) THEN IF(X_CURRENT.GE.X_POINTS(K1).AND.X_CURRENT.LE.X_POINTS(K1+1)) & GOTO 20 ENDIF cd IF(K1.LT.1.OR.K1.GT.N_POINTS-1) THEN IF(X_CURRENT.LE.X_POINTS(1)) THEN K1 = 1 ELSEIF(X_CURRENT.GE.X_POINTS(N_POINTS)) THEN K1= N_POINTS - 1 ELSE K1 = 1 K2 = N_POINTS 1 CONTINUE IF(K2.GT.K1+1) THEN K = (K1+K2)/2 IF(X_POINTS(K).GT.X_CURRENT) THEN K2 = K ELSE K1 = K ENDIF GOTO 1 ENDIF ENDIF cd ENDIF C c----Calculate value of function at the given point 20 CONTINUE B = Y_POINTS(K1) A = (Y_POINTS(K1+1)-Y_POINTS(K1))/(X_POINTS(K1+1)-X_POINTS(K1)) DX = X_CURRENT-X_POINTS(K1) Y_CURRENT = A*DX+B IF(X_CURRENT.LT.X_POINTS(1)) THEN Y_CURRENT = AMAX1(0.1*Y_POINTS(1), & AMIN1(10.0*Y_POINTS(1),Y_CURRENT)) ELSE IF(X_CURRENT.GT.X_POINTS(N_POINTS)) THEN Y_CURRENT = AMAX1(0.1*Y_POINTS(N_POINTS), & AMIN1(10.0*Y_POINTS(N_POINTS),Y_CURRENT)) ENDIF RETURN END