! File:         submodule_sph_particles_constructor_bin.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 (sph_particles) constructor_bin

  !************************************************
  !
  !# This SUBMODULE contains the implementation
  !  of the constructor of TYPE [[particles]]
  !  from an |id| binary file.
  !
  !  FT 01.03.2022
  !
  !************************************************


  IMPLICIT NONE


  CONTAINS


  MODULE PROCEDURE construct_particles_bin

    !**************************************************
    !
    !# Constructs a [[particles]] object from an |id|
    !  binary file
    !
    !  FT 01.03.2022
    !
    !**************************************************


    USE tensor,         ONLY: n_sym4x4
    USE utility,        ONLY: compute_g4, determinant_sym4x4, &
                              spacetime_vector_norm_sym4x4,zero, one
    USE input_output,   ONLY: read_options
    USE options,        ONLY: eos_str
    USE units,          ONLY: set_units
    USE pwp_EOS,        ONLY: select_EOS_parameters
    USE sph_variables,  ONLY: allocate_SPH_memory, deallocate_SPH_memory


    IMPLICIT NONE


    INTEGER:: a

    DOUBLE PRECISION:: g4(n_sym4x4), det, v_norm
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp1
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp2
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp3
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp4
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp5
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp6
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp7
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp8
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp9
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: tmp10

    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: sq_detg4


    PRINT *, "** Assigning the ID to the particles..."
    PRINT *

    CALL set_units('NSM')
    CALL read_options

    CALL parts% read_sphincs_dump_print_formatted( namefile_bin= namefile, &
                                                   save_data   = .TRUE. )

    ALLOCATE( tmp1 ( parts% npart ) )
    ALLOCATE( tmp2 ( parts% npart ) )
    ALLOCATE( tmp3 ( parts% npart ) )
    ALLOCATE( tmp4 ( parts% npart ) )
    ALLOCATE( tmp5 ( parts% npart ) )
    ALLOCATE( tmp6 ( parts% npart ) )
    ALLOCATE( tmp7 ( parts% npart ) )
    ALLOCATE( tmp8 ( parts% npart ) )
    ALLOCATE( tmp9 ( parts% npart ) )
    ALLOCATE( tmp10( parts% npart ) )

    CALL id% read_id_particles( parts% npart, &
                                parts% pos( 1, : ), &
                                parts% pos( 2, : ), &
                                parts% pos( 3, : ), &
                                parts% lapse, &
                                parts% shift_x, &
                                parts% shift_y, &
                                parts% shift_z, &
                                parts% g_xx, &
                                parts% g_xy, &
                                parts% g_xz, &
                                parts% g_yy, &
                                parts% g_yz, &
                                parts% g_zz, &
                                tmp4, &
                                tmp5, &
                                tmp6, &
                                tmp7, &
                                tmp8, &
                                tmp9, &
                                tmp10 )

    DEALLOCATE( tmp1  )
    DEALLOCATE( tmp2  )
    DEALLOCATE( tmp3  )
    DEALLOCATE( tmp4  )
    DEALLOCATE( tmp5  )
    DEALLOCATE( tmp6  )
    DEALLOCATE( tmp7  )
    DEALLOCATE( tmp8  )
    DEALLOCATE( tmp9  )
    DEALLOCATE( tmp10 )

    ALLOCATE( sq_detg4( parts% npart ) )

    DO a= 1, parts% npart, 1

      CALL compute_g4( parts% lapse(a), &
            [parts% shift_x(a), parts% shift_y(a), parts% shift_z(a)], &
            [parts% g_xx(a), parts% g_xy(a), parts% g_xz(a), &
             parts% g_yy(a), parts% g_yz(a), parts% g_zz(a)], &
             g4(1:n_sym4x4) )

      CALL determinant_sym4x4( g4(1:n_sym4x4), det )
      IF( ABS(det) < 1D-10 )THEN
          PRINT *, "** ERROR! The determinant of the spacetime metric is " &
                   // "effectively 0 at particle ", a
          PRINT *
          STOP
      ELSEIF( det > 0 )THEN
          PRINT *, "** ERROR! The determinant of the spacetime metric is " &
                   // "positive at particle ", a
          PRINT *
          STOP
      ENDIF
      sq_detg4(a)= SQRT(-det)

      IF( parts% theta(a) < zero )THEN
        PRINT *, "** ERROR! The computing frame generalized Lorentz factor ", &
                 "is negative at particle  ", a
        PRINT *, " * Its value is ", parts% theta(a)
        PRINT *, " * Stopping.."
        PRINT *
        STOP
      ELSEIF( parts% theta(a) == zero )THEN
        PRINT *, "** ERROR! The computing frame generalized Lorentz factor ", &
                 "is zero at particle ", a
        PRINT *, " * Its value is ", parts% theta(a)
        PRINT *, " * Stopping.."
        PRINT *
        STOP
      ENDIF

      v_norm= zero
      CALL spacetime_vector_norm_sym4x4( g4, parts% v(:,a), v_norm )
      IF( v_norm > zero )THEN
        PRINT *, "** ERROR! The computing frame particle 4-velocity is ", &
                 "spacelike at particle ", a
        PRINT *, " * Its norm is ", v_norm
        PRINT *, " * Stopping.."
        PRINT *
        STOP
      ELSEIF( v_norm == zero )THEN
        PRINT *, "** ERROR! The computing frame particle 4-velocity is ", &
                 "null at particle ", a
        PRINT *, " * Its norm is ", v_norm
        PRINT *, " * Stopping.."
        PRINT *
        STOP
      ENDIF

    ENDDO

    parts% nstar_sph= parts% nlrf_sph*sq_detg4*parts% theta

    CALL allocate_sph_memory

    !CALL OMP_SET_NUM_THREADS(1)

    CALL select_EOS_parameters( eos_str )

    CALL parts% test_recovery( parts% npart,       &
                               parts% pos,         &
                               parts% nlrf_sph,    &
                               parts% u_sph,       &
                               parts% pressure_sph, &
                               parts% v(1:3,:),    &
                               parts% theta,       &
                               parts% nstar_sph )

    DEALLOCATE( sq_detg4 )
    CALL deallocate_sph_memory

  END PROCEDURE construct_particles_bin


END SUBMODULE constructor_bin