particles Derived Type

  TYPE:: particles
  !! TYPE representing a |sph| particle distribution


    PRIVATE


    INTEGER:: npart
    !! Total particle number
    INTEGER:: n_matter
    !! Number of matter objects in the physical system
    INTEGER, DIMENSION(:), ALLOCATABLE:: npart_i
    !! Array storing the particle numbers for the matter objects
    INTEGER, DIMENSION(:), ALLOCATABLE:: npart_fin
    !! Array storing the index of the last particle for each matter objects
    INTEGER:: distribution_id
    !! Identification number for the particle distribution
    INTEGER:: call_flag= 0
    !# Flag that is set different than 0 if the SUBROUTINE
    !  compute_and_print_sph_variables is called
    LOGICAL:: cold_system
    !# `.TRUE.` if the system is at zero temperature (no thermal component);
    !  `.FALSE.` otherwise
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: barycenter
    !# Array storing the centers of mass of the matter objects
    DOUBLE PRECISION, DIMENSION(4):: barycenter_system
    !# Array storing the center of mass of the **entire particle distribution**

    INTEGER, DIMENSION(:), ALLOCATABLE:: baryon_density_index
    !# Array storing the indices to use with [[particles:baryon_density]]
    !  to sort the elements of [[particles:baryon_density]] in increasing
    !  order

    !
    !-- Hydro variables on the particles
    !

    !> 2-D array storing the particle positions
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: pos
    !& 1-D array storing the baryon mass density in the fluid frame
    !  \([\mathrm{kg}\,\mathrm{m}^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: baryon_density
    !> 1-D array storing the energy density
    !  \([\mathrm{kg}\,c^2\,\mathrm{m}^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: energy_density
    !> 1-D array storing the specific internal energy \([c^2]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: specific_energy
    !& 1-D array storing the specific internal energy \([c^2]\) computed using
    !  formula (9) in Read et al., Phys.Rev.D79:124032,2009,
    !  [arXiv:0812.2163][https://arxiv.org/abs/0812.2163]{:target="_blank"}
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: u_sph
    !> 1-D array storing the pressure \([\mathrm{kg}\,c^2\,\mathrm{m}^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: pressure
    !& 1-D array storing the pressure in code units
    !  \([\mathrm{amu}\,c^2\,\mathrm{L_\odot}^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: pressure_sph
    !& 1-D array storing the x component of the fluid 3-velocity wrt
    !  the Eulerian observer \([c]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: v_euler_x
    !& 1-D array storing the y component of the fluid 3-velocity wrt
    !  the Eulerian observer \([c]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: v_euler_y
    !& 1-D array storing the z component of the fluid 3-velocity wrt
    !  the Eulerian observer \([c]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: v_euler_z

    !
    !-- Arrays to store the electron fraction Ye as a function of the
    !-- baryon number density for beta-equilibrated EoS at T~0,
    !-- imported from the CompOSE database's and software's files
    !

    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nb_table
    !# Array storing the values of the baryon number density in the |compose|
    !  table. @todo ADD UNITS
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: Ye_table
    !! Array storing the values of the electron fraction in the |compose| table

    !
    !-- Spacetime fields
    !

    !> Array storing the values of the lapse function on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: lapse
    !& Array storing the values of the x component of the shift vector
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: shift_x
    !& Array storing the values of the y component of the shift vector
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: shift_y
    !& Array storing the values of the z component of the shift vector
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: shift_z
    !& Array storing the values of the xx component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_xx
    !& Array storing the values of the xy component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_xy
    !& Array storing the values of the xz component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_xz
    !& Array storing the values of the xz component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_yy
    !& Array storing the values of the yz component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_yz
    !& Array storing the values of the zz component of the spatial metric
    !  on the particles
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: g_zz

    !
    !-- SPH fields
    !

    !& 1-D array storing baryon density in the local rest frame
    !  \([\mathrm{baryon}\, (L_\odot)^{-3}]\), computed directly from
    !  the |lorene| density
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nlrf
    !& 1-D array storing baryon density in the local rest frame
    !  \([\mathrm{baryon}\, (L_\odot)^{-3}]\), computed from the kernel
    !  interpolated proper baryon number density [[particles:nstar_sph]]
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nlrf_sph
    !> 1-D array storing the baryon number per particle
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nu
    !& 1-D array storing the SPH estimate of the proper baryon number density
    !  \([\mathrm{baryon}\, (L_\odot)^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nstar
    !& 1-D array storing the particle number density
    !  \([\mathrm{particle}\, (L_\odot)^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: particle_density
    !& 1-D array storing the SPH estimate of the proper baryon number density,
    !  from kernel interpolation \([\mathrm{baryon}\, (L_\odot)^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nstar_sph
    !& 1-D array storing the SPH estimate of the particle number density, from
    !  kernel interpolation \([\mathrm{particle}\, (L_\odot)^{-3}]\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: particle_density_sph
    !> 1-D array storing the enthalpy \([??]\) @todo add units
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: enthalpy
    !> 2-D array storing the coordinate fluid 4-velocity \([c]\)
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: v
    !> 1-D array storing the generalized Lorentz factor
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: Theta
    !> 1-D array storing the electron fraction
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: Ye
    !> 1-D array storing the smoothing length \(L_\odot\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: h
    !> 1-D array storing the particle volumes \(L_\odot^3\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: pvol
    !> 1-D array storing the particle masses \(M_\odot\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: masses
    !> 1-D array storing the surfaces of the matter objects
    TYPE(surface),  DIMENSION(:), ALLOCATABLE:: surfaces
    !& Ratio between baryonic masses of the matter objects and the maximum
    !  baryonic mass among them @warning always \(< 1\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: mass_ratios
    !& Ratio between baryonic masses of the matter objects and the total
    !  baryonic mass of all the matter objects @warning always \(< 1\)
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: mass_fractions
    !> Total baryon number
    DOUBLE PRECISION:: nbar_tot
    !> Baryon numbers of the matter objects
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nbar_i
    !& Ratio between the max and min of the baryon number per particle, over
    !  all the matter objects
    DOUBLE PRECISION:: nuratio
    !& Desired ratio between the max and min of the baryon number per particle,
    !  over all the matter objects. **Only used when redistribute_nu is .TRUE.**
    !  @warning Almost deprecated
    DOUBLE PRECISION:: nu_ratio_des
    !> Baryon number ratios on the matter objects
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: nuratio_i

    DOUBLE PRECISION:: adm_mass
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: adm_linear_momentum_i
    !# Estimate of the \(\mathrm{ADM}\) linear momentum of each matter object,
    !  computed from the canonical momentum per baryon on the particles
    DOUBLE PRECISION, DIMENSION(3)               :: adm_linear_momentum_fluid
    !# Estimate of the \(\mathrm{ADM}\) linear momentum of the fluid computed
    !  from the canonical momentum per baryon on the particles

    !
    !-- Strings
    !

    CHARACTER(LEN= 50):: sphincs_id_particles
    !! String containing the name of the particles parameter file

    CHARACTER(LEN=:), ALLOCATABLE:: compose_path
    !# String storing the local path to the directory containing the
    !  CompOSE |eos|
    CHARACTER(LEN=:), ALLOCATABLE:: compose_filename
    !# String storing the subpath of compose_path to the CompOSE file with
    !  .beta extension

    !
    !-- Equations of state
    !

    TYPE(eos), DIMENSION(:), ALLOCATABLE:: all_eos
    !# Array of TYPE [[eos]] containing the |eos| information for all the
    !  matter objects

    !
    !-- Procedure pointers
    !

    !PROCEDURE(post_process_sph_id_int), POINTER, NOPASS:: post_process_sph_id
    PROCEDURE(), POINTER, NOPASS:: post_process_sph_id
    !# Pointer to a procedure that post_process the |sph| |id|; for example,
    !  correct for the residual ADM linear momentum.

    !
    !-- Steering variables
    !

    !> `.TRUE.` if the object is empty, `.FALSE.` if it's not empty
    LOGICAL:: empty_object
    !& `.TRUE.` if the binary files for SPHINCS_BSSN are to be exported,
    !  `.FALSE.` otherwise
    LOGICAL, PUBLIC:: export_bin
    !& `.TRUE.` if the ID in the formatted files is to be on the xy plane only,
    !  `.FALSE.` otherwise
    LOGICAL, PUBLIC:: export_form_xy
    !& `.TRUE.` if the ID in the formatted files is to be on the x axis only,
    !  `.FALSE.` otherwise
    LOGICAL, PUBLIC:: export_form_x
    !& `.TRUE.` if the threshold on the baryon mass density should e applied
    !  when placing particles on lattices, `.FALSE.` otherwise
    LOGICAL:: use_thres
    !& `.TRUE.` if the baryon number per particle should be reassigned, trying
    !  to obtain a baryon number ratio no larger than nu_ratio,
    !  when placing particles on lattices; `.FALSE.` otherwise
    LOGICAL:: redistribute_nu
    !& `.TRUE.` if the baryon number per particle should be corrected to account
    !  for the total baryon masses of the stars, `.FALSE.` otherwise
    LOGICAL:: correct_nu
    !& `.TRUE.` if the electron fraction \(Y_e\) should be read from the CompOSE
    !  table with extension.beta, `.FALSE.` otherwise
    !  @todo Chamge name of this variable to assign_Ye_compose. Check that
    !        the used EOS is indeed the one used to read \(Y_e\)
    LOGICAL:: compose_eos
    !& `.TRUE.` if the particle positions on ellipsoidal surfaces are randomized
    !  in the \(\phi\) direction, `.FALSE.` otherwise
    LOGICAL:: randomize_phi
    !& `.TRUE.` if the particle positions on ellipsoidal surfaces are randomized
    !  in the \(\theta\) direction, `.FALSE.` otherwise
    LOGICAL:: randomize_theta
    !& `.TRUE.` if the particle positions on ellipsoidal surfaces are randomized
    !  in the \(r\) direction, `.FALSE.` otherwise
    LOGICAL:: randomize_r
    !& `.TRUE.` if the Artificial Pressure Method (APM) has to be applied to the
    !  particles, `.FALSE.` otherwise
    LOGICAL, DIMENSION(:), ALLOCATABLE:: apm_iterate
    !& `.TRUE.` to allow the particles to move where the density is 0 during the
    !  Artificial Pressure Method (APM) iteration. This can be useful when the
    !  system has an irregular geometry, as, for example, an ejecta
    !  `.FALSE.` otherwise
    LOGICAL, DIMENSION(:), ALLOCATABLE:: use_atmosphere
    !& `.TRUE.` if the baryon number per particle \(\nu\) has to be read from
    !  the formatted file containing the particle positions, `.FALSE.` otherwise
    LOGICAL:: read_nu
    !& `.TRUE.` if the particles on star 2 should be the reflection of the
    !  particles on star 1 with respect to the \(yz\) plane, only if the baryon
    !  masses of the stars differe less than \(0.2\%\); `.FALSE.` otherwise
    LOGICAL:: reflect_particles_x

    !
    !-- Timers
    !

    !> Timer that times how long it takes to place particles on the stars
    TYPE(timer), PUBLIC:: placer_timer
    !& Timer that times how long it takes to check if there are multiple
    !  particles at the same positions
    TYPE(timer), PUBLIC:: same_particle_timer
    !& Timer that times how long it takes to perform the APM on the matter
    !  objects
    TYPE(timer), DIMENSION(:), ALLOCATABLE, PUBLIC:: apm_timers
    !& Timer that times how long it takes to import the \(\texttt{|lorene|}\) ID
    !  at the particle positions
    TYPE(timer), PUBLIC:: importer_timer
    !& Timer that times how long it takes to compute the SPH variables at the
    !  particle positions
    TYPE(timer), PUBLIC:: sph_computer_timer


    CONTAINS


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

    PROCEDURE:: read_particles_formatted_file
    !# Read particle positions and nu from a formatted file with the following
    !  format:
    !
    !   1. The first line must contain the integer number of objects for the
    !      system, and the particle numbers on each object;
    !      each column separated by one tab.
    !   2. The other lines must contain at least 3 columns with the x, y, z
    !      coordinates of the particle positions. An optional 4th column can
    !      contain the values of nu on the particles. If the 4th column is
    !      not present, nu will be roughly estimated using the density read
    !      from the ID and the average volume per particle; the latter is
    !      computed by computing the average particle separation along the
    !      z direction.

    PROCEDURE:: place_particles_lattice
    !! Places particles on a single lattice that surrounds both stars

    PROCEDURE:: place_particles_ellipsoidal_surfaces
    !! Places particles on ellipsoidal surfaces on one star

    PROCEDURE, NOPASS:: impose_equatorial_plane_symmetry
    !# Reflects the positions of the particles on a matter object with respect
    !  to the \(xy\) plane

    PROCEDURE:: get_object_of_particle
    !# Returns the number of the matter object asociated with the particle
    !  number given as input. Example: give number \(n\) as input; this
    !  particle number corresponds to a particle on matter object \(m\).
    !  This functions returns \(m\).

    PROCEDURE, NOPASS:: perform_apm
    !! Performs the Artificial Pressure Method (APM) on one star's particles

  !  GENERIC:: reshape_sph_field => reshape_sph_field_1d_ptr, &
  !                                 reshape_sph_field_2d_ptr
  !  !# GENERIC PROCEDURE, overloded to reallocate 1d and 2d arrays
  !  PROCEDURE:: reshape_sph_field_1d_ptr => reshape_sph_field_1d
  !  !! Reallocates a 1d array
  !  PROCEDURE:: reshape_sph_field_2d_ptr => reshape_sph_field_2d
  !  !! Reallocates a 2d array

    PROCEDURE:: allocate_particles_memory
    !! Allocates memory for the [[particles]] member arrays

    PROCEDURE:: deallocate_particles_memory
    !! Deallocates memory for the [[particles]] member arrays

    PROCEDURE:: compute_sph_hydro
    !# Computes the hydro fields on a section of the particles specified as
    !  input.
    !  First, computes the |sph| pressure starting from the |sph| baryon mass
    !  density, and the specific internal
    !  energy. The pressure is computed differently for different |eos|, and
    !  for cold and hot systems.
    !  Then computes the enthalpy and the sound speed accordingly.

    PROCEDURE:: read_compose_composition
    !! Reads the \(Y_e(n_b)\) table in the CompOSE file with extension .beta

    PROCEDURE:: compute_Ye
    !# Interpates linearly the electron fraction \(Y_e\) at the particle
    !  densities; that is, assigns \(Y_e\) at the particle positions

    PROCEDURE:: test_recovery
    !# Computes the conserved variables from the physical ones, and vice versa,
    !  to test that the recovered physical variables are the same to those
    !  computed from the |id|. @todo add reference for recovery

    !PROCEDURE:: compute_adm_momentum
    !# Computes an estimate of the \(\mathrm{ADM}\) linear momentum using
    !  the canonical momentum per baryon on the particles


    !
    !-- PUBLIC SUBROUTINES
    !

    PROCEDURE, PUBLIC:: analyze_hydro
    !# Scans the hydro fields taken from \(\texttt{|lorene|}\) to look
    !  for negative or zero values

    PROCEDURE, PUBLIC:: compute_and_print_sph_variables
    !# Computes the SPH variables at the particle positions, and optionally
    !  prints them to a binary file to be read by \(\texttt{SPHINCS_BSSN}\)
    !  and \(\texttt{splash}\), and to a formatted file to be read by
    !  \(\texttt{gnuplot}\), by calling
    !  [[particles:print_formatted_id_particles]]

    PROCEDURE, PUBLIC:: read_sphincs_dump_print_formatted
    !# Reads the binary ID file printed by
    !  [[particles:compute_and_print_sph_variables]], and prints it
    !  to a formatted file. @todo use this procedure in a second constructor

    PROCEDURE, PUBLIC:: print_formatted_id_particles
    !! Prints the |sph| |id| to a formatted file

    PROCEDURE, PUBLIC:: print_summary
    !! Prints the SPH ID to a formatted file

    PROCEDURE, PUBLIC:: is_empty
    !# Returns `.TRUE` if the [[particles]] object is empty, `.FALSE` otherwise
    !  @warning experimental, not actively used in the code yet

    !PROCEDURE, PUBLIC:: write_lorene_bns_id_dump

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

    PROCEDURE, PUBLIC:: get_npart
    !! Returns [[particles:npart]]
    PROCEDURE, PUBLIC:: get_npart_i
    !! Returns the number of particles on the object `i_matter`
    PROCEDURE, PUBLIC:: get_n_matter
    !! Returns [[particles:n_matter]]
    PROCEDURE, PUBLIC:: get_nuratio
    !! Returns [[particles:nuratio]]
    PROCEDURE, PUBLIC:: get_nuratio_i
    !! Returns the baryon number ratio on the object `i_matter`
    PROCEDURE, PUBLIC:: get_pos
    !! Returns [[particles:pos]]
    PROCEDURE, PUBLIC:: get_vel
    !! Returns [[particles:v]]
    PROCEDURE, PUBLIC:: get_nstar
    !! Returns [[particles:nstar]]
    PROCEDURE, PUBLIC:: get_nstar_sph
    !! Returns [[particles:nstar_sph]]
    PROCEDURE, PUBLIC:: get_nlrf
    !! Returns [[particles:nlrf]]
    PROCEDURE, PUBLIC:: get_nlrf_sph
    !! Returns [[particles:nlrf_sph]]
    PROCEDURE, PUBLIC:: get_nu
    !! Returns [[particles:nu]]
    PROCEDURE, PUBLIC:: get_u
    !! Returns [[particles:specific_energy]]
    PROCEDURE, PUBLIC:: get_u_sph
    !! Returns [[particles:u_sph]]
    PROCEDURE, PUBLIC:: get_pressure
    !! Returns [[particles:pressure]]
    PROCEDURE, PUBLIC:: get_pressure_sph
    !! Returns [[particles:pressure_sph]]
    PROCEDURE, PUBLIC:: get_theta
    !! Returns [[particles:theta]]
    PROCEDURE, PUBLIC:: get_h
    !! Returns [[particles:h]]
    PROCEDURE, PUBLIC:: get_lapse
    !! Returns [[particles:lapse]]
    PROCEDURE, PUBLIC:: get_shift
    !! Returns ([[particles:shift_x]],[[particles:shift_y]],[[particles:shift_z]])
    PROCEDURE, PUBLIC:: get_g3
    !! Returns ([[particles:g_xx]],[[particles:g_xy]],[[particles:g_xz]],[[particles:g_yy]],[[particles:g_yz]],[[particles:g_zz]])
    PROCEDURE, PUBLIC:: get_compose_eos
    !! Returns [[particles:compose_eos]]
    PROCEDURE, PUBLIC:: get_eos_id
    !! Returns the |eos| identifier for matter object `i_matter`

    FINAL:: destruct_particles
    !! Finalizer (Destructor) of [[particles]] object

  END TYPE particles