!!----------------------------------------------------------------------------
!! $Id: compute_local_energy.f90,v 1.12 1997/11/12 02:03:48 svmcguir Exp $ 
!! subroutine compute_local_energy
!!    This subroutine will compute the local energy 
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------
real*8 function compute_local_energy(t,p,max_params,params,payoff,eta)
use PATHS
use input_params
implicit none
integer,intent(IN)                             :: max_params
real*8,dimension(max_params),intent(IN)        :: params
type (PATH),dimension(:,:),intent(OUT)         :: eta
real*8,dimension(:,:),intent(IN)               :: payoff
integer                                        :: t,p

real*8                                         :: potential_energy 
real*8                                         :: kinetic_energy 

interface
   real*8 function compute_lpotential_energy(t,p,max_params,params,payoff,eta)
   use PATHS
   use input_params
   implicit none
   integer,intent(IN)                             :: max_params
   real*8,dimension(max_params),intent(IN)        :: params
   type (PATH),dimension(:,:),intent(OUT)         :: eta
   real*8,dimension(:,:),intent(IN)               :: payoff
   integer                                        :: t,p
   end function

   real*8 function compute_lkinetic_energy(t,p,max_params,params,payoff,eta)
   use PATHS
   use input_params
   implicit none
   integer,intent(IN)                             :: max_params
   real*8,dimension(max_params),intent(IN)        :: params
   type (PATH),dimension(:,:),intent(OUT)         :: eta
   real*8,dimension(:,:),intent(IN)               :: payoff
   integer                                        :: t,p
   end function
   
end interface

   potential_energy=compute_lpotential_energy(t,p,max_params,params,payoff,eta)
   kinetic_energy  =compute_lkinetic_energy(t,p,max_params,params,payoff,eta)

   compute_local_energy = -params(P_C1)*potential_energy                      &
                          -params(P_C2)*kinetic_energy
return
end function compute_local_energy

!!----------------------------------------------------------------------------
!! function compute_plotential_energy 
!!    This subroutine will compute the local potential energy 
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------
real*8 function compute_lpotential_energy(t,p,max_params,params,payoff,eta)
use PATHS
use input_params
implicit none
integer,intent(IN)                             :: max_params
real*8,dimension(max_params),intent(IN)        :: params
type (PATH),dimension(:,:),intent(OUT)         :: eta
real*8,dimension(:,:),intent(IN)               :: payoff
integer                                        :: t,p

integer                                        :: j,k, stock
real*8                                         :: sum , prod

compute_lpotential_energy = 0.0d0
stock = p
sum = 0.0d0
do k = 1 , int(params(P_NTIME))
   if ( (eta(k,stock)%spin<>(0.0d0)).and.(payoff(k,stock)<>(0.0d0))) then
      prod = 1.0d0
      do j = 1 , k - 1
         prod = prod* (1.0d0-eta(j,stock)%spin)
      end do !j
      sum = sum + prod*eta(k,stock)%spin * payoff(k,stock)
   end if
end do !k
compute_lpotential_energy = compute_lpotential_energy + sum

return
end function compute_lpotential_energy
!!----------------------------------------------------------------------------
!! function compute_kinetic_energy 
!!    This subroutine will compute the change in the local potential energy 
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------

real*8 function compute_lkinetic_energy(t,p,max_params,params,payoff,eta)
use PATHS
use input_params
implicit none
integer,intent(IN)                             :: max_params
real*8,dimension(max_params),intent(IN)        :: params
type (PATH),dimension(:,:),intent(OUT)         :: eta
real*8,dimension(:,:),intent(IN)               :: payoff
integer                                        :: t,p

integer                                        :: i, j, k, KI, KJ, JI, JJ
real*8                                         :: sum
integer                                        :: nnl

   compute_lkinetic_energy = 0.0d0

!! For local energy changes only loop over the NN list of (t,p) termi

   do nnl = 1 , eta(t,p)%nnn 
      j = eta(t,p)%nni(nnl)
      k = eta(t,p)%nnj(nnl)
      do i = 1, eta(j,k)%nnn 
         JI = eta(j,k)%nni(i)
         JJ = eta(j,k)%nnj(i)
         compute_lkinetic_energy=compute_lkinetic_energy +                  &
              (eta(j,k)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
      end do !i
   end do !nnl
! now we must change the energy of the t,p point too... since this is not in the
! nn list
   j = t
   k = p
   do i = 1, eta(j,k)%nnn 
      JI = eta(j,k)%nni(i)
      JJ = eta(j,k)%nnj(i)
      compute_lkinetic_energy=compute_lkinetic_energy +                  &
           (eta(j,k)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
   end do !i

   compute_lkinetic_energy=compute_lkinetic_energy / eta(j,k)%nnn

return
end function compute_lkinetic_energy
