MODULE globals
  DOUBLE PRECISION, PARAMETER :: dCloseEnoughToZero = 1D-12
  DOUBLE PRECISION, PARAMETER :: normSquareCloseToZero = 1D-1
  DOUBLE PRECISION, PARAMETER :: ulinks = 0
  DOUBLE PRECISION, PARAMETER :: urechts = 10D0
  INTEGER :: nLimit = 4000
  DOUBLE PRECISION, PARAMETER :: c = 2.99792458D8
  DOUBLE PRECISION, PARAMETER :: hbarc = 197.327053D-9
  DOUBLE PRECISION, PARAMETER :: hbar = hbarc/c
  DOUBLE PRECISION, PARAMETER :: g = 9.80665
  DOUBLE PRECISION, PARAMETER :: mcsq = 939.56563D6
  DOUBLE PRECISION, PARAMETER :: m = mcsq/(c*c)
  DOUBLE PRECISION :: h ! = 0.1D0 ! automatisch
  DOUBLE PRECISION, PARAMETER :: E0 = (m*g**2*hbar**2/2D0)**(1D0/3D0)
  CONTAINS
  INTEGER FUNCTION getSpanIndex(z) RESULT(result)
    DOUBLE PRECISION z
    result = (z - ulinks)/h
  END FUNCTION
END MODULE globals
PROGRAM neutron
 USE globals
 IMPLICIT NONE
 CALL clearQCache()
 h = (urechts - ulinks)/DBLE(nLimit)
 CALL findAllResults(0D0, 30D0, 0.1D0)
END
SUBROUTINE findAllResults(aepsillinks, epsilrechts, roh0)
 USE globals
 IMPLICIT NONE
 INTEGER i
 INTEGER j
 DOUBLE PRECISION integrateUsquareBySimpson
 DOUBLE PRECISION, intent(in) :: aepsillinks, epsilrechts, roh0
 DOUBLE PRECISION findLowestEpsilon
 DOUBLE PRECISION epsil, epsillinks
 DOUBLE PRECISION epsils(0:3)
 DOUBLE PRECISION u
 DOUBLE PRECISION findEpsilonWithUAbove
 DOUBLE PRECISION r
 LOGICAL allOrthogonal
 epsillinks = aepsillinks
 OPEN(11, FILE="L", STATUS="REPLACE")
 DO i = 0, 3, 1
  epsil = findLowestEpsilon(epsillinks, epsilrechts, roh0)
  IF (DABS(u(0, epsil, roh0)) .lt. dCloseEnoughToZero) THEN
    epsils(i) = epsil
    WRITE(*,*) "Eigenenergie: ", epsil*E0*1D12, " peV"      ! <-- (un)comment to use
    CALL writeNormU(epsil, roh0)                             ! <-- (un)comment to use
  ELSE
    !WRITE(*,*) "???", epsil, DABS(u(0, epsil, roh0))
    EXIT
  END IF
  ! u(0, epsillinks) > 0 sicherstellen
  epsillinks = findEpsilonWithUAbove(epsil, epsilrechts, roh0)
  IF (epsilrechts .LT. epsillinks) THEN
    EXIT
  END IF
  !h = h/2
  !nLimit = nLimit*2
  !WRITE(*,*) "===="
 END DO
 CLOSE(11)
 WRITE(*,*) "(siehe Datei 'L' für die normierten Wellenfunktionen)"
 allOrthogonal = .TRUE.
 DO i = 0, 3, 1
   DO j = 0, 3, 1
     r = integrateUsquareBySimpson(ulinks, urechts, 0.001D0, epsils(i), epsils(j), roh0)
     r = r/DSQRT(integrateUsquareBySimpson(ulinks, urechts, 0.001D0, epsils(i), epsils(i)))
     r = r/DSQRT(integrateUsquareBySimpson(ulinks, urechts, 0.001D0, epsils(j), epsils(j)))
     IF (i .EQ. j) THEN
       IF (r .LT. normSquareCloseToZero) THEN
         allOrthogonal = .FALSE.
       END IF
     ELSE ! should be orthogonal
       IF (r .GE. normSquareCloseToZero) THEN
         allOrthogonal = .FALSE.
       END IF
     END IF
   END DO
 END DO
 IF(.NOT. allOrthogonal) THEN
   WRITE(*,*) "keine orthonormierten Ergebnisse"
 END IF
