SUBROUTINE ACT_RHM(TK, XIL, XGK, XPY, S, T, U, ILMENITE, 1 GEIKIELITE, HEMATITE, PYROPHANITE, LSTOP) C C Program to calculate the activities of ilmenite, geikielite, C Hematite and pyrophanite according to the solution model of C Ghiorso and Sack (1990, Appendix) - extention of Ghiorso (1990) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DOUBLE PRECISION ILMENITE, ILMENITE_PURE LOGICAL LSTOP, CONVERGED, DEBUG EXTERNAL ORDER DIMENSION WORK(27), AUXVAR(29) C P = 1.0 HIL = 17000.0 SIL = 0.0 HGK = 17000.0 SGK = 0.0 HPY = 17000.0 SPY = 0.0 DWILHM = -19000.0 DWGKHM = -19000.0 DWPYHM = -19000.0 WILHM = 5000.0 WGKHM = 5000.0 WPYHM = 5000.0 DWILGK = 12000.0 DWGKIL = 8000.0 DWILPY = 17000.0 DWPYIL = 17000.0 DWGKPY = 8000.0 DWPYGK = 12000.0 WILGKDIS = 25000.0 WILGKORD = 5000.0 WILPYDIS = 2200.0 WILPYORD = 2200.0 WGKPYDIS = 25000.0 WGKPYORD = 5000.0 C AUXVAR( 1) = TK AUXVAR( 2) = P AUXVAR( 3) = XIL AUXVAR( 4) = XGK AUXVAR( 5) = XPY AUXVAR( 6) = HIL AUXVAR( 7) = SIL AUXVAR( 8) = HGK AUXVAR( 9) = SGK AUXVAR(10) = HPY AUXVAR(11) = SPY AUXVAR(12) = DWILHM AUXVAR(13) = DWGKHM AUXVAR(14) = DWPYHM AUXVAR(15) = WILHM AUXVAR(16) = WGKHM AUXVAR(17) = WPYHM AUXVAR(18) = DWILGK AUXVAR(19) = DWGKIL AUXVAR(20) = DWILPY AUXVAR(21) = DWPYIL AUXVAR(22) = DWGKPY AUXVAR(23) = DWPYGK AUXVAR(24) = WILGKDIS AUXVAR(25) = WILGKORD AUXVAR(26) = WILPYDIS AUXVAR(27) = WILPYORD AUXVAR(28) = WGKPYDIS AUXVAR(29) = WGKPYORD C DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY C WORK(1) = XIL*0.9 ! Initial guess for s WORK(2) = XGK*0.9 ! Initial guess for t WORK(3) = XPY*0.9 ! Initial guess for u ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) S = WORK(1) T = WORK(2) U = WORK(3) IF(MODE.NE.0) THEN WRITE(*,'(1X,''Mode (VARMET), s, t, u = '',I3, 1X, 3G12.6)') 2 MODE, S, T, U IF((S.LT.-XIL).OR.(S.GT.XIL).OR. 1 (T.LT.-XGK).OR.(T.GT.XGK).OR. 2 (U.LT.-XPY).OR.(U.GT.XPY)) THEN WRITE(*,*) 'Abort from ACT_RHM.' LSTOP = .TRUE. RETURN END IF END IF C RTLNHM = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - 4.0*DWPYHM*U**2*XPY - 4.0*DGPY*U**2*XPY 3 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 4 - 4.0*DWILPY*S**2*XPY - 4.0*DWILHM*S**2*XPY + WILHM*XIL**2 5 + WILHM*XGK*XIL - WILGKDIS*XGK*XIL + WGKHM*XGK*XIL 6 - 4.0*DWPYIL*U**2*XIL - 4.0*DWPYHM*U**2*XIL 7 - 4.0*DWGKIL*T**2*XIL - 4.0*DWGKHM*T**2*XIL 8 - 4.0*DWILHM*S**2*XIL - 4.0*DGIL*S**2*XIL + WGKHM*XGK**2 9 - 4.0*DWPYHM*U**2*XGK - 4.0*DWPYGK*U**2*XGK A - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK B - 4.0*DWILHM*S**2*XGK - 4.0*DWILGK*S**2*XGK C - S*U*WILPYORD + S*U*WILPYDIS - S*T*WILGKORD + S*T*WILGKDIS D - T*U*WGKPYORD + T*U*WGKPYDIS + 2.0*DWPYHM*U**2 E + 3.0*DGPY*U**2 + DWPYGK*T*U + DWGKPY*T*U + DWPYIL*S*U F + DWILPY*S*U + 2.0*DWGKHM*T**2 + 3.0*DGGK*T**2 + DWILGK*S*T G + DWGKIL*S*T + 2.0*DWILHM*S**2 + 3.0*DGIL*S**2 HEMATITE = (1.0-XIL-XGK-XPY)**2*EXP(RTLNHM/(8.3143*TK)) RTLNIL = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - WPYHM*XPY + WILPYDIS*XPY - WILHM*XPY 3 - 4.0*DWPYHM*U**2*XPY - 4.0*DGPY*U**2*XPY 4 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 5 - 4*DWILPY*S**2*XPY - 4.0*DWILHM*S**2*XPY 6 + 4.0*DWILPY*S*XPY + 4.0*DWILHM*S*XPY + WILHM*XIL**2 7 + WILHM*XGK*XIL - WILGKDIS*XGK*XIL + WGKHM*XGK*XIL 8 - 2.0*WILHM*XIL - 4.0*DWPYIL*U**2*XIL - 4.0*DWPYHM*U**2*XIL 9 - 4.0*DWGKIL*T**2*XIL - 4.0*DWGKHM*T**2*XIL A - 4.0*DWILHM*S**2*XIL - 4.0*DGIL*S**2*XIL B + 4.0*DWILHM*S*XIL + 4.0*DGIL*S*XIL + WGKHM*XGK**2 C - WILHM*XGK + WILGKDIS*XGK - WGKHM*XGK D - 4.0*DWPYHM*U**2*XGK - 4.0*DWPYGK*U**2*XGK E - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK F - 4.0*DWILHM*S**2*XGK - 4.0*DWILGK*S**2*XGK G + 4.0*DWILHM*S*XGK + 4.0*DWILGK*S*XGK - S*U*WILPYORD H + U*WILPYORD + S*U*WILPYDIS - U*WILPYDIS + WILHM I - S*T*WILGKORD + T*WILGKORD + S*T*WILGKDIS -T*WILGKDIS J - T*U*WGKPYORD + T*U*WGKPYDIS + 2.0*DWPYIL*U**2 K + 4.0*DWPYHM*U**2 + 3.0*DGPY*U**2 + DWPYGK*T*U + DWGKPY*T*U L + DWPYIL*S*U + DWILPY*S*U - DWPYIL*U - DWILPY*U M + 2.0*DWGKIL*T**2 + 4.0*DWGKHM*T**2 + 3.0*DGGK*T**2 N + DWILGK*S*T + DWGKIL*S*T - DWILGK*T - DWGKIL*T O + 4.0*DWILHM*S**2 + 5.0*DGIL*S**2 - 4.0*DWILHM*S P - 6.0*DGIL*S + DGIL ILMENITE = 0.25*(XIL+S)*(XIL+S+XGK+T+XPY+U) 1 * EXP(RTLNIL/(8.3143*TK)) RTLNGK = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - WPYHM*XPY + WGKPYDIS*XPY - WGKHM*XPY 3 - 4.0*DWPYHM*U**2*XPY - 4.0*DGPY*U**2*XPY 4 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 5 + 4.0*DWGKPY*T*XPY + 4.0*DWGKHM*T*XPY - 4.0*DWILPY*S**2*XPY 6 - 4.0*DWILHM*S**2*XPY + WILHM*XIL**2 + WILHM*XGK*XIL 7 - WILGKDIS*XGK*XIL + WGKHM*XGK*XIL - WILHM*XIL 8 + WILGKDIS*XIL - WGKHM*XIL - 4.0*DWPYIL*U**2*XIL 9 - 4.0*DWPYHM*U**2*XIL - 4.0*DWGKIL*T**2*XIL A - 4.0*DWGKHM*T**2*XIL + 4.0*DWGKIL*T*XIL + 4.0*DWGKHM*T*XIL B - 4.0*DWILHM*S**2*XIL - 4.0*DGIL*S**2*XIL + WGKHM*XGK**2 C - 2.0*WGKHM*XGK - 4.0*DWPYHM*U**2*XGK - 4.0*DWPYGK*U**2*XGK D - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK E + 4.0*DWGKHM*T*XGK + 4.0*DGGK*T*XGK - 4.0*DWILHM*S**2*XGK F - 4.0*DWILGK*S**2*XGK - S*U*WILPYORD + S*U*WILPYDIS G - S*T*WILGKORD + S*WILGKORD + S*T*WILGKDIS - S*WILGKDIS H - T*U*WGKPYORD + U*WGKPYORD + T*U*WGKPYDIS - U*WGKPYDIS I + WGKHM + 4.0*DWPYHM*U**2 + 2.0*DWPYGK*U**2 + 3.0*DGPY*U**2 J + DWPYGK*T*U + DWGKPY*T*U + DWPYIL*S*U + DWILPY*S*U K - DWPYGK*U - DWGKPY*U + 4.0*DWGKHM*T**2 + 5.0*DGGK*T**2 L + DWILGK*S*T + DWGKIL*S*T - 4.0*DWGKHM*T - 6.0*DGGK*T M + 4.0*DWILHM*S**2 + 2.0*DWILGK*S**2 + 3.0*DGIL*S**2 N - DWILGK*S - DWGKIL*S + DGGK GEIKIELITE = 0.25*(XGK+T)*(XIL+S+XGK+T+XPY+U) 1 * EXP(RTLNGK/(8.3143*TK)) C RTLNPY = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - 2.0*WPYHM*XPY - 4.0*DWPYHM*U**2*XPY 3 - 4.0*DGPY*U**2*XPY + 4.0*DWPYHM*U*XPY + 4.0*DGPY*U*XPY 4 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 5 - 4.0*DWILPY*S**2*XPY - 4.0*DWILHM*S**2*XPY 6 + WILHM*XIL**2 + WILHM*XGK*XIL - WILGKDIS*XGK*XIL 7 + WGKHM*XGK*XIL - WPYHM*XIL + WILPYDIS*XIL - WILHM*XIL 8 - 4.0*DWPYIL*U**2*XIL - 4.0*DWPYHM*U**2*XIL 9 + 4.0*DWPYIL*U*XIL + 4.0*DWPYHM*U*XIL - 4.0*DWGKIL*T**2*XIL A - 4.0*DWGKHM*T**2*XIL - 4.0*DWILHM*S**2*XIL B - 4.0*DGIL*S**2*XIL + WGKHM*XGK**2 - WPYHM*XGK C + WGKPYDIS*XGK - WGKHM*XGK - 4.0*DWPYHM*U**2*XGK D - 4.0*DWPYGK*U**2*XGK + 4.0*DWPYHM*U*XGK + 4.0*DWPYGK*U*XGK E - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK F - 4.0*DWILHM*S**2*XGK - 4.0*DWILGK*S**2*XGK + WPYHM G - S*U*WILPYORD + S*WILPYORD + S*U*WILPYDIS - S*WILPYDIS H - S*T*WILGKORD + S*T*WILGKDIS - T*U*WGKPYORD + T*WGKPYORD I + T*U*WGKPYDIS - T*WGKPYDIS + 4.0*DWPYHM*U**2 J + 5.0*DGPY*U**2 + DWPYGK*T*U + DWGKPY*T*U + DWPYIL*S*U K + DWILPY*S*U - 4.0*DWPYHM*U - 6.0*DGPY*U + 2.0*DWGKPY*T**2 L + 4.0*DWGKHM*T**2 + 3.0*DGGK*T**2 + DWILGK*S*T + DWGKIL*S*T M - DWPYGK*T - DWGKPY*T+2*DWILPY*S**2 + 4.0*DWILHM*S**2 N + 3.0*DGIL*S**2 - DWPYIL*S - DWILPY*S + DGPY PYROPHANITE = 0.25*(XPY+U)*(XIL+S+XGK+T+XPY+U) 1 * EXP(RTLNPY/(8.3143*TK)) C C Calculations for pure ilmenite C AUXVAR( 3) = 1.0 AUXVAR( 4) = 0.0 AUXVAR( 5) = 0.0 WORK(1) = 0.9 ! Initial guess for s WORK(2) = 0.0 WORK(3) = 0.0 ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) S_PURE = WORK(1) IF(MODE.NE.0) WRITE(*,*) 'Mode (VARMET) = ', MODE RTLNIL_PURE = DGIL*(1.0-S_PURE)**2 ILMENITE_PURE = 0.25*(1.0+S_PURE)**2 1 * EXP(RTLNIL_PURE/(8.3143*TK)) ILMENITE = ILMENITE/ILMENITE_PURE C C Calculations for pure geikielite C AUXVAR( 3) = 0.0 AUXVAR( 4) = 1.0 AUXVAR( 5) = 0.0 WORK(1) = 0.0 WORK(2) = 0.9 ! Initial guess for t WORK(3) = 0.0 ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) T_PURE = WORK(2) IF(MODE.NE.0) WRITE(*,*) 'Mode (VARMET) = ', MODE RTLNGK_PURE = DGGK*(1.0-T_PURE)**2 GEIKIELITE_PURE = 0.25*(1.0+T_PURE)**2 1 * EXP(RTLNGK_PURE/(8.3143*TK)) GEIKIELITE = GEIKIELITE/GEIKIELITE_PURE C C Calculations for pure pyrophanite C AUXVAR( 3) = 0.0 AUXVAR( 4) = 0.0 AUXVAR( 5) = 1.0 WORK(1) = 0.0 WORK(2) = 0.0 WORK(3) = 0.9 ! Initial guess for u ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) U_PURE = WORK(3) IF(MODE.NE.0) WRITE(*,*) 'Mode (VARMET) = ', MODE RTLNPY_PURE = DGPY*(1.0-U_PURE)**2 PYROPHANITE_PURE = 0.25*(1.0+U_PURE)**2 1 * EXP(RTLNPY_PURE/(8.3143*TK)) PYROPHANITE = PYROPHANITE/PYROPHANITE_PURE C LSTOP = .FALSE. RETURN END SUBROUTINE ORDER(B, P0, G, AUXVAR, LSTOP) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION AUXVAR(29), B(3), G(3) LOGICAL LSTOP C TK = AUXVAR( 1) P = AUXVAR( 2) XIL = AUXVAR( 3) XGK = AUXVAR( 4) XPY = AUXVAR( 5) HIL = AUXVAR( 6) SIL = AUXVAR( 7) HGK = AUXVAR( 8) SGK = AUXVAR( 9) HPY = AUXVAR(10) SPY = AUXVAR(11) DWILHM = AUXVAR(12) DWGKHM = AUXVAR(13) DWPYHM = AUXVAR(14) WILHM = AUXVAR(15) WGKHM = AUXVAR(16) WPYHM = AUXVAR(17) DWILGK = AUXVAR(18) DWGKIL = AUXVAR(19) DWILPY = AUXVAR(20) DWPYIL = AUXVAR(21) DWGKPY = AUXVAR(22) DWPYGK = AUXVAR(23) WILGKDIS = AUXVAR(24) WILGKORD = AUXVAR(25) WILPYDIS = AUXVAR(26) WILPYORD = AUXVAR(27) WGKPYDIS = AUXVAR(28) WGKPYORD = AUXVAR(29) C S = B(1) T = B(2) U = B(3) XHM = 1.0 - XIL - XGK - XPY IF((S.LT.0.0).OR.(S.GT.XIL).OR.(T.LT.0.0).OR.(T.GT.XGK).OR. 1 (U.LT.0.0).OR.(U.GT.XPY)) THEN LSTOP = .TRUE. RETURN END IF C IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WILPYDIS-WILHM)*XIL*XPY 1 + (-WPYHM+WGKPYDIS-WGKHM)*XGK*XPY + (WPYHM+DGPY)*XPY 2 + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 3 + (2.0*DWGKPY+2.0*DWGKHM)*T**2*XPY 4 + (2.0*DWILPY+2.0*DWILHM)*S**2*XPY - WILHM*XIL**2 5 + (-WILHM+WILGKDIS-WGKHM)*XGK*XIL + (WILHM+DGIL)*XIL 6 + (2.0*DWPYIL+2.0*DWPYHM)*U**2*XIL 7 + (2.0*DWGKIL+2.0*DWGKHM)*T**2*XIL 8 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL - WGKHM*XGK**2 9 + (WGKHM+DGGK)*XGK A + (2.0*DWPYHM+2.0*DWPYGK)*U**2*XGK B + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK C + (2.0*DWILHM+2.0*DWILGK)*S**2*XGK D + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) E + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) F + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) G + (-2.0*DWPYHM-3.0*DGPY)*U**2 H + (-2.0*DWGKHM-3.0*DGGK)*T**2 I + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILPY+2.0*DWILHM)*S*XPY 1 + 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 2 + 2.0*(2.0*DWILHM+2.0*DWILGK)*S*XGK 3 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 4 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 5 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 2.0*(2.0*DWGKPY+2.0*DWGKHM)*T*XPY 1 + 2.0*(2.0*DWGKIL+2.0*DWGKHM)*T*XIL 2 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 3 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 4 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 5 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYIL+2.0*DWPYHM)*U*XIL 2 + 2.0*(2.0*DWPYHM+2.0*DWPYGK)*U*XGK 3 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 4 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 5 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T+XPY-U) 1 *LOG((XIL-S+XGK-T+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T+XPY+U) 1 *LOG((XIL+S+XGK+T+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = WILPYDIS*XIL*XPY + WGKPYDIS*XGK*XPY + DGPY*XPY 1 + 2.0*DGPY*U**2*XPY + 2.0*DWGKPY*T**2*XPY 2 + 2.0*DWILPY*S**2*XPY + WILGKDIS*XGK*XIL + DGIL*XIL 3 + 2.0*DWPYIL*U**2*XIL + 2.0*DWGKIL*T**2*XIL 4 + 2.0*DGIL*S**2*XIL + DGGK*XGK + 2.0*DWPYGK*U**2*XGK 5 + 2.0*DGGK*T**2*XGK + 2.0*DWILGK*S**2*XGK 6 + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 7 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 8 + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 9 - 3.0*DGPY*U**2 - 3.0*DGGK*T**2 - 3.0*DGIL*S**2 G(1) = 4.0*DWILPY*S*XPY + 4.0*DGIL*S*XIL + 4.0*DWILGK*S*XGK 1 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 2 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 - 6.0*DGIL*S G(2) = 4.0*DWGKPY*T*XPY + 4.0*DWGKIL*T*XIL + 4.0*DGGK*T*XGK 1 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 2 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 - 6.0*DGGK*T G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYIL*U*XIL + 4.0*DWPYGK*U*XGK 1 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 2 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 4 - 6.0*DGPY*U C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T+XPY-U) 1 *LOG((XIL-S+XGK-T+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T+XPY+U) 1 *LOG((XIL+S+XGK+T+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WGKPYDIS-WGKHM)*XGK*XPY 1 + (WPYHM+DGPY)*XPY + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (2.0*DWGKPY+2.0*DWGKHM)*T**2*XPY - WGKHM*XGK**2 3 + (WGKHM+DGGK)*XGK + (2.0*DWPYHM+2.0*DWPYGK)*U**2*XGK 4 + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 5 + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 6 + (-2.0*DWPYHM-3.0*DGPY)*U**2 7 + (-2.0*DWGKHM-3.0*DGGK)*T**2 G(1) = 0.0 G(2) = 2.0*(2.0*DWGKPY+2.0*DWGKHM)*T*XPY 1 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 2 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYHM+2.0*DWPYGK)*U*XGK 2 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK-T+XPY-U)*LOG((XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK+T+XPY+U)*LOG((XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WILPYDIS-WILHM)*XIL*XPY 1 + (WPYHM+DGPY)*XPY + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (2.0*DWILPY+2.0*DWILHM)*S**2*XPY - WILHM*XIL**2 3 + (WILHM+DGIL)*XIL + (2.0*DWPYIL+2.0*DWPYHM)*U**2*XIL 4 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL 6 + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 7 + (-2.0*DWPYHM-3.0*DGPY)*U**2 8 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILPY+2.0*DWILHM)*S*XPY 1 + 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 2 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 0.0 G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYIL+2.0*DWPYHM)*U*XIL 2 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XPY-U)*LOG((XIL-S+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XPY+U)*LOG((XIL+S+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK P0 = - WILHM*XIL**2 + (-WILHM+WILGKDIS-WGKHM)*XGK*XIL 1 + (WILHM+DGIL)*XIL + (2.0*DWGKIL+2.0*DWGKHM)*T**2*XIL 2 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL - WGKHM*XGK**2 3 + (WGKHM+DGGK)*XGK + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 4 + (2.0*DWILHM+2.0*DWILGK)*S**2*XGK 5 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 6 + (-2.0*DWGKHM-3.0*DGGK)*T**2 7 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 1 + 2.0*(2.0*DWILHM+2.0*DWILGK)*S*XGK 2 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 2.0*(2.0*DWGKIL+2.0*DWGKHM)*T*XIL 1 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 2 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T)*LOG((XIL-S+XGK-T)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T)*LOG((XIL+S+XGK+T)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = WGKPYDIS*XGK*XPY + DGPY*XPY + 2.0*DGPY*U**2*XPY 1 + 2.0*DWGKPY*T**2*XPY + DGGK*XGK + 2.0*DWPYGK*U**2*XGK 2 + 2.0*DGGK*T**2*XGK + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 - 3.0*DGPY*U**2 - 3.0*DGGK*T**2 G(1) = 0.0 G(2) = 4.0*DWGKPY*T*XPY + 4.0*DGGK*T*XGK - 6.0*DGGK*T 1 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYGK*U*XGK - 6.0*DGPY*U 1 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK-T+XPY-U)*LOG((XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK+T+XPY+U)*LOG((XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGPY = HPY - TK*SPY P0 = WILPYDIS*XIL*XPY + DGPY*XPY + 2.0*DGPY*U**2*XPY 1 + 2.0*DWILPY*S**2*XPY + DGIL*XIL + 2.0*DWPYIL*U**2*XIL 2 + 2.0*DGIL*S**2*XIL + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 - 3.0*DGPY*U**2 - 3.0*DGIL*S**2 G(1) = 4.0*DWILPY*S*XPY + 4.0*DGIL*S*XIL - 6.0*DGIL*S 1 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) G(2) = 0.0 G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYIL*U*XIL - 6.0*DGPY*U 1 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XPY-U)*LOG((XIL-S+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XPY+U)*LOG((XIL+S+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK P0 = WILGKDIS*XGK*XIL + DGIL*XIL + 2.0*DWGKIL*T**2*XIL 1 + 2.0*DGIL*S**2*XIL + DGGK*XGK + 2.0*DGGK*T**2*XGK 2 + 2.0*DWILGK*S**2*XGK - 3.0*DGGK*T**2 - 3.0*DGIL*S**2 3 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(1) = 4.0*DGIL*S*XIL + 4.0*DWILGK*S*XGK - 6.0*DGIL*S 1 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(2) = 4.0*DWGKIL*T*XIL + 4.0*DGGK*T*XGK - 6.0*DGGK*T 1 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(3) = 0.0 C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T)*LOG((XIL-S+XGK-T)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T)*LOG((XIL+S+XGK+T)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL P0 = - WILHM*XIL**2 + (WILHM+DGIL)*XIL 1 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL 2 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 1 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 0.0 G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*LOG((XIL-S)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGGK = HGK - TK*SGK P0 = - WGKHM*XGK**2 + (WGKHM+DGGK)*XGK 1 + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 2 + (-2.0*DWGKHM-3.0*DGGK)*T**2 G(1) = 0.0 G(2) = 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 1 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*LOG((XGK-T)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (WPYHM+DGPY)*XPY 1 + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (-2.0*DWPYHM-3.0*DGPY)*U**2 G(1) = 0.0 G(2) = 0.0 G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*LOG((XPY-U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN P0 = 0.0 G(1) = 0.0 G(2) = 0.0 G(3) = 0.0 C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL P0 = DGIL*(1.0-S**2) G(1) = - 2.0*DGIL*S G(2) = 0.0 G(3) = 0.0 C P0 = P0 + 8.3143*TK*(1.0+S)*LOG((1.0+S)/2.0) G(1) = G(1) + 8.3143*TK*LOG((1.0+S)/2.0) P0 = P0 + 8.3143*TK*(1.0-S)*LOG((1.0-S)/2.0) G(1) = G(1) - 8.3143*TK*LOG((1.0-S)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGGK = HGK - TK*SGK P0 = DGGK*(1.0-T**2) G(1) = 0.0 G(2) = - 2.0*DGGK*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*(1.0+T)*LOG((1.0+T)/2.0) G(2) = G(2) + 8.3143*TK*LOG((1.0+T)/2.0) P0 = P0 + 8.3143*TK*(1.0-T)*LOG((1.0-T)/2.0) G(2) = G(2) - 8.3143*TK*LOG((1.0-T)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGPY = HPY - TK*SPY P0 = DGPY*(1.0-U**2) G(1) = 0.0 G(2) = 0.0 G(3) = - 2.0*DGPY*U C P0 = P0 + 8.3143*TK*(1.0+U)*LOG((1.0+U)/2.0) G(3) = G(3) + 8.3143*TK*LOG((1.0+U)/2.0) P0 = P0 + 8.3143*TK*(1.0-U)*LOG((1.0-U)/2.0) G(3) = G(3) - 8.3143*TK*LOG((1.0-U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN STOP 'Absurd logic error in call to min(s,t,u) function.' C END IF C P0 = P0/1000.0 G(1) = G(1)/1000.0 G(2) = G(2)/1000.0 G(3) = G(3)/1000.0 C LSTOP = .FALSE. RETURN END SUBROUTINE VARMET(SUB, N, B, P0, G, X, C, T, BMAT, NB, AUXVAR, 1 ITMAX, DEBUG, MODE) C C Algorithm 21. Variable metric method. (page 159-160) C Nash, J.C. Compact Numerical Methods for Computers: Linear C Algebra and Function Minimisation. John Wiley and Sons, C New York, 227 pages (1979) C C SUB(B, P0, G, AUXVAR, LSTOP) SUBROUTINE which returns the value C of the function (P0) and its gradient (G) at the "point" B C subject to values of auxillary "known" parameters contained C in the array AUXVAR. LSTOP is returned with a value of .TRUE. C if the function or its gradient cannot be evaluated at B. C N Number of parameters in the function C B Parameter values C ON INPUT: Initial guess C ON OUTPUT: Final values of the function parameters which C minimize the function defined in S C P0 Function value at the minimum C G Gradient of the function evaluated at the minimum C X Working array of length N C C Working array of length N C T Working array of length N C BMAT Woorking square matrix of size N*N C NB Row dimension of BMAT C AUXVAR User supplied auxillary variables passed to S() in order C to evaluate the function at its gradient C ITMAX Maximum number of allowed calls to S() C DEBUG Logical variable, set to .TRUE. in DEBUG mode C MODE Status variable, On exit Mode equals C 0 when computation is successful C 1 when initial point is infeasible C 2 when iteration limit is exceeded C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION AUXVAR(*), B(N), G(N), X(N), C(N), T(N), BMAT(NB, N) INTEGER COUNT DOUBLE PRECISION K LOGICAL LSTOP, DEBUG, CONVERGED, RESTART C C Step 0 Enter n, the number of parameters in the function C Enter b[i], i=1,2,...,n, the starting values for the C parameters C Let ig = 0 (graient evaluation counter) C Let ifn = 0 (function evaluation counter) C Let w = 0.2, let tolerance = 0.0001 to define acceptable C point strategy C IG = 0 IFN = 0 W = 0.2 TOLERANCE = 0.0001 MODE = 0 CONVERGED = .FALSE. RESTART = .TRUE. C C Step 1 Compute P0 = S(b). If this is not possible, stop, as the C initial point is infeasible. Let ifn = ifn + 1 C Step 2 Compute g(b), the gradient at the current point. C Let ig = ig + 1 C CALL SUB(B, P0, G, AUXVAR, LSTOP) IF(LSTOP) THEN MODE = 1 RETURN END IF IFN = IFN + 1 IG = IG + 1 DO WHILE(.NOT.CONVERGED) C C Step 3 Set BMAT to the unit matrix 1n. C Let ilast = ig to mark point at which B most recently C set to unity C IF(RESTART) THEN DO I=1,N DO J=1,N BMAT(I,J) = 0.0 END DO BMAT(I,I) = 1.0 END DO ILAST = IG RESTART = .FALSE. END IF C C The top of the main iteration loop. To follow the progress of the C minimisation ig, ifn, and P0 could be printed here. C IF(DEBUG) THEN WRITE(*,'(1X,''P0, IG, IFN:'', G12.6,2I10)') P0, IG, IFN END IF C DO I=1,N X(I) = B(I) C(I) = G(I) END DO C C Step 5 Form t = -Bg and t^T g C D1 = 0.0 DO I=1,N S = 0.0 DO J=1,N S = S - BMAT(I,J)*G(J) END DO T(I) = S D1 = D1 - S*G(I) END DO C C Step 6 If D1>0, goto step 7 to perform linear search. Otherwise, C search direction is uphill, so if ilast=ig, goto step 18, C since the algorithm is taken to have converged. Otherwise, C goto step 3 to restart iteration. C IF(D1.GT.0.0) THEN C C Step 7 Let k=1, to intialise step length C K = 1.0 C C Step 8 Let count = 0 C LSTOP = .TRUE. DO WHILE(LSTOP) COUNT = 0 DO I=1,N B(I) = X(I) + K*T(I) IF(B(I).EQ.X(I)) COUNT = COUNT + 1 END DO C C Step 9 Convergence test C If Count < n, goto step 10 to evaluate the function C Otherwise, if ilast=ig, goto step 18, since the C algorithm is taken to have converged. Otherwise goto C step 3 to restart iteration C IF(COUNT.LT.N) THEN C C Step 10 Compute P=S(b), let ifn=ifn+1. If the function cannot be C computed, let k = w*k, then goto step 8 to repeat the C adjustment of parameters C CALL SUB(B, P, G, AUXVAR, LSTOP) IFN = IFN + 1 IF(IFN.GT.ITMAX) THEN MODE = 2 RETURN END IF IF(LSTOP) THEN K = W*K ELSE C C Step 11 Acceptance test on new point C If P < P0 - D1*k*tolerance, goto 12 to continue iteration. C Otherwise, let k=w*k, then goto step 8 C IF(P.LT.(P0-D1*K*TOLERANCE)) THEN C C Step 12 Let p0=P to save the function value. C Compute g at b, the gradient at the new point. C Let ig = ig + 1 C P0 = P IG = IG + 1 C C Step 13 Prepare for matrix update. C D1 = 0.0 DO I=1,N T(I) = K*T(I) C(I) = G(I) - C(I) D1 = D1 + T(I)*C(I) END DO C C Step 14 If D1 <= 0, goto step 3 to restart iteration since positive C definiteness of matrix cannot be guaranteed and there is C a danger of a division by zero C IF(D1.LE.0.0) THEN RESTART = .TRUE. ELSE C C Step 15 Computation of BMATy, y^T BMAT y C D2 = 0.0 DO I=1,N S = 0.0 DO J=1,N S = S + BMAT(I,J)*C(J) END DO X(I) = S D2 = D2 + S*C(I) END DO C C Step 16 Let D2 = 1 + D2/D1. Note this is (15.36) times d1. C D2 = 1.0 + D2/D1 DO I=1,N DO J=1,N BMAT(I,J) = BMAT(I,J) - (T(I)*X(J) 1 + X(I)*T(J) 2 - D2*T(I)*T(J)) / D1 END DO END DO C C This completes the update. C C Note that there are other ways to organise the update. The C rounding errors involved in some cases may cause different C schemes to have very different convergence patterns. So far no C one scheme seems to be better than the others, so the most direct C implementation has been chosen here. C C Step 17 Goto step 4 for another iteration. C END IF ELSE K = W*K LSTOP = .TRUE. END IF END IF ELSE IF(IG.EQ.ILAST) THEN CONVERGED = .TRUE. LSTOP = .FALSE. ELSE RESTART = .TRUE. LSTOP = .FALSE. END IF END DO ELSE IF(IG.EQ.ILAST) THEN CONVERGED = .TRUE. ELSE RESTART = .TRUE. END IF END DO C C Step 18 PRINT IG, IFN, P0, b[i], i=1,2,...,n. C STOP. C IF(DEBUG) THEN WRITE(*,'(1X,''IG, IFN, P0:'',2I10,G12.6)') IG, IFN, P0 END IF C RETURN END