PROGRAM EINSTEIN IMPLICIT NONE REAL*8 :: Z, E, J1, J2, J1A, J2A, TIME, J1G, J2G INTEGER :: I, J, K, L, EX_P INTEGER :: ZINT, EINT Z = 0.0 E = 0.0 EX_P = - 8.0 ZINT = 999 EINT = 7 OPEN(1, FILE='test.dat', STATUS='UNKNOWN') 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 ! DO L = 1, 10 Z = 0.0 EX_P = - 8.0 DO I = 1, EINT DO J = 1, 9 E = J * 10.0**EX_P DO K = 1, ZINT CALL EINSTEIN_G(Z, E, J1, J2) ! CALL EINSTEIN_A(Z, E, J1A, J2A) CALL EINSTEIN_GUO (Z, E, J1G, J2G) ! CALL EINSTEIN_NEW (Z, E, J1G, J2G) IF (J1 > TINY(0.1) .OR. J1A > TINY(0.1) .OR. J1A > TINY(0.1) ) THEN ! WRITE (1,12) -LOG10(E), Z, (J1-J1A), (J2-J2A) ! WRITE (1, '(4F15.4)') -LOG10(E), Z, (J1G-J1)/J1*100, (J2G-J2)/J2*100 ! WRITE (1, 12) -LOG10(E), Z, (J1G-J1)/J1*100, (J2G-J2)/J2*100 ! WRITE (1,12) -LOG10(E), Z, LOG(J1A), -LOG(ABS(J2A)) ! WRITE (1,12) -LOG10(E), Z, LOG(J1), -LOG(ABS(J2)) ! WRITE (1,12) -LOG10(E), Z, (J1-J1A)/J1*100, (J2-J2A)/J2*100 ! WRITE (1,12) -LOG10(E), Z, (J1-J1A)/J1*100, (J2-J2A)/J2*100 WRITE (1,12) -LOG10(E), Z, (J1-J1G)/J1*100, (J2-J2G)/J2*100 ELSE WRITE (1,12) -LOG10(E), Z, -999.0, -999.0 END IF Z = Z + 0.01 END DO Z = 0.01 END DO EX_P = EX_P + 1 END DO ! WRITE (*,*) L ! END DO CALL CPU_TIME(TIME) WRITE (*,*) TIME 12 FORMAT(F12.8,',', F15.3,',', F20.4,',', F30.4) END PROGRAM !*********************************************************************| SUBROUTINE EINSTEIN_NEW (Z1, R, J1, J2) 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_NEW !*********************************************************************| SUBROUTINE EINSTEIN_G (Z1, R, J1, J2) IMPLICIT NONE REAL*8, PARAMETER :: GAMMA= 0.5772156649 REAL*8, PARAMETER :: PI= 3.14159265358979 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 END IF J1 = I1 J2 = I2 END SUBROUTINE EINSTEIN_G !*********************************************************************| SUBROUTINE F_1 (Z, E, F1) ! O.K 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 IMPLICIT NONE INTEGER K REAL*8 Z, E, F1, F2 CALL F_1(Z-K,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) 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 REAL*8, PARAMETER :: gamma = 0.5772156649 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) IMPLICIT NONE REAL*8, PARAMETER :: GAMMA= 0.5772156649 REAL*8, PARAMETER :: PI= 3.14159265358979 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*(D1-E1)) 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) !*********************************************************************| END IF J1 = I1 J2 = I2 !*********************************************************************| END SUBROUTINE EINSTEIN_A