!---------------------------------------------------------
! Input: the real+imaginary parts of a wave function
!
! Output: Expectation values
!---------------------------------------------------------
subroutine     observables(psi,V,Nx,xmin,dx, avX,avP,avX2,avP2,avE,DX2,DP2)
 implicit none
 integer                  :: Nx
 real   (8)               :: dx,xmin
 complex(8), dimension(Nx):: psi
 real   (8), dimension(Nx):: V
 complex(8), allocatable  :: obs(:),x(:)
 complex(8)               :: integral
 real   (8)               :: avX,avX2,avP,avP2,avV,avE,DX2,DP2,norm
 integer                  :: i
 complex(8), parameter    :: IMU=(0.0D0,1.0D0)
!---------------------------------------------------------
 ALLOCATE(obs(Nx),x(Nx))
!---------------------------------------------------------
!Define values of x:
 do i = 1, Nx
  x(i) = xmin + (i-1)*dx
 end do
!---------------------------------------------------------
!Norm of psi:
 obs  = psi * DCONJG(psi)
 call integrate(obs,dx,Nx,integral)
 norm = DBLE(integral)
 psi  = psi / norm
!---------------------------------------------------------
!<x>
 obs = psi * x * DCONJG(psi)
 call integrate(obs,dx,Nx,integral)
 avX = DBLE(integral)
!---------------------------------------------------------
!<p>
 obs(1)  = DCONJG(psi(1 )) * (psi(2  ) - psi(1   ))/dx
 do i    = 2, Nx - 1
  obs(i) = DCONJG(psi(i )) * (psi(i+1) - psi(i -1))/(2.0D0*dx)
 end do
 obs(Nx) = DCONJG(psi(Nx)) * (psi(Nx ) - psi(Nx-1))/dx
 call integrate(obs,dx,Nx,integral)
 avP     = DIMAG(integral) ! multiply by -i
!---------------------------------------------------------
!<x^2>
 obs = psi * x * x * DCONJG(psi)
 call integrate(obs,dx,Nx,integral)
 avX2= DBLE(integral)
!---------------------------------------------------------
!<p^2>
 obs(1 ) = DCONJG(psi(1 ))*(2.0D0*psi(1 )-5.0D0*psi(2   )+4.0D0*psi(3   )-psi(4   ))/(dx*dx)
 do i=2,Nx-1
  obs(i) = DCONJG(psi(i ))*(psi(i+1)-2.0D0*psi(i)+psi(i-1))/(dx*dx)
 enddo
 obs(Nx) = DCONJG(psi(Nx))*(2.0D0*psi(Nx)-5.0D0*psi(Nx-1)+4.0D0*psi(Nx-2)-psi(Nx-3))/(dx*dx)
 call integrate(obs,dx,Nx,integral)
 avP2=-DBLE(integral)
!---------------------------------------------------------
!<V(x)>
 obs = psi * V * DCONJG(psi)
 call integrate(obs,dx,Nx,integral)
 avV = DBLE(integral)
!---------------------------------------------------------
!<E>
 avE = 0.5D0 * avP2 + avV
!---------------------------------------------------------
!Uncertainties: DX2, DP2:
 DX2 = avX2  - avX  * avX
 DP2 = avP2  - avP  * avP
!---------------------------------------------------------
 DEALLOCATE(obs,x)
end subroutine observables
!===========================================================
!
!Simpson's rule to integrate psi(x).
!For n intervals of width dx (n even)
!Simpson's rule is:
!int(f(x)dx) = 
! (dx/3)*(f(x_0)+4 f(x_1)+2 f(x_2)+...+4 f(x_{n-1})+f(x_n))
!
!Input:   Discrete values of function psi(Nx)
!         Integration step dx
!Returns: Integral(psi(x) dx)
!===========================================================
!---------------------------------------------------------
subroutine     integrate(psi,dx,N,integral)
 implicit none
 integer    :: N
 complex(8) :: psi(N), integral, int
 real   (8) :: dx
 integer    :: Nx, i
!---------------------------------------------------------
!Nx must be odd:
 Nx = N; if(MOD(N,2) == 0) Nx = N-1
!---------------------------------------------------------
!zeroth order point:
 i    = 1
 int  = psi(i)
!---------------------------------------------------------
!odd  order points (i=k+1 is even):
 do i = 2,Nx-1,2
  int = int + 4.0D0*psi(i)
 enddo
!---------------------------------------------------------
!even order points:
 do i=3,Nx-2,2
  int = int + 2.0D0*psi(i)
 enddo
!---------------------------------------------------------
!last point:
 i    = Nx
 int  = int + psi(i)
!---------------------------------------------------------
!measure normalization:
 int  = int*dx/3.0D0
!---------------------------------------------------------
!If N is even, we add the last point:
 if(MOD(N,2) == 0) int  = int + psi(N)*dx
!---------------------------------------------------------
!----- final result:
 integral = int

end subroutine integrate
!---------------------------------------------------------
