PROGRAM EINSTEIN
	! This program shows how to calculate the Einstain integrals. The application 
	!	of its subroutine can be found in the example, missouri.f90.
	! To test this program, simply RUN it and provide the values of z and E.
	!
	! 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. 132(?)
	!
	! Written by JUNKE GUO, 07/05/2005

		implicit none
		real :: z, E, J10, J1, J20, F2, S, J2, E1
		real, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651
		integer :: k, n
 
		write(*,*) 'Input the values of z and E:'
		read(*,*) z, E

		! 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(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
		! for integer z
		else
			! for z=0
			if (anint(z) == 0) then
				J1 = 1 - E
				J2 = -1 + E - E*log(E)
			else
				! for z=1				
				J1 = -1 - log(E) + E
				J2 = 1 - E + E*log(E) - 0.5*log(E)**2
			end if

			! for z=2,3,4,...
			if (anint(z) > 1) then 
				E1 = E/(1-E)
				do k = 2, anint(z)
					J1 = 1.0/(k-1.0)*E/E1**k - k/(k-1.0)*J1
					J2 = E/E1**k*log(E)/(k-1.0) - k/(k-1.0)*J2 + J1/(k-1.0)
				end do
			end if
		end if

		write(*,*) 'z=',z, 'E=',E
		write(*,*) 'J1=',J1, 'J2=', J2
				
	CONTAINS
		REAL FUNCTION F1(z,E)
			implicit none
			real, intent(in) :: z, E
			real :: S1, E1
			integer :: k

			E1 = E/(1.0-E)
			S1 = -1.0/(1.0-z)*E1**(1.0-z)
			do k = 2, ceiling(z) + 4
				S1 = S1 + (-1)**k/(k-z)*E1**(k-z)
			end do
			
			F1 = E/E1**z - z*S1
		END FUNCTION F1

	END PROGRAM EINSTEIN