END
! Voraussetzung: u(0,alinks)>0, u(0,arechts)<0
DOUBLE PRECISION FUNCTION findLowestEpsilon(alinks, arechts, roh0) RESULT (prevEpsil)
 USE globals
 IMPLICIT NONE
 INTEGER i
 DOUBLE PRECISION u
 DOUBLE PRECISION findEpsilon
 DOUBLE PRECISION roh0
 DOUBLE PRECISION, intent(in) :: alinks, arechts
 DOUBLE PRECISION :: epsillinks, epsilrechts
 DOUBLE PRECISION :: epsil
 DOUBLE PRECISION findEpsilonWithUBelow
 epsillinks = alinks
 epsilrechts = arechts
 epsil = 0D0
 DO i = 0, 5000, 1
   prevEpsil = epsil
   epsil = findEpsilon(epsillinks, epsilrechts, roh0)
   IF (DABS(u(0, epsil, roh0)) .ge. dCloseEnoughToZero) THEN ! wir nehmen an, es existiert gar keine passende Wellenfunktion mehr.
     EXIT
   END IF
   ! u(0, epsilrechts) < 0 sicherstellen
   epsilrechts = findEpsilonWithUBelow(epsillinks, epsil, roh0)
   IF (epsilrechts .LE. epsillinks) THEN
     EXIT
   END IF
 END DO
END
DOUBLE PRECISION FUNCTION findEpsilonWithUBelow(epsillinks, aepsil, roh0) RESULT(result) ! <--
  DOUBLE PRECISION, intent(in) :: epsillinks, aepsil
  result = aepsil
  result = result - 0.00001D0
  DO WHILE (u(0, result, roh0) .GE. 0D0 .AND. result .GE. epsillinks)
    result = result - 0.00001D0
  END DO
END
DOUBLE PRECISION FUNCTION findEpsilonWithUAbove(aepsil, epsilrechts, roh0) RESULT(result) ! -->
  DOUBLE PRECISION, intent(in) :: aepsil, epsilrechts
  result = aepsil
  result = result + 0.00001D0
  DO WHILE (u(0, result, roh0) .LE. 0D0 .AND. result .LE. epsilrechts)
    result = result + 0.00001D0
  END DO
END

! Voraussetzung: u(0,alinks)>0, u(0,arechts)<0
DOUBLE PRECISION FUNCTION findEpsilon(alinks, arechts, roh0) RESULT (epsil)
 USE globals
 IMPLICIT NONE 
 DOUBLE PRECISION roh0
 DOUBLE PRECISION u
 DOUBLE PRECISION, intent(in) :: alinks, arechts
 DOUBLE PRECISION :: epsillinks, epsilrechts
 DOUBLE PRECISION :: prevLinks, prevRechts
 INTEGER i
 epsillinks = alinks
 epsilrechts = arechts
 DO i = 0, 10000, 1
  prevLinks = epsillinks
  prevRechts = epsilrechts
  epsil = (epsilrechts + epsillinks)/2D0
  !WRITE(*,*) "findEpsilon is at epsil", epsil, epsillinks, epsilrechts
  IF (DABS(u(0, epsil, roh0)) .lt. dCloseEnoughToZero) THEN
   !WRITE(*,*) "OK", epsil, DABS(u(0, epsil, roh0))
   exit
  ELSE IF (u(0, epsillinks, roh0)*u(0, epsil, roh0) .LT. 0D0) THEN
   epsilrechts = epsil
  ELSE 
   epsillinks = epsil
  END IF
  IF(epsillinks .EQ. prevLinks .AND. epsilrechts .EQ. prevRechts) THEN
    EXIT
  END IF
  !WRITE (*,*) epsillinks, epsilrechts
 END DO
 IF(i .GT. 10000) THEN
   !WRITE(*,*) "(aufgegeben)", epsil, DABS(u(0, epsil, roh0))
 END IF
END

!w(n+1) = epsil - (roh(n+1) - roh0)
!Q(n+1) = (1D0 + h*h/12D0*w(n+1))*u(n+1)
!Q(n+1) = 12D0*u(n) - 10D0*Q(n) - Q(n-1)
!Q(n) = 12D0*u(n+1) - 10D0*Q(n+1) - Q(n+2)

DOUBLE PRECISION FUNCTION w(n, epsil, roh0) RESULT(result)
 USE globals
 IMPLICIT NONE
 INTEGER n
 DOUBLE PRECISION epsil, roh0
 DOUBLE PRECISION roh
 result = epsil - (roh(n, roh0) - roh0)
END

DOUBLE PRECISION FUNCTION roh(n, roh0) RESULT(result)
 USE globals
 IMPLICIT NONE
 INTEGER n
 DOUBLE PRECISION roh0
 result = roh0 + n*h
END

