!!----------------------------------------------------------------------------
!! $Id: compute_energy.f90,v 1.12 1997/11/12 02:03:48 svmcguir Exp $
!! subroutine compute_energy
!!    This subroutine will compute the energy of the entire system 
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------
real*8 function compute_energy(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

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

interface
   real*8 function compute_potential_energy(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
   end function

   real*8 function compute_kinetic_energy(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
   end function

end interface

   potential_energy = compute_potential_energy(max_params,params,payoff,eta)
   kinetic_energy   = compute_kinetic_energy(max_params,params,payoff,eta)

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

!!----------------------------------------------------------------------------
!! function compute_potential_energy 
!!    This subroutine will initialize the spins
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------
real*8 function compute_potential_energy(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                                        :: j,k, stock
real*8                                         :: sum , prod
   compute_potential_energy = 0.0d0
   do stock  = 1 , int(params(P_NPATH))
      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_potential_energy = compute_potential_energy + sum
   end do !stock

return
end function compute_potential_energy
!!----------------------------------------------------------------------------
!! function compute_kinetic_energy 
!!    This subroutine will initialize the spins
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!----------------------------------------------------------------------------

real*8 function compute_kinetic_energy(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                                        :: i, j, k, KI, KJ, JI, JJ
real*8                                         :: sum

   compute_kinetic_energy = 0.0d0

   do k = 1, int(params(P_NPATH))
      do j = 1, int(params(P_NTIME))
         do i = 1, eta(j,k)%nnn 
            JI = eta(j,k)%nni(i)
            JJ = eta(j,k)%nnj(i)
            compute_kinetic_energy=compute_kinetic_energy +                &
                 (eta(j,k)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0) &
                 / eta(j,k)%nnn

         end do !i
      end do !j
   end do !k

return
end function compute_kinetic_energy
