math.f90

Ground state of the time-independent Gross-Pitaevskii equation

!**********************************************************************
!
! Mathematical routines
!
!**********************************************************************

subroutine GaussHermite(n,x,w)
  !
  ! Calculates the abscissas and weight factors for the
  ! Gauss-Hermite quadrature
  !
  implicit none
  integer, parameter :: dp=kind(1.d0), itMax=20
  real(dp), parameter :: pi=3.1415926535897932d0, epsilon=1.d-14

  integer, intent(in) :: n
  real(dp), dimension(n), intent(out) :: x,w

  integer :: i,j,jj
  real(dp) :: pimq,nmh,xold,xnew,Ha,Hb,Hc,Hp,fact
  logical :: isOdd 

  pimq = pi**(-0.25d0)
  fact = sqrt(real(2*n,dp))
  nmh = 1.d0/sqrt(real(n,dp))

  if(2*(n/2) /= n) then
     isOdd = .true.
  else
     isOdd = .false.
  end if

  do j = n/2+1,n
     ! initial guess
     if(j == n/2+1) then
        if(isOdd) then
           xold = 0.d0
        else
        xold = 2.8259d-2 + 0.9778d0*nmh
        end if
     elseif(isOdd .and. j == n/2+2) then
        xold = 2.4866d-2 + 2.0863d0*nmh
     elseif(j == n) then
        xold = x(n-1) + 0.5d0 + 1.2635d0*nmh
     else
        xold = 2.d0*x(j-1)-x(j-2)
     end if

     !*** find zero by Newton's method

     do i=1,itMax
        ! calculate value of the Hermite polynomial
        Ha = 0.d0
        Hb = pimq
        do jj = 1,n
           Hc = xold*sqrt(2./real(jj,dp))*Hb &
                & - sqrt(real(jj-1,dp)/real(jj,dp))*Ha
           Ha = Hb
           Hb = Hc
        end do
        
        Hp = fact*Ha
        xnew = xold - Hc/Hp
        
        if(abs(xnew-xold) <= epsilon) exit
        xold = xnew
     end do

     if(i > itMax) then
        write(*,*) "Abscissa ",j," not found in GaussHermite"
        stop
     end if

     x(j) = xnew
     x(n-j+1) = -xnew
     w(j) = 2.d0/Hp**2
     w(n-j+1) = w(j)

  end do

  return
end subroutine GaussHermite


subroutine harmonic(n,nx,x,psi)
  !
  ! Calculates the nth eigenfunction of the harmonic oscillator
  ! for k=1, m=1, hbar=1
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: pi=3.1415926535897932d0

  integer, intent(in) :: n,nx
  real(dp), dimension(nx), intent(in) :: x
  real(dp), dimension(nx), intent(out) :: psi

  integer :: i
  real(dp) :: norm
  real(dp), external :: He,factorial

  norm = 1.d0/sqrt(2.d0**n*factorial(n)*sqrt(pi))

  do i=1,nx
     psi(i) = He(n,x(i))*exp(-0.5d0*x(i)**2)
  end do

  psi=psi*norm

return
end subroutine harmonic


real(kind(1.d0)) function He(n,x)
  !
  ! Calculates the Hermite polynomial H_n(x) by recursion
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n
  real(dp), intent(in) :: x

  integer :: i
  real(dp) :: a,b,c

  ! Verify input value
  if(n < 0) then
     write(*,*) 'Negative value n in function He(n,x), n=',n
     stop
  end if

  a = 1.d0
  b = 2.d0*x

  if(n==0) then
     c = a
  elseif(n==1) then
     c = b
  else
     do i=1,n-1
        c = 2.d0*(x*b-real(i,dp)*a)
        a = b
        b = c
     end do
  end if

  He = c

  return
end function He


real(kind(1.d0)) function factorial(n)
  !
  ! Calculates the factorial of integer n
  !
  ! Maximum value of n is 170
  !
  implicit none
  integer, parameter :: dp=kind(1.d0), nmax=170

  integer, intent(in) :: n

  integer :: j
  integer, save :: ncalc=0
  real(dp), dimension(0:nmax), save :: a(0:nmax)

  a(0) = 1.d0

  if (n < 0 .or. n > nmax) then
     write(*,*) 'Incorrect value for function factorial, n=',n 
     stop 
  elseif(n > ncalc) then
     do j = ncalc+1,n
        a(j) = real(j,dp)*a(j-1)
     end do
     ncalc = n
  endif
  factorial = a(n)         

  return
end function factorial


real(kind(1.d0)) function integrate3D(nx,ny,nz,dx,dy,dz,f)
  !
  ! Three-dimensional integral
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: nx,ny,nz
  real(dp), intent(in) :: dx,dy,dz
  real(dp), dimension(nx,ny,nz), intent(in) :: f

  integer :: k
  real(dp), dimension(nz) :: integrand

  real(dp), external :: integrate1D,integrate2D

  do k = 1,nz
     integrand(k) = integrate2D(nx,ny,dx,dy,f(1,1,k))
  end do

  integrate3D = integrate1D(nz,dz,integrand)

  return
end function integrate3D


real(kind(1.d0)) function integrate2D(nx,ny,dx,dy,f)
  !
  ! Two-dimensional integral
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: nx,ny
  real(dp), intent(in) :: dx,dy
  real(dp), dimension(nx,ny), intent(in) :: f

  integer :: j
  real(dp), dimension(ny) :: integrand

  real(dp), external :: integrate1D

  do j = 1,ny
     integrand(j) = integrate1D(nx,dx,f(1,j))
  end do

  integrate2D = integrate1D(ny,dy,integrand)

  return
end function integrate2D


real(kind(1.d0)) function integrate1D(nx,dx,f)
  !
  ! One-dimensional integral (trapezoidal rule)
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: nx
  real(dp), intent(in) :: dx
  real(dp), dimension(nx), intent(in) :: f

  integrate1D = (0.5d0*(f(1) + f(nx)) + sum(f(2:nx-1)))*dx

  return
end function integrate1D