DOUBLE PRECISION RECURSIVE FUNCTION Qvonu(n, epsil, roh0) RESULT(result)
 USE globals
 IMPLICIT NONE
 INTEGER n
 DOUBLE PRECISION epsil, roh0
 DOUBLE PRECISION u, w  
 result = (1D0 + h*h/12D0*w(n, epsil, roh0))*u(n, epsil, roh0)
END
SUBROUTINE clearQCache()
 IMPLICIT NONE
 COMMON previousResults, previousResultKnowns, previousEpsils
 DOUBLE PRECISION previousResults(0:9999)
 LOGICAL previousResultKnowns(0:9999)
 DOUBLE PRECISION previousEpsils(0:9999)
 INTEGER i
 DO i=0,9999,1
   previousResultKnowns(i) = .FALSE.
   previousResults(i) = 42D0
   previousEpsils(i) = 0D0
 END DO
END
DOUBLE PRECISION RECURSIVE FUNCTION Q(n, epsil, roh0) RESULT(result)
 USE globals
 IMPLICIT NONE
 INTEGER n
 DOUBLE PRECISION epsil, roh0
 DOUBLE PRECISION u, Qvonu
 COMMON previousResults, previousResultKnowns, previousEpsils
 DOUBLE PRECISION previousResults(0:9999)
 LOGICAL previousResultKnowns(0:9999)
 DOUBLE PRECISION previousEpsils(0:9999)
 IF (previousResultKnowns(n) .AND. previousEpsils(n) .EQ. epsil) THEN
   result = previousResults(n)
 ELSE
   previousEpsils(n) = epsil
   IF (n .eq. nLimit) THEN
      result = Qvonu(n, epsil, roh0)
   ELSE IF (n .eq. nLimit - 1) THEN
      result = Qvonu(n, epsil, roh0)
   ELSE 
      result = 12D0*u(n+1, epsil, roh0) - 10D0*Q(n+1, epsil, roh0) - Q(n+2, epsil, roh0)
   END IF
   previousResults(n) = result
   previousResultKnowns(n) = .TRUE.
 END IF
END
DOUBLE PRECISION RECURSIVE FUNCTION u(n, epsil, roh0) RESULT(result)
 USE globals
 IMPLICIT NONE
 INTEGER, intent(in) :: n
 DOUBLE PRECISION, intent(in) :: epsil, roh0
 DOUBLE PRECISION Q, w
 IF (n .eq. nLimit) THEN
  result = 0D0
 ELSE IF (N .eq. nLimit - 1) THEN
  result = 1D-10
 ELSE 
  result = Q(n, epsil, roh0)/(1D0 + h*h/12D0*w(n, epsil, roh0))
 END IF
END

DOUBLE PRECISION FUNCTION integrateUsquareBySimpson(left, right, dz, aepsil, bepsil, roh0) RESULT(result)
  USE globals
  DOUBLE PRECISION, intent(in) :: left, right, dz, aepsil, bepsil, roh0
  DOUBLE PRECISION u
  INTEGER count
  DOUBLE PRECISION factor
  DOUBLE PRECISION z
  DOUBLE PRECISION ua(0:99999)
  INTEGER i
  count = (right - left)/dz
  DO i = getSpanIndex(left), getSpanIndex(right), 1
    ua(i) = u(i, aepsil, roh0)
  END DO
  ! int f(x) dx = dz/3 * (f(x_0) + 4 f(x_1) + 2 f(x_2) + 4 f(x_3) + ... + 4 f(x_n-1) + f(x_n))
  result = ua(getSpanIndex(left))*u(getSpanIndex(left), bepsil, roh0)
  DO i = 1,count - 1,1
    IF (MOD(i, 2) .EQ. 0) THEN
      factor = 2D0
    ELSE
      factor = 4D0
    END IF
    z = left + dz*i
    result = result + factor*(ua(getSpanIndex(z))*u(getSpanIndex(z), bepsil, roh0))
  END DO
  result = result + ua(getSpanIndex(right))*u(getSpanIndex(right), bepsil, roh0)
  result = result*dz/3
END

SUBROUTINE writeNormU(epsil, roh0)
  USE globals
  IMPLICIT NONE
  DOUBLE PRECISION roh
  DOUBLE PRECISION u
  DOUBLE PRECISION epsil, roh0
  DOUBLE PRECISION length
  DOUBLE PRECISION integrateUsquareBySimpson
  INTEGER i
  length = DSQRT(integrateUsquareBySimpson(ulinks, urechts, 0.001D0, epsil, epsil, roh0))
  DO i = 0, nLimit - 1, 1
    WRITE (11,*) roh(i, roh0), u(i, epsil, roh0)/length
  END DO
  WRITE (11,*) ""
  WRITE (11,*) ""
END

