Subroutine calculateFWHVariables(Q, Qdot, isSolid, Qi, QiDot, Lij, LijDot)
use VariableConversion, only: Pressure, PressureDot
use FWHDefinitions, only: rho0, P0, c0, U0, M0
use Utilities, only: AlmostEqual
implicit none
real(kind=RP), dimension(NCONS), intent(in) :: Q ! horses variables array
real(kind=RP), dimension(NCONS), intent(in) :: Qdot ! horses time derivatives array
logical, intent(in) :: isSolid
real(kind=RP), dimension(NDIM), intent(out) :: Qi ! fwh Qi array, related with the acoustic pressure thickness
real(kind=RP), dimension(NDIM), intent(out) :: Qidot
real(kind=RP), dimension(NDIM,NDIM), intent(out) :: Lij ! fwh Lij tensor: related with the acoustic pressure loading
real(kind=RP), dimension(NDIM,NDIM), intent(out) :: LijDot
!local variables
real(kind=RP) :: P, pDot
real(kind=RP), dimension(NDIM,NDIM) :: Pij ! fwh perturbation stress tensor
! real(kind=RP), dimension(NDIM:NDIM) :: tau
integer :: i, j, ii, jj
P = Pressure(Q)
pDot = PressureDot(Q,Qdot)
Pij = 0.0_RP
LijDot = 0.0_RP
do i=1,NDIM
Pij(i,i) = P - P0
!pressure derivative of LijDot
LijDot(i,i) = pDot
end do
!TODO use the stress tensor and the time derivative for Lij and LijDot respectively
! call getStressTensor(Q, U_x, U_y, U_z, tau)
! Pij = Pij - tau
! LijDot = LijDot - tauDot
! set values for solid (impermeable) surface
Qi(:) = -rho0*U0(:)
Qidot = 0.0_RP
Lij = Pij
! Lij = 0.0_RP
!calculate terms for permeable surface
if (.not. isSolid) then
Qi(:) = Qi(:) + Q(2:4)
! convert to complete velocity instead of perturbation velocity
QiDot(:) = QiDot(:) + Qdot(2:4)
do j = 1, NDIM
jj = j + 1
do i = 1, NDIM
! one index is added since rhoV1 = Q(2), rhoV2 = Q(3) ...
ii = i + 1
Lij(i,j) = Lij(i,j) + (Q(ii) - Q(1)*U0(i))*(Q(jj)/Q(1))
LijDot(i,j) = LijDot(i,j) + ( Qdot(ii) - Q(ii)/Q(1)*Qdot(1) )/Q(1) * Q(jj) + &
(Q(ii)/Q(1) - U0(i)) * Qdot(jj)
end do
end do
end if
End Subroutine calculateFWHVariables
