FWHSurfaceIntegral Function

   Function FWHSurfaceIntegral(self, f, isSolid) result(Pacc)

       use FWHDefinitions, only: rho0, P0, c0, U0, M0
       use VariableConversion, only: Pressure, PressureDot
       use fluiddata, only: dimensionless
       implicit none

       class(ObserverSourcePairClass)                      :: self
       class(Face), intent(in)                             :: f
       logical, intent(in)                                 :: isSolid
       real(kind=RP),dimension(3)                          :: Pacc  ! acoustic pressure values
       ! logical, intent(in)                                 :: isSolid, interpolate
       ! integer, intent(in), optional                       :: bufferPosition

       ! local variables
       integer                                             :: i, j  ! face indexes
       real(kind=RP), dimension(NDIM)                      :: Qi,QiDot, n
       real(kind=RP), dimension(NDIM,NDIM)                 :: Lij, LijDot
       type(NodalStorage_t), pointer                       :: spAxi, spAeta
       real(kind=RP)                                       :: Pt, Pl
       real(kind=RP)                                       :: LR, MR, UmMr,LdotR, LM
       ! integer                                             :: storePosition

       ! Initialization
       Pt = 0.0_RP
       Pl = 0.0_RP
       spAxi  => NodalStorage(f % Nf(1))
       spAeta => NodalStorage(f % Nf(2))


       associate( Q => f % storage(self % elementSide) % Q )
           associate( Qdot => f % storage(self % elementSide) % Qdot )

    !           **********************************
    !           Computes the surface integral
    !              I = \int vec{f}·vec{n} * vec{g}·vec{r} dS
    !           **********************************
    !
                do j = 0, f % Nf(2) ;    do i = 0, f % Nf(1)

                   n = f % geom % normal(:,i,j) * self % normalCorrection
                   call calculateFWHVariables(Q(:,i,j), Qdot(:,i,j), isSolid, Qi, QiDot, Lij, LijDot)

                   LR = dot_product(matmul(Lij, n(:)), self%reUnitVect(:,i,j))
                   MR = dot_product(M0(:), self % reUnitVect(:,i,j))
                   UmMr = 1 - MR
                   LdotR = dot_product(matmul(LijDot, n(:)), self%reUnitVect(:,i,j))
                   LM = dot_product(matmul(Lij, n(:)), M0(:))

                   ! loading term integrals
                   Pl = Pl +  dot_product(matmul(LijDot,n(:)),self%reUnitVect(:,i,j)) / (self%reStar(i,j) * c0) * &
                             spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   Pl = Pl +  dot_product(matmul(Lij,n(:)),self%reStarUnitVect(:,i,j)) / (self%reStar(i,j)**2) * &
                             spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + LdotR / ( c0 * self % re(i,j) * (UmMr**2) ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + (LR - LM) / ( (self % re(i,j) * UmMr)**2 ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                   ! Pl = Pl + (LR * (MR - (dimensionless % Mach**2))) / ( (self % re(i,j)**2) * (UmMr**3) ) * &
                   !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)

                   ! thickness term integrals, only for permeable surfaces
                   if (.not. isSolid) then
                       Pt = Pt + (1 - dot_product(M0(:),self%reUnitVect(:,i,j))) * dot_product(QiDot(:),n(:)) / (self%reStar(i,j)) * &
                                 spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       Pt = Pt -  dot_product(U0(:),self%reStarUnitVect(:,i,j)) * dot_product(Qi(:),n(:)) / (self%reStar(i,j)**2) * &
                                 spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       ! Pt = Pt + dot_product(QiDot(:),n(:)) / (self%reStar(i,j) * (UmMr**2)) * &
                       !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)
                       ! Pt = Pt + (dot_product(Qi(:), n(:)) * c0 * (MR - (dimensionless % Mach**2))) / ( (self % re(i,j)**2) * (UmMr**3) ) * &
                       !           spAxi % w(i) * spAeta % w(j) * f % geom % jacobian(i,j)

                   end if  
                end do          ;    end do
           end associate
       end associate

       Pt = Pt / (4.0_RP * PI)
       Pl = Pl / (4.0_RP * PI)

      ! get total acoustic pressure as the sum of the two components (the quadrapol terms are being ignored)
       Pacc = (/Pt, Pl, Pt+Pl/)

   End Function FWHSurfaceIntegral