!********************************************************************* ! !********************************************************************* PROGRAM EINSTEIN ! ! This subroutine is accoding to the references: ! A. Roland, P. Mewis and U. Zanke (2004).Discussion on "Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! ! Remarks: Take care of floating point precision truncation of your compiler when you turn on optmization flags ... ! To aggressive optimization makes certain compiler do not execute the benchmark loop as the compilers that the arguments do not change ... ! ! Written by ARON ROLAND, 26/12/2004, Vrdnik, Yugoslawija ! IMPLICIT NONE REAL*8 :: Z, E, J1, J2, J1A, J2A, TIME1, TIME2, J1G, J2G, J1N ,J2N REAL*8 :: LOGERR1, LOGERR2 INTEGER :: I, J, K, L, EX_P INTEGER :: ZINT, EINT , LAUF, ITEST INTEGER, SAVE :: APPROX LOGICAL :: TECPLOT = .FALSE. LOGICAL :: LRAW = .TRUE. LOGICAL :: LTEC = .FALSE. LOGICAL :: LBENCH = .FALSE. LOGICAL :: LCALC = .TRUE. OPEN(1, FILE='result.dat', STATUS='UNKNOWN') IF (TECPLOT) THEN WRITE(1,*) 'TITLE = "EINSTEIN.DAT"' WRITE(1,*) 'Variables = "-LOG(A/H)" "Z" "LOG(I1)" "LOG(-I2)"' WRITE(1,*) 'ZONE T=','"1"', ' I =',9*EINT*ZINT END IF WRITE (*,*) 'Calculation (0) or Bechmarking (1) ?' READ (*,*) ITEST IF (ITEST == 0) THEN LBENCH = .FALSE. LCALC = .TRUE. ELSE LBENCH = .TRUE. LCALC = .FALSE. END IF WRITE (*,*) 'Tecplot (0) or Rawout (1) ?' READ (*,*) ITEST IF (ITEST == 0) THEN LRAW = .FALSE. LTEC = .TRUE. ELSE LRAW = .TRUE. LTEC = .FALSE. END IF IF (LBENCH) THEN WRITE (*,*) 'Which formula you wan´t to use ? Type integer given in brackets' WRITE (*,*) 'GUO RECCURENCE FORMULA (1), 1996' WRITE (*,*) 'GUO & JULIEN recent best approx. (2), 2005' WRITE (*,*) 'ROLAND et Al. explicit approx. (3), 2005' WRITE (*,*) 'GUO RECCURENCE + ROLAND et Al. (4), 2005' WRITE (*,*) 'GUO CLOSURE TO DISCUSSION (5), 2006' READ (*,*) APPROX ELSE WRITE (*,*) 'Which formula you wan´t to use ? Type integer given in brackets' WRITE (*,*) 'GUO RECCURENCE FORMULA (1), 1996' WRITE (*,*) 'ROLAND et Al. explicit approx. (2), 2005' WRITE (*,*) 'GUO RECCURENCE + ROLAND et Al. (3), 2005' WRITE (*,*) 'GUO CLOSURE TO DISCUSSION (4), 2006' READ (*,*) APPROX END IF CALL CPU_TIME(TIME1) IF (LCALC) THEN Z = 0.0 EX_P = - 8.0 ZINT = 999 EINT = 7 ! DO I = 1, EINT ! WRITE (*,*) I, EINT ! DO J = 1, 9 ! E = J * 10.0**EX_P E = 0.001 DO K = 1, ZINT CALL EINSTEIN_G(Z, E, J1, J2) SELECT CASE (APPROX) CASE (1) CALL EINSTEIN_GUO (Z, E, J1N, J2N) CASE (2) CALL EINSTEIN_A (Z, E, J1N, J2N) CASE (3) CALL EINSTEIN_A_G (Z, E, J1N, J2N) CASE (4) CALL EINSTEIN_GUO3 (Z, E, J1N, J2N) CASE DEFAULT STOP END SELECT IF (J1 > TINY(1.0) .AND. J1N > TINY(1.0) .AND. J1N > TINY(1.0) ) THEN IF (LRAW) THEN IF ((J1-J1N)/J1*100 .GT. 0.0) THEN LOGERR1 = LOG( (J1-J1N)/J1*100 )) ELSE LOGERR1 = -LOG( ABS( (J1-J1N)/J1*100 ) ) END IF IF ((J2-J2N)/J2*100 .GT. 0.0) THEN LOGERR2 = LOG( (J2-J2N)/J2*100 ) ELSE LOGERR2 = -LOG( ABS( (J2-J2N)/J2*100 ) ) END IF ! WRITE (1, '(7F15.8)') Z, (J1-J1N)/J1*100, (J2-J2N)/J2*100, LOG(J1), LOG(J1N), -LOG(ABS(J2)), -LOG(ABS(J2N)) WRITE (1, '(3F15.8)') Z, LOGERR1, LOGERR2 END IF IF (LTEC) THEN WRITE (1,12) -LOG10(E), Z, (J1-J1N)/J1*100, (J2-J2N)/J2*100 END IF ELSE WRITE (1, '(8F15.4)') Z, 0., 0. END IF Z = Z + 0.01 END DO Z = 0.0 ! END DO EX_P = EX_P + 1 ! END DO END IF IF (LBENCH) THEN DO L = 1, 10 Z = 0.0 EX_P = - 8.0 ZINT = 999 EINT = 7 DO I = 1, EINT DO J = 1, 9 E = J * 10.0**EX_P ! E = 0.1 DO K = 1, ZINT SELECT CASE (INT(APPROX)) CASE (1) CALL EINSTEIN_GUO (Z, E, J1N, J2N) CASE (2) CALL EINSTEIN_G (Z, E, J1N, J2N) CASE (3) CALL EINSTEIN_A (Z, E, J1N, J2N) CASE (4) CALL EINSTEIN_A_G (Z, E, J1N, J2N) CASE (5) CALL EINSTEIN_GUO3 (Z, E, J1N, J2N) CASE DEFAULT STOP END SELECT Z = Z + 0.01 END DO Z = 0.01 END DO EX_P = EX_P + 1 END DO WRITE (*,*) L END DO END IF CALL CPU_TIME(TIME2) WRITE (*,*) 'THE CALCULATION TIME IS', TIME2-TIME1, '[sec]' 12 FORMAT(F12.8,',', F15.3,',', F20.4,',', F30.4) 13 FORMAT(2F15.4, 8F25.5) END PROGRAM !********************************************************************* ! !********************************************************************* ! SUBROUTINE EINSTEIN_A_G (Z1, R, J1, J2) ! ! This subroutine is accoding to the references: ! A. Roland, P. Mewis and U. Zanke (2004).Discussion on "Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! ! Written by ARON ROLAND, 2/04/2005, Darmstadt, Germany ! IMPLICIT NONE REAL*8, PARAMETER :: GAMMA= 0.5772156649 REAL*8, PARAMETER :: PI= 3.14159265358979 REAL*8 Z, Z1, R, J1, J2 Z = (NINT(Z1 * 100.0))/100.0 IF (ABS(Z-INT(Z)) > 1.E-3) THEN ! Z is not an integer ! IF (R < 10**-3 .OR. Z < 3.0) THEN CALL EINSTEIN_A(Z, R, J1, J2) ELSE CALL EINSTEIN_GUO (Z, R, J1, J2) END IF END IF END SUBROUTINE EINSTEIN_A_G !********************************************************************* ! !********************************************************************* ! ! This subroutine is an extension to the optimal approximation which reduces the error for the Integral J2 ! A. Roland, P. Mewis and U. Zanke (2004).Discussion on "Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! ! Written by ARON ROLAND, 2/04/2005, Darmstadt, Germany ! SUBROUTINE EINSTEIN_A_G (Z1, R, J1OUT, J2OUT) IMPLICIT NONE REAL*8, PARAMETER :: GAMMA= 0.5772156649 REAL*8, PARAMETER :: PI= 3.14159265358979 REAL*8 Z, Z1, R, J1, J2, J1OUT, J2OUT Z = (NINT(Z1 * 100.0))/100.0 IF (ABS(Z-INT(Z)) > 1.E-3) THEN ! Z is not an integer ! IF (R < 10**-3 .OR. Z < 3.0) THEN CALL EINSTEIN_A(Z, R, J1, J2) ELSE CALL EINSTEIN_GUO (Z, R, J1, J2) END IF END IF END SUBROUTINE EINSTEIN_A_G !********************************************************************* ! !********************************************************************* SUBROUTINE EINSTEIN_G (Z1, R, J1, J2) ! ! This subroutine is accoding to the references: ! Guo, J. and Julien, P.Y. (2004)."Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! ! Written by ARON ROLAND, on the basis of JUNKE GUO'S MATHLAB CODE ! IMPLICIT NONE REAL*8, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651 REAL*8 R, Z, EFAK, I1, I2, I2_1, A, B, C, D, E, J1, J2, & & F_G, L, S, R1, R2, R3, R4, Z1, Z2, Z3, Z4, Z1_1, Z1_2, & & A1, B1, C1, D1, S1, E1, I1_1, S2, A2, B2, Z1_3, E2, D2, & & I1_N, PHI, F1, F2 L = LOG(R) Z = (NINT(Z1 * 100.0))/100.0 IF (ABS(Z-INT(Z)) > 1.E-3) THEN ! Z is not an integer ! CALL F_1(Z, R, F1) ! WRITE (*,*) 'F1', F1 I1 = Z*PI/SIN(Z*PI) - F1 PHI = PI**2./6*Z/(1+Z)**0.7152 CALL F_2(Z, R, F2) ! WRITE (*,*) 'F2', F2 I2 = Z*PI/SIN(Z*PI)*(PI*COTAN(Z*PI)-1-1./Z+PHI)-F2 ELSE I1 = 0.0 I2 = 0.0 END IF J1 = I1 J2 = I2 END SUBROUTINE EINSTEIN_G !*********************************************************************| SUBROUTINE F_1 (Z, E, F1) ! O.K ! This subroutine is accoding to the references: ! Guo, J. and Julien, P.Y. (2004)."Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 IMPLICIT NONE INTEGER K REAL*8 Z, E, F1 F1 = ((1-E)**Z)/(E**(Z-1)) DO K = 1, 50 F1 = F1+Z *(-1)**K/(Z-K)*(E/(1-E))**(K-Z) END DO END SUBROUTINE F_1 !*********************************************************************| SUBROUTINE F_2 (Z, E, F2) ! O.K ! This subroutine is accoding to the references: ! Guo, J. and Julien, P.Y. (2004)."Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 IMPLICIT NONE INTEGER K REAL*8 Z, E, F1, F2 CALL F_1(Z,E,F1) F2 = F1 *(LOG(E)+1./(Z-1)) DO K = 1, 50 CALL F_1(Z-K,E,F1) F2 = F2 + Z*(-1)**K*F1/(Z-K)/(Z-K-1) END DO END SUBROUTINE F_2 !********************************************************************* ! !********************************************************************* SUBROUTINE EINSTEIN_GUO (Z, R, J_1, J_2) ! ! This subroutine is accoding to the references: ! Guo, J, Wu, B, and Molinas, A.: ! Analytical Integration of Einstein’s integrals Sediment transport Integrals ! International Conference on Reservoir Sedimentation ! Fort Collins, CO., Sep.9-13, 1996 ! ! Written by ARON ROLAND, 03/03/2005 IMPLICIT NONE REAL*8 I_1, I_2, J_1, J_2, J_10, J_20, E, Ro, Z_1, Z1, R, R1, R2, R3, R4 REAL*8 F1, F2, F1b, phi, Ro1, J_1_1, F3, EFAK, A1, Z, B1, J_2_1, L REAL*8, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651 REAL*8, DIMENSION (5) :: J1_S INTEGER n, m, k, I L = LOG(R) Z = (NINT(Z * 100.0))/100.0 IF ((Z - INT(Z)) >= 0.001) THEN IF (Z .LT. 1.0 ) THEN F3 = (1-gamma)-log(2-Z)+1/(1-Z)+1/(2*(2-Z))+1./(24*(2-Z)**2) J_1 = (Z*pi)/(sin(Z*pi))-R**(1-Z)/(1-Z) J_2 = -Z*pi*F3/sin(Z*pi)-(R**(1-Z)*L)/(1-Z)+ R**(1-Z)/(1-Z)**2 END IF IF (Z .GT. 1.0 ) THEN Z1 = Z-INT(Z) J_1 = (Z1*pi)/(sin(Z1*pi))-R**(1-Z1)/(1-Z1) J1_S(1) = J_1 DO I = 1, INT(Z) Z_1 = I*1.0+(Z-INT(Z)) A1 = 1/(Z_1-1) J_1_1 = A1*((1-R)**Z_1)/(R**(Z_1-1))-(Z_1*A1)*J_1 J_1 = J_1_1 J1_S(I+1) = J_1 END DO F3 = (1-gamma)-log(2-Z1)+1./(1-Z1)+1./(2*(2-Z1))+1./(24*(2-Z1)**2) J_2 = -Z1*pi*F3/sin(Z1*pi)-(R**(1-Z1)*L)/(1-Z1)+R**(1-Z1)/(1-Z1)**2 DO I = 1, INT(Z) Z_1 = I*1.0+(Z-INT(Z)) A1 = 1./(Z_1-1) B1 = R**(Z_1-1)/(1-R)**Z_1 J_2_1 = A1*(L/B1-Z_1*J_2+J1_S(I+1)) J_2 = J_2_1 END DO END IF END IF END !********************************************************************* ! !********************************************************************* SUBROUTINE EINSTEIN_A (Z1, R, J1, J2) ! ! This subroutine is accoding to the references: ! ! This subroutine is accoding to the references: ! A. Roland, P. Mewis and U. Zanke (2004).Discussion on "Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! ! Written by ARON ROLAND, 03/03/2005 IMPLICIT NONE REAL*8, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651 REAL*8 R, Z, EFAK, I1, I2, I2_1, A, B, C, D, E, J1, J2, & & F_G, L, S, R1, R2, R3, R4, Z1, Z2, Z3, Z4, Z1_1, Z1_2, & & A1, B1, C1, D1, S1, E1, I1_1, S2, A2, B2, Z1_3, E2, D2, & & I1_N, PHI, F1, F2 L = LOG(R) Z = (NINT(Z1 * 100.0))/100.0 IF ((Z - INT(Z)) >= 0.001) THEN Z1_1 = Z-1 Z1_2 = Z-2 Z1_3 = Z-3 Z2 = 2-Z Z3 = 3-Z Z4 = 4-Z F_G = (1-GAMMA)-LOG(ABS(Z4))+1/Z3+1./(2*Z4)+1./(12*Z4**2) A = 1./Z1_1 A1 = 1./Z1_2 A2 = 1./Z1_3 B = R**Z1_1/((1-R)**Z) B1 = R**Z1_2/((1-R)**Z1_1) B2 = R**Z1_3/((1-R)**Z1_2) C = Z/(Z1_1) C1 = Z1_1/Z1_2 S1 = SIN(Z1_2*PI) S2 = SIN(Z1_3*PI) D1 = Z1_2*PI/S1 D2 = Z1_3*PI/S2 E1 = R**Z3/Z3 E2 = R**Z4/Z4 I1 = A/B-C*(A1/B1-Z1_1*A1*(A2/B2-(Z1_2*A2)*(D2-E2))) I1_1 =(A1/B1-C1*(D1-E1)) I2_1 = -(Z1_2*PI*F_G)/S1-R**Z3*L/Z3+R**Z3/Z3**2 I2 = A*(L/B-Z*(A1*(L/B1-Z1_1*I2_1+I1_1))+I1) J1 = I1 J2 = I2 ELSE J1 = 0.0 J2 = 0.0 END IF END SUBROUTINE EINSTEIN_A !********************************************************************* ! !********************************************************************* SUBROUTINE EINSTEIN_GUO3(z, E, J1, J2) ! ! This subroutine is accoding to the references: ! Guo, J. and Julien, P.Y. (2004)."Efficient algorithm for computing ! Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201 ! Guo, J. and Julien, P.Y. (2006). Closure of the discussion, JHE, 132(?) ! ! Written by JUNKE GUO, 07/05/2005 ! implicit none REAL*8 :: z, E REAL*8 :: J1, J2 REAL*8 :: J10, J20, F1, F2, S, E1 REAL*8, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651 integer :: k, n ! for non-integer z if (abs(z-anint(z)) > 0.005) then J10 = z*pi/sin(z*pi) ! J1 from zero to 1 J1 = J10 - F1(z,E) ! J1 from E to 1 ! J2 from 0 to 1 n = ceiling(real(z)) + 2 J20 = 1-gamma - log(n-z) - 0.5/(n-z) + 1.0/12.0/(n-z)**2 - 1.0/120.0/(n-z)**4 do k = 1, n J20 = J20 + 1.0/(k-z) end do J20 = -J10*J20 ! J2 from E to 1 F2 = F1(z,E)*(log(E) + 1.0/(z-1.0)) S = -1.0/(z-1.0)/(z-2.0)*F1(z-1.0,E) do k = 2, n S = S + (-1)**k/(z-k)/(z-k-1.0)*F1(z-k,E) end do F2 = F2 + z*S J2 = J20 - F2 else J1 = 0.0 J2 = 0.0 end if END SUBROUTINE REAL*8 FUNCTION F1(z,E) implicit none REAL*8 :: z, E REAL*8 :: S1, E1 integer :: k E1 = E/(1.0-E) S1 = -1.0/(1.0-z)*E1**(1.0-z) do k = 2, ceiling(real(z)) + 4 S1 = S1 + (-1)**k/(k-z)*E1**(k-z) end do F1 = E/E1**z - z*S1 END FUNCTION F1 !********************************************************************* ! !*********************************************************************