! File:         submodule_bns_fuka_interpolate.f90
! Authors:      Francesco Torsello (FT)
!************************************************************************
! Copyright (C) 2020-2023 Francesco Torsello                            *
!                                                                       *
! This file is part of SPHINCS_ID                                       *
!                                                                       *
! SPHINCS_ID is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by  *
! the Free Software Foundation, either version 3 of the License, or     *
! (at your option) any later version.                                   *
!                                                                       *
! SPHINCS_ID is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of        *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
! GNU General Public License for more details.                          *
!                                                                       *
! You should have received a copy of the GNU General Public License     *
! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>.   *
! The copy of the GNU General Public License should be in the file      *
! 'COPYING'.                                                            *
!************************************************************************

SUBMODULE (bns_fuka) interpolate

  !****************************************************
  !
  !# Implementation of the methods of TYPE [[bnsfuka]] that
  !  interpolate the data on a lattice, at a particle
  !  position
  !
  !  FT 28.06.2022
  !
  !****************************************************


  USE utility,  ONLY: zero, one


  IMPLICIT NONE


  CONTAINS


  !-------------------!
  !--  SUBROUTINES  --!
  !-------------------!


  MODULE PROCEDURE interpolate_fuka_id_particles

    !****************************************************
    !
    !# Stores the hydro |id| in the arrays needed to
    !  compute the |sph| |id|
    !
    !  FT 28.06.2022
    !
    !****************************************************

    USE constants,  ONLY: MSun, amu
    USE utility,    ONLY: two
    USE numerics,   ONLY: trilinear_interpolation

    IMPLICIT NONE

    LOGICAL, PARAMETER:: debug= .FALSE.
    INTEGER:: a, star
    DOUBLE PRECISION:: zp
    !DOUBLE PRECISION, &
    !  DIMENSION(this% nx_lattice, this% ny_lattice, this% nz_lattice, 3):: coords

    !$OMP PARALLEL DO DEFAULT( NONE ) &
    !$OMP             SHARED( n, this, x, y, z, lapse, &
    !$OMP                     shift_x, shift_y, shift_z, &
    !$OMP                     g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, &
    !$OMP                     baryon_density, energy_density, &
    !$OMP                     specific_energy, pressure, &
    !$OMP                     u_euler_x, u_euler_y, u_euler_z ) &
    !$OMP             PRIVATE( a, star, zp )
    DO a= 1, n, 1

      IF( (this% center(1,1) - this% radii(1,1) <= x(a)) &
          .AND. &
          (x(a) <= this% center(1,1) + this% radii(1,2)) )THEN

        star= 1

      ELSEIF( (this% center(2,1) - this% radii(2,1) <= x(a)) &
              .AND. &
              (x(a) <= this% center(2,1) + this% radii(2,2)) )THEN

        star= 2

      ELSE

        star= -1

      ENDIF

      IF( star == -1 )THEN
        baryon_density(a) = zero
        specific_energy(a)= zero
        pressure(a)       = zero
        u_euler_x(a)      = zero
        u_euler_y(a)      = zero
        u_euler_z(a)      = zero
        energy_density(a) = zero
        g_xx(a)           = zero
        g_yy(a)           = zero
        g_zz(a)           = zero
        g_xy(a)           = zero
        g_xz(a)           = zero
        g_yz(a)           = zero
        lapse(a)          = zero
        shift_x(a)        = zero
        shift_y(a)        = zero
        shift_z(a)        = zero
        CYCLE
      ENDIF

      !coords= this% star_lattice(star)% coords

      zp= z(a)

      baryon_density(a) = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% mass_density, &
                                debug= .FALSE. ) &
                                *MSun/amu

      specific_energy(a)= trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% specific_energy )

      pressure(a)       = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% pressure ) &
                                *MSun/amu

      u_euler_x(a)      = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% v_eul_x )
      u_euler_y(a)      = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% v_eul_y )
      u_euler_z(a)      = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% v_eul_z )

      IF( baryon_density(a) == zero )THEN
        specific_energy(a)= zero
        u_euler_x(a)      = zero
        u_euler_y(a)      = zero
        u_euler_z(a)      = zero
      ENDIF

      energy_density(a) = baryon_density(a)*(one + specific_energy(a))

      g_xx(a)           = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% g_xx )

      g_yy(a)= g_xx(a)
      g_zz(a)= g_xx(a)
      g_xy(a)= zero
      g_xz(a)= zero
      g_yz(a)= zero

      lapse(a)          = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% lapse )
      shift_x(a)        = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% shift_x )
      shift_y(a)        = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% shift_y )
      shift_z(a)        = trilinear_interpolation( x(a), y(a), zp, &
                                this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                this% star_lattice(star)% coords, &
                                this% star_lattice(star)% shift_z )

    ENDDO
    !$OMP END PARALLEL DO

  END PROCEDURE interpolate_fuka_id_particles


  MODULE PROCEDURE interpolate_fuka_id_mass_b

    !****************************************************
    !
    !# Stores the hydro |id| in the arrays needed to
    !  compute the baryon mass, storing it to variables
    !  (not arrays as the others SUBROUTINES in
    !  the [[bnsfuka@interpolate]] SUBMODULE).
    !
    !  FT 28.06.2022
    !
    !****************************************************

    USE tensor,    ONLY: jxx, jxy, jxz, jyy, jyz, jzz
    USE constants, ONLY: MSun, amu
    USE numerics,   ONLY: trilinear_interpolation

    IMPLICIT NONE

    INTEGER:: star
    DOUBLE PRECISION:: zp, veuler_x, veuler_y, veuler_z
    !DOUBLE PRECISION, &
    !  DIMENSION(this% nx_lattice, this% ny_lattice, this% nz_lattice, 3):: coords

    IF( (this% center(1,1) - this% radii(1,1) <= x) &
        .AND. &
        (x <= this% center(1,1) + this% radii(1,2)) )THEN

      star= 1

    ELSEIF( (this% center(2,1) - this% radii(2,1) <= x) &
            .AND. &
            (x <= this% center(2,1) + this% radii(2,2)) )THEN

      star= 2

    ELSE

      star= -1

    ENDIF

    IF( star == -1 )THEN
      baryon_density= zero
      g(jxx)        = zero
      g(jyy)        = zero
      g(jzz)        = zero
      g(jxy)        = zero
      g(jxz)        = zero
      g(jyz)        = zero
      gamma_euler   = zero
      RETURN
    ENDIF

    !coords= this% star_lattice(star)% coords

    zp= z

    baryon_density= trilinear_interpolation( x, y, zp, &
                          this% nx_lattice, this% ny_lattice, this% nz_lattice,&
                          this% star_lattice(star)% coords, &
                          this% star_lattice(star)% mass_density )!*MSun/amu

    g(jxx)= trilinear_interpolation( x, y, zp, &
                  this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                  this% star_lattice(star)% coords, &
                  this% star_lattice(star)% g_xx )
    g(jyy)= g(jxx)
    g(jzz)= g(jxx)
    g(jxy)= zero
    g(jxz)= zero
    g(jyz)= zero

    veuler_x= trilinear_interpolation( x, y, zp, &
                     this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                     this% star_lattice(star)% coords, &
                     this% star_lattice(star)% v_eul_x )
    veuler_y= trilinear_interpolation( x, y, zp, &
                     this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                     this% star_lattice(star)% coords, &
                     this% star_lattice(star)% v_eul_y )
    veuler_z= trilinear_interpolation( x, y, zp, &
                     this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                     this% star_lattice(star)% coords, &
                     this% star_lattice(star)% v_eul_z )

    ! See eq.(7.3.13) in Alcubierre, "Introduction to 3+1 Numerical Relativity"
    ! The following formula assumes a conformally flat metric in Cartesian
    ! coordinates
    gamma_euler= one/SQRT( one - g(jxx)*( veuler_x*veuler_x &
                                        + veuler_y*veuler_y &
                                        + veuler_z*veuler_z ) );

  END PROCEDURE interpolate_fuka_id_mass_b


  !-----------------!
  !--  FUNCTIONS  --!
  !-----------------!


  MODULE PROCEDURE interpolate_fuka_mass_density

    !***********************************************
    !
    !# Returns the mass density at the point
    !  given as argument, in units of
    !  \(M_\odot/L_\odot^3\).
    !
    !  FT 28.06.2022
    !
    !***********************************************

    USE timing,    ONLY: timer
    USE numerics,  ONLY: trilinear_interpolation
    USE utility,   ONLY: spherical_from_cartesian, two


    IMPLICIT NONE

    INTEGER:: star
    DOUBLE PRECISION:: zp

    !TYPE(timer):: dbg_timer

    !dbg_timer= timer("dbg_timer")
    !CALL dbg_timer% start_timer()

    IF( (this% center(1,1) - this% radii(1,1) <= x) &
        .AND. &
        (x <= this% center(1,1) + this% radii(1,2)) )THEN

      star= 1

    ELSEIF( (this% center(2,1) - this% radii(2,1) <= x) &
            .AND. &
            (x <= this% center(2,1) + this% radii(2,2)) )THEN

      star= 2

    ELSE

      star= -1

    ENDIF

    IF( star == -1 )THEN
      res= zero
      RETURN
    ENDIF

    zp= z

    res= trilinear_interpolation( x, y, zp, &
                                  this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                  this% star_lattice(star)% coords, &
                                  this% star_lattice(star)% mass_density )
    !CALL dbg_timer% stop_timer()
    !CALL dbg_timer% print_timer( 2 )
    !STOP
  END PROCEDURE interpolate_fuka_mass_density


  MODULE PROCEDURE interpolate_fuka_spatial_metric

    !***********************************************
    !
    !# Returns the spatial metric.
    !
    !  FT 28.06.2022
    !
    !***********************************************

    USE constants, ONLY: pi
    USE numerics,  ONLY: trilinear_interpolation
    USE utility,   ONLY: spherical_from_cartesian, two

    IMPLICIT NONE

    INTEGER:: star
    DOUBLE PRECISION:: zp

    IF( (this% center(1,1) - this% radii(1,1) <= x) &
        .AND. &
        (x <= this% center(1,1) + this% radii(1,2)) )THEN

      star= 1

    ELSEIF( (this% center(2,1) - this% radii(2,1) <= x) &
            .AND. &
            (x <= this% center(2,1) + this% radii(2,2)) )THEN

      star= 2

    ELSE

      star= -1

    ENDIF

    IF( star == -1 )THEN
      res= zero
      RETURN
    ENDIF

    zp= z
    res= trilinear_interpolation( x, y, zp, &
                                  this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                  this% star_lattice(star)% coords, &
                                  this% star_lattice(star)% g_xx )

  END PROCEDURE interpolate_fuka_spatial_metric


  MODULE PROCEDURE interpolate_fuka_pressure

    !***********************************************
    !
    !# Returns the spatial metric.
    !
    !  FT 28.06.2022
    !
    !***********************************************

    USE constants, ONLY: pi
    USE numerics,  ONLY: trilinear_interpolation
    USE utility,   ONLY: spherical_from_cartesian, two

    IMPLICIT NONE

    INTEGER:: star
    DOUBLE PRECISION:: zp

    IF( (this% center(1,1) - this% radii(1,1) <= x) &
        .AND. &
        (x <= this% center(1,1) + this% radii(1,2)) )THEN

      star= 1

    ELSEIF( (this% center(2,1) - this% radii(2,1) <= x) &
            .AND. &
            (x <= this% center(2,1) + this% radii(2,2)) )THEN

      star= 2

    ELSE

      star= -1

    ENDIF

    IF( star == -1 )THEN
      res= zero
      RETURN
    ENDIF

    zp= z
    res= trilinear_interpolation( x, y, zp, &
                                  this% nx_lattice, this% ny_lattice, this% nz_lattice, &
                                  this% star_lattice(star)% coords, &
                                  this% star_lattice(star)% pressure )

  END PROCEDURE interpolate_fuka_pressure


  MODULE PROCEDURE is_hydro_positive_interpolation

    !************************************************
    !
    !# Return 1 if the energy density is nonpositive
    !  or if the specific energy is nonpositive,
    !  or if the pressure is nonpositive
    !  at the specified point; return 0 otherwise
    !
    !  FT 28.06.2022
    !
    !************************************************

    IMPLICIT NONE

    LOGICAL:: outside_y, outside_z!, outside_x

    !outside_x=.NOT.( &
    !            (
    !              this% center(1,1) - this% radii(1,1) <= x &
    !              .AND. &
    !              x <= this% center(1,1) + this% radii(1,2) &
    !            ) &
    !            .OR. &
    !            ( &
    !              this% center(2,1) - this% radii(2,1) <= x &
    !              .AND. &
    !              x <= this% center(2,1) + this% radii(2,2) &
    !            ) &
    !          )

    outside_y=.NOT.( &
                ( &
                  this% center(1,2) - this% radii(1,3) <= y &
                  .AND. &
                  y <= this% center(1,2) + this% radii(1,4) &
                ) &
                .OR. &
                ( &
                  this% center(2,2) - this% radii(2,3) <= y &
                  .AND. &
                  y <= this% center(2,2) + this% radii(2,4) &
                ) &
              )

    outside_z=.NOT.( &
                ( &
                  this% center(1,3) - this% radii(1,5) <= z &
                  .AND. &
                  z <= this% center(1,3) + this% radii(1,6) &
                ) &
                .OR. &
                ( &
                  this% center(2,3) - this% radii(2,5) <= z &
                  .AND. &
                  z <= this% center(2,3) + this% radii(2,6) &
                ) &
              )

    IF( this% read_mass_density( x, y, z ) <= zero &
        !.OR. outsize_x &
        .OR. outside_y &
        .OR. outside_z &
    )THEN
      res= .FALSE.
    ELSE
      res= .TRUE.
    ENDIF

  END PROCEDURE is_hydro_positive_interpolation


END SUBMODULE interpolate