! File:         module_utility.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'.                                                            *
!************************************************************************

MODULE utility

  !***********************************************************************
  !
  !# This MODULE contains useful miscellaneous PROCEDURES and variables
  !
  !***********************************************************************


  USE matrix,     ONLY: determinant_4x4_matrix
  USE constants,  ONLY: MSun, G_Msun, c_light2


  IMPLICIT NONE


  INTEGER, PARAMETER:: flag$sph= -37457
  INTEGER, PARAMETER:: flag$tpo= -6543

  !
  !-- Identifiers for the supported equations of state
  !
  INTEGER, PARAMETER:: eos$poly        = 1
  INTEGER, PARAMETER:: eos$pwpoly      = 2
  INTEGER, PARAMETER:: eos$tabu$compose= 3


  DOUBLE PRECISION, PARAMETER:: zero        = 0.D0
  DOUBLE PRECISION, PARAMETER:: one         = 1.D0
  DOUBLE PRECISION, PARAMETER:: two         = 2.D0
  DOUBLE PRECISION, PARAMETER:: three       = 3.D0
  DOUBLE PRECISION, PARAMETER:: four        = 4.D0
  DOUBLE PRECISION, PARAMETER:: five        = 5.D0
  DOUBLE PRECISION, PARAMETER:: seven       = 7.D0
  DOUBLE PRECISION, PARAMETER:: ten         = 10.D0
  DOUBLE PRECISION, PARAMETER:: km2m        = ten*ten*ten
  DOUBLE PRECISION, PARAMETER:: g2kg        = one/(ten*ten*ten)
  DOUBLE PRECISION, PARAMETER:: MSun_geo    = G_Msun/c_light2/ &
                                               (ten*ten*ten*ten*ten)
  !# Msun_geo = 1.47662503825040 km
  !  see https://einsteintoolkit.org/thornguide/EinsteinBase/HydroBase/documentation.html
  DOUBLE PRECISION, PARAMETER:: MeV2amuc2    = one/931.49432
  !! Conversion factor from \(\mathrm{MeV}\) to \(\mathrm{amu}*c^2\)
  DOUBLE PRECISION, PARAMETER:: km2Msun_geo  = one/MSun_geo
  DOUBLE PRECISION, PARAMETER:: density_si2cu= (MSun_geo*km2m)**3/(MSun*g2kg)
  !# Conversion factor for the baryon mass density, from SI units to code
  !  units, but NOT measured in units of \(m_0c^2\)
  !
  !  `dens_si2cu`\(=\simeq\dfrac{(1477\mathrm{m})^3}{2*10^30\mathrm{kg}}=1.6186541582311746851140226630074\times 10^{-21}\dfrac{\mathrm{m}^3}{\mathrm{kg}}\)

  INTEGER:: itr
  !! Iterator for loops
  INTEGER:: itr3
  !! Iterator for loops
  INTEGER:: itr4
  !! Iterator for loops
  INTEGER:: ios
  !! Variable to store the state of I/O
  INTEGER:: cnt= 0
  !! Counter

  !
  !-- Variables to print progress on screen
  !

  INTEGER:: perc
  !! Percentage of execution time (integer)
  DOUBLE PRECISION:: perc2
  !! Percentage of execution time (double)
  CHARACTER, PARAMETER:: creturn= ACHAR(13)
  !! Carriage return

  LOGICAL:: file_exists
  !! `.TRUE.` if a given file exists; `.FALSE.` otherwise
  LOGICAL:: show_progress
  !# `.TRUE.` if loop progress is to be printed to standard output;
  !  `.FALSE.` otherwise

  CHARACTER(LEN=:), ALLOCATABLE:: err_msg
  !! String storing error messages

  CHARACTER(LEN=500):: hostname
  !# String storing the name of the host machine

  CHARACTER(LEN=10):: version
  !# String storing the version of |sphincsid|

  !
  !-- Variables used to set the run_id
  !

  CHARACTER(8)  :: date
  !! Date when the run starts
  CHARACTER(10) :: time
  !! Time when the run starts
  CHARACTER(5)  :: zone
  !! Place where the run runs
  INTEGER, DIMENSION(8) :: values
  !# An integer array of 8 elements described below:
  !
  !   1. The year as a 4-digit integer
  !   2. The month as an integer from 1 to 12
  !   3. The day of the moneth as an integer from 1 to 31
  !   4. The time difference, in minutes, with respect to UTC
  !      (Coordinated Universal Time)
  !   5. The hour of the day as an integer from 1 to 23
  !   6. The minutes of the hour as an integer from 1 to 59
  !   7. The second of the minute as an integer from 0 to 60
  !   8. The millisecond of the second as an integer from 0 to 999
  CHARACTER(LEN=19):: run_id
  !! Identification string for the run
  CHARACTER(LEN=19):: end_time
  !! Time when the run ends


  INTEGER, PARAMETER:: max_length= 100
  !! Maximum length for strings
  INTEGER, PARAMETER:: max_n_id= 50
  ! Maximum number of physical systems
  INTEGER, PARAMETER:: max_n_parts= 250
  !! Maximum number of particle distributions

  INTEGER:: n_id
  !! Number of physical systems to set up
  INTEGER:: ref_lev
  !! Number of refinement levels
  INTEGER:: constraints_step
  !! Export the constraints every constraints_step-th step

  INTEGER, PARAMETER:: test_int= - 112
  INTEGER, DIMENSION(max_n_id, max_n_parts):: placer= test_int
  !# Matrix storing the information on how to place particles for each bns
  !  object. Row i contains information about the i^th bns object.


  DOUBLE PRECISION:: numerator_ratio_dx
  !# Numerator of the rational ratio between the large grid spacing and the
  !  medium one,equal to the ratio between the medium grid spacing nd the small
  !  one. Not used in this PROGRAM, but needed since the PROGRAM reads the same
  !  parameter file as the convergence_test PROGRAM
  DOUBLE PRECISION:: denominator_ratio_dx
  !# Denominator of the rational ratio between the large grid spacing and the
  !  medium one,equal to the ratio between the medium grid spacing nd the small
  !  one. Not used in this PROGRAM, but needed since the PROGRAM reads the same
  !  parameter file as the convergence_test PROGRAM

  ! Logical variables to steer the execution
  LOGICAL:: export_bin, export_form, export_form_xy, export_form_x, &
            use_eos_from_id, &
            compute_constraints, export_constraints_xy, &
            export_constraints_x, export_constraints, &
            export_constraints_details, compute_parts_constraints, &
            one_lapse, zero_shift, run_sph, run_spacetime, estimate_length_scale

  CHARACTER(LEN=max_length), DIMENSION(max_length):: filenames= "0"
  !! Array of strings storing the names of the |id| files
  CHARACTER(LEN=max_length), DIMENSION(50):: eos_filenames= "use_id"
  !! Array of strings storing the names of the files containing the |eos|
  CHARACTER(LEN=max_length):: common_path
  !# String storing the path to the directory where the |id| files
  !  are stored. Local and global paths can be used
  CHARACTER(LEN=max_length):: common_eos_path
  !# String storing the path to the directory where the |eos| files
  !  are stored. Local and global paths can be used
  CHARACTER(LEN=max_length):: sph_path
  !# String storing the local path to the directory where the
  !  SPH output is to be saved
  CHARACTER(LEN=max_length):: spacetime_path
  !# String storing the local path to the directory where the
  !  spacetime output is to be saved


  TYPE surface
    LOGICAL:: is_known= .FALSE.
    DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE:: points
    !# Array containing the coordinates of the matter objects' surfaces
    !  The first index runs over the stars; the second and third run over
    !  the surface points (number of points for \(\theta\) and \(\varphi\));
    !  the fourth index runs over the Cartesian coordinates of the points
  END TYPE surface


  TYPE tabu_eos
    DOUBLE PRECISION, PUBLIC, DIMENSION(:,:), ALLOCATABLE:: table_eos
    !# Array containing the tabulated |eos|.
    !  Its size is determined when reading the tabulated |eos| from file
    !  The first index runs over the field, the second over its values.
  END TYPE tabu_eos


  CONTAINS


  SUBROUTINE read_sphincs_id_parameters()

    !***********************************************
    !
    !# Read the parameters to steer SPHINCS_ID
    !
    !  FT
    !
    !***********************************************

    IMPLICIT NONE

    INTEGER:: stat
    INTEGER, PARAMETER:: unit_parameters= 17

    CHARACTER(LEN=:), ALLOCATABLE:: sphincs_id_parameters_namefile
    CHARACTER(LEN=100):: msg

    ! Namelist containing parameters read from sphincs_id_parameters.par
    ! by the SUBROUTINE read_sphincs_id_parameters of this PROGRAM
    NAMELIST /sphincs_id_parameters/ &
              n_id, common_path, filenames, &
              common_eos_path, use_eos_from_id, eos_filenames, &
              placer, export_bin, export_form, export_form_xy, &
              export_form_x, export_constraints_xy, &
              export_constraints_x, compute_constraints, &
              export_constraints, export_constraints_details, &
              constraints_step, compute_parts_constraints, &
              numerator_ratio_dx, denominator_ratio_dx, ref_lev, &
              one_lapse, zero_shift, show_progress, &
              run_sph, run_spacetime, sph_path, spacetime_path, &
              estimate_length_scale

    use_eos_from_id= .TRUE.

    sphincs_id_parameters_namefile= 'sphincs_id_parameters.dat'

    INQUIRE( FILE= sphincs_id_parameters_namefile, EXIST= file_exists )
    IF( file_exists )THEN

      OPEN( unit_parameters, FILE= sphincs_id_parameters_namefile, &
            STATUS= 'OLD' )

    ELSE

      PRINT*
      PRINT*,'** ERROR: ', sphincs_id_parameters_namefile, " file not found!"
      PRINT*
      STOP

    ENDIF

    READ( UNIT= unit_parameters, NML= sphincs_id_parameters, IOSTAT= stat, &
          IOMSG= msg )

    IF( stat /= 0 )THEN
      PRINT *, "** ERROR: Error in reading ", sphincs_id_parameters_namefile, &
               ". The IOSTAT variable is ", stat, &
               "The error message is", msg
      STOP
    ENDIF

    CLOSE( UNIT= unit_parameters )

    DO itr= 1, max_length, 1
      IF( TRIM(filenames(itr)).NE."0" )THEN
        cnt= cnt + 1
      ENDIF
    ENDDO
    IF( cnt.NE.n_id )THEN
      PRINT *, "** ERROR! The number of file names is", cnt, &
               "and n_id=", n_id, ". The two should be the same."
      PRINT *
      STOP
    ENDIF

   !DO itr= 1, n_id, 1
   !  DO itr2= 1, max_n_parts, 1
   !    IF( placer( itr, itr2 ) == test_int )THEN
   !      PRINT *
   !      PRINT *, "** ERROR! The array placer does not have ", &
   !               "enough components to specify all the desired ", &
   !               "particle distributions. Specify the ", &
   !               "components in file sphincs_id_particles.par"
   !      PRINT *
   !      STOP
   !    ENDIF
   !  ENDDO
   !ENDDO

  END SUBROUTINE read_sphincs_id_parameters


  SUBROUTINE test_status( io_stat, io_msg, opt_msg )

    !***********************************************
    !
    !# Test if a status variable is 0 or not
    !
    !  FT 17.09.2020
    !
    !***********************************************

    IMPLICIT NONE

    INTEGER,            INTENT(IN)           :: io_stat
    !! Status variable
    CHARACTER(LEN=100), INTENT(IN)           :: io_msg
    !! Status message
    CHARACTER(LEN=*),   INTENT(IN), OPTIONAL :: opt_msg
    !! Optional status message

    IF( io_stat > 0 )THEN

      PRINT *
      PRINT *, "***** ERROR! IOSTAT > 0. ", &
               "The error message is: ", io_msg
      IF( PRESENT( opt_msg ) )THEN
        PRINT *, opt_msg
      ENDIF
      PRINT *
      STOP

    ENDIF

  END SUBROUTINE test_status


  PURE FUNCTION is_finite_number( x ) RESULT( res )

    !***********************************************
    !
    !# Test if a double precision is a finite number
    !
    !  FT 11.02.2022
    !
    !***********************************************

    USE, INTRINSIC:: IEEE_ARITHMETIC, ONLY: IEEE_IS_FINITE

    DOUBLE PRECISION, INTENT(IN):: x
    LOGICAL:: res

    res= (.NOT.ISNAN(x)) .AND. IEEE_IS_FINITE(x)

  END FUNCTION is_finite_number


  SUBROUTINE scan_1d_array_for_nans( array_size, array, name )

    !***********************************************
    !
    !# Test if a double precision is a finite number
    !
    !  FT 11.02.2022
    !
    !***********************************************

    INTEGER,                        INTENT(IN):: array_size
    DOUBLE PRECISION, DIMENSION(:), INTENT(IN):: array
    CHARACTER(LEN=*),               INTENT(IN):: name

    INTEGER:: i

    !$OMP PARALLEL DO DEFAULT( NONE ) &
    !$OMP             SHARED( array, array_size, name ) &
    !$OMP             PRIVATE( i )
    DO i= 1, array_size, 1

      IF(.NOT.is_finite_number(array(i)))THEN
        PRINT *, "** ERROR! The array ", name, "does not contain a finite", &
                 " number at position ", i
        PRINT *, " * ", name, "(", i, ")=", array(i)
        PRINT *
        STOP
      ENDIF

    ENDDO
    !$OMP END PARALLEL DO

  END SUBROUTINE scan_1d_array_for_nans


  SUBROUTINE scan_3d_array_for_nans( size_1, size_2, size_3, array, name )

    !***********************************************
    !
    !# Test if a double precision is a finite number
    !
    !  FT 11.02.2022
    !
    !***********************************************

    INTEGER,                            INTENT(IN):: size_1
    INTEGER,                            INTENT(IN):: size_2
    INTEGER,                            INTENT(IN):: size_3
    DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN):: array
    CHARACTER(LEN=*),                   INTENT(IN):: name

    INTEGER:: i, j, k

    !$OMP PARALLEL DO DEFAULT( NONE ) &
    !$OMP             SHARED( array, size_1, size_2, size_3, name ) &
    !$OMP             PRIVATE( i, j, k )
    DO k= 1, size_3, 1
      DO j= 1, size_2, 1
        DO i= 1, size_1, 1

          IF(.NOT.is_finite_number(array(i,j,k)))THEN
            PRINT *, "** ERROR! The array ", name, "does not contain a ", &
                     "finite number at position ", "(", i, ",", j, ",", k, ")"
            PRINT *, " * ", name, "(", i, ",", j, ",", k, ")=", array(i,j,k)
            PRINT *
            STOP
          ENDIF

        ENDDO
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO

  END SUBROUTINE scan_3d_array_for_nans


  PURE SUBROUTINE compute_g4( lapse, shift, g3, g4 )

    !***********************************************
    !
    !# Computes the spacetime metric from lapse,
    !  shift and spatial metric
    !
    !  FT 27.11.2020
    !
    !  Generalized to not be bound to the mesh
    !
    !  FT 07.02.2022
    !
    !***********************************************

    USE tensor,     ONLY: itt, itx, ity, itz, ixx, ixy, &
                          ixz, iyy, iyz, izz, jxx, jxy, jxz, &
                          jyy, jyz, jzz, jx, jy, jz

    IMPLICIT NONE

    DOUBLE PRECISION,                INTENT(IN)   :: lapse
    !! Lapse function
    DOUBLE PRECISION, DIMENSION(3),  INTENT(IN)   :: shift
    !! Contravariant shift vector
    DOUBLE PRECISION, DIMENSION(6),  INTENT(IN)   :: g3
    !! Covariant spatial metric
    DOUBLE PRECISION, DIMENSION(10), INTENT(INOUT):: g4
    !! Covariant spacetime metric

    g4(itt)= - lapse*lapse + g3(jxx)*shift(jx)*shift(jx)     &
                           + g3(jxy)*shift(jx)*shift(jy)*two &
                           + g3(jxz)*shift(jx)*shift(jz)*two &
                           + g3(jyy)*shift(jy)*shift(jy)     &
                           + g3(jyz)*shift(jy)*shift(jz)*two &
                           + g3(jzz)*shift(jz)*shift(jz)

    g4(itx)= g3(jxx)*shift(jx) + g3(jxy)*shift(jy) + g3(jxz)*shift(jz)
    g4(ity)= g3(jxy)*shift(jx) + g3(jyy)*shift(jy) + g3(jyz)*shift(jz)
    g4(itz)= g3(jxz)*shift(jx) + g3(jyz)*shift(jy) + g3(jzz)*shift(jz)

    g4(ixx)= g3(jxx)
    g4(ixy)= g3(jxy)
    g4(ixz)= g3(jxz)
    g4(iyy)= g3(jyy)
    g4(iyz)= g3(jyz)
    g4(izz)= g3(jzz)

  END SUBROUTINE compute_g4


  SUBROUTINE compute_tpo_metric( g4, lapse, shift, g3 )

    !***********************************************
    !
    !# Computes the lapse,shift and spatial metric
    !  from the covariant spacetime metric
    !
    !  FT 12.04.2022
    !
    !***********************************************

    USE matrix,     ONLY: invert_3x3_matrix
    USE tensor,     ONLY: itt, itx, ity, itz, ixx, ixy, &
                          ixz, iyy, iyz, izz, jxx, jxy, jxz, &
                          jyy, jyz, jzz, jx, jy, jz, n_sym4x4, n_sym3x3

    IMPLICIT NONE

    DOUBLE PRECISION, DIMENSION(n_sym4x4),  INTENT(IN) :: g4
    !! Covariant spacetime metric
    DOUBLE PRECISION,                       INTENT(OUT):: lapse
    !! Lapse function
    DOUBLE PRECISION, DIMENSION(3),         INTENT(OUT):: shift
    !! Contravariant shift vector
    DOUBLE PRECISION, DIMENSION(n_sym3x3),  INTENT(OUT):: g3
    !! Covariant spatial metric

    DOUBLE PRECISION, DIMENSION(3,3):: gmat
    DOUBLE PRECISION, DIMENSION(3,3):: gmat_inv

    g3(jxx)= g4(ixx)
    g3(jxy)= g4(ixy)
    g3(jxz)= g4(ixz)
    g3(jyy)= g4(iyy)
    g3(jyz)= g4(iyz)
    g3(jzz)= g4(izz)

    gmat(1,1)= g3(jxx)
    gmat(1,2)= g3(jxy)
    gmat(1,3)= g3(jxz)
    gmat(2,1)= g3(jxy)
    gmat(2,2)= g3(jyy)
    gmat(2,3)= g3(jyz)
    gmat(3,1)= g3(jxz)
    gmat(3,2)= g3(jyz)
    gmat(3,3)= g3(jzz)

    CALL invert_3x3_matrix( gmat, gmat_inv )

    shift(jx)= gmat_inv(jx,jx)*g4(itx) + gmat_inv(jx,jy)*g4(ity) &
             + gmat_inv(jx,jz)*g4(itz)
    shift(jy)= gmat_inv(jy,jx)*g4(itx) + gmat_inv(jy,jy)*g4(ity) &
             + gmat_inv(jy,jz)*g4(itz)
    shift(jz)= gmat_inv(jz,jx)*g4(itx) + gmat_inv(jz,jy)*g4(ity) &
             + gmat_inv(jz,jz)*g4(itz)

    lapse= SQRT( shift(jx)*g4(itx) + shift(jy)*g4(ity) + shift(jz)*g4(itz) &
                 - g4(itt) )

  END SUBROUTINE compute_tpo_metric


  SUBROUTINE invert_sym4x4( A, iA )

    !****************************************************************
    !
    !# Invert a \(4\times 4\) symemtric matrix stored as a \(10\)-vector
    !  @note The inverse (and the adjugate) of a symmetric matrix
    !        is symmetric
    !
    !  FT 25.04.2022
    !
    !****************************************************************

    USE tensor,               ONLY: n_sym4x4
    !USE metric_on_particles,  ONLY: gvec2mat, mat2gvec
    USE matrix,               ONLY: invert_4x4_matrix

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: A(:)
    !# The \(4\times 4\) symmetric matrix, given as a 10-vector.
    !  The first 3 components run over the numbers of grid points
    !  along each axis. The fourth index runs over the number of
    !  independent components of the \(4\times 4\) symmetric matrix.
    DOUBLE PRECISION, INTENT(OUT):: iA(n_sym4x4)
    !# Inverse of the \(4\times 4\) symmetric matrix, given as input.

    DOUBLE PRECISION:: Amat(4,4)
    !! The \(4\times 4\) symmetric matrix as a matrix
    DOUBLE PRECISION:: iAmat(4,4)
    !! The inverse of the \(4\times 4\) symmetric matrix, as a matrix

    IF( SIZE(A) /= n_sym4x4 )THEN
      PRINT *, "** ERROR in determinant_sym4x4_grid in MODULE utility.", &
               " This subroutine needs a symmetric matrix with 10 components,",&
               " and a ", SIZE(A), "component matrix was given instead."
      STOP
    ENDIF

    !CALL gvec2mat(A,Amat)
    !Amat= [[A(1),A(2),A(3),A(4)], &
    !       [A(2),A(5),A(6),A(7)], &
    !       [A(3),A(6),A(8),A(9)], &
    !       [A(4),A(7),A(9),A(10)]]
    Amat(1,1)= A(1)
    Amat(1,2)= A(2)
    Amat(1,3)= A(3)
    Amat(1,4)= A(4)
    Amat(2,1)= A(2)
    Amat(2,2)= A(5)
    Amat(2,3)= A(6)
    Amat(2,4)= A(7)
    Amat(3,1)= A(3)
    Amat(3,2)= A(6)
    Amat(3,3)= A(8)
    Amat(3,4)= A(9)
    Amat(4,1)= A(4)
    Amat(4,2)= A(7)
    Amat(4,3)= A(9)
    Amat(4,4)= A(10)

    CALL invert_4x4_matrix(Amat,iAmat)

    !CALL mat2gvec(iA,iAmat)
    iA= [iAmat(1,1),iAmat(1,2),iAmat(1,3),iAmat(1,4), &
         iAmat(2,2),iAmat(2,3),iAmat(2,4), &
         iAmat(3,3),iAmat(3,4), &
         iAmat(4,4)]

  END SUBROUTINE invert_sym4x4


  SUBROUTINE invert_sym3x3( A, iA )

    !****************************************************************
    !
    !# Invert a \(3\times 3\) symemtric matrix stored as a \(6\)-vector
    !  @note The inverse (and the adjugate) of a symmetric matrix
    !        is symmetric
    !
    !  FT 25.04.2022
    !
    !****************************************************************

    USE tensor, ONLY: n_sym3x3
    USE matrix, ONLY: invert_3x3_matrix

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: A(:)
    !# The \(3\times 3\) symmetric matrix, given as a \(6\)-vector.
    DOUBLE PRECISION, INTENT(OUT):: iA(n_sym3x3)
    !# Inverse of the \(3\times 3\) symmetric matrix, given as input.

    DOUBLE PRECISION:: Amat(3,3)
    !! The \(3\times 3\) symmetric matrix as a matrix
    DOUBLE PRECISION:: iAmat(3,3)
    !! The inverse of the \(3\times 3\) symmetric matrix, as a matrix

    IF( SIZE(A) /= n_sym3x3 )THEN
      PRINT *, "** ERROR in determinant_sym3x3_grid in MODULE utility.", &
               " This subroutine needs a symmetric matrix with 6 components,",&
               " and a ", SIZE(A), "component matrix was given instead."
      STOP
    ENDIF

    CALL vec2mat_sym3x3(A,Amat)

    CALL invert_3x3_matrix(Amat,iAmat)

    CALL mat2vec_sym3x3(iAmat,iA)

  END SUBROUTINE invert_sym3x3


  SUBROUTINE determinant_sym4x4( A, det )

    !****************************************************************
    !
    !# Compute the determinant of a \(4\times 4\) symmetric matrix
    !  field, given as a 10-vector
    !
    !  FT 26.01.2022
    !
    !  Generalized to not be bound to the mesh
    !
    !  FT 07.02.2022
    !
    !****************************************************************

    USE tensor, ONLY: itt, itx, ity, itz, ixx, ixy, ixz, iyy, iyz, izz, n_sym4x4

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: A(:)
    !# The \(4\times 4\) symmetric matrix, given as a 10-vector.
    !  The first 3 components run over the numbers of grid points
    !  along each axis. The fourth index runs over the number of
    !  independent components of the \(4\times 4\) symmetric matrix.
    DOUBLE PRECISION, INTENT(OUT):: det
    !! Determinant of the \(4\times 4\) symmetric matrix

    IF( SIZE(A) /= n_sym4x4 )THEN
      PRINT *, "** ERROR in determinant_sym4x4_grid in MODULE utility.", &
               " This subroutine needs a symmetric matrix with 10 components,",&
               " and a ", SIZE(A), "component matrix was given instead."
      STOP
    ENDIF

    det=   A(itt)*(A(ixx)*(A(iyy)*A(izz) &
         - A(iyz)* A(iyz)) &
         + A(ixy)*(A(iyz)* A(ixz) &
         - A(ixy)* A(izz)) &
         + A(ixz)*(A(ixy)* A(iyz) &
         - A(iyy)* A(ixz))) &
         - A(itx)*(A(itx)*(A(iyy)*A(izz) &
         - A(iyz)* A(iyz)) &
         + A(ixy)*(A(iyz)* A(itz) &
         - A(ity)* A(izz)) &
         + A(ixz)*(A(ity)* A(iyz) &
         - A(iyy)* A(itz))) &
         + A(ity)*(A(itx)*(A(ixy)*A(izz) &
         - A(iyz)* A(ixz)) &
         + A(ixx)*(A(iyz)* A(itz) &
         - A(ity)* A(izz)) &
         + A(ixz)*(A(ity)* A(ixz) &
         - A(ixy)* A(itz))) &
         - A(itz)*(A(itx)*(A(ixy)*A(iyz) &
         - A(iyy)* A(ixz)) &
         + A(ixx)*(A(iyy)* A(itz) &
         - A(ity)* A(iyz)) &
         + A(ixy)*(A(ity)* A(ixz) &
         - A(ixy)* A(itz)))

  END SUBROUTINE determinant_sym4x4


  SUBROUTINE determinant_sym3x3( A, det )

    !****************************************************************
    !
    !# Compute the determinant of a \(3\times 3\) symmetric matrix
    !  field, given as a 6-vector, at a given grid point
    !
    !  FT 26.03.2021
    !
    !  Generalized to not be bound to the mesh
    !
    !  FT 07.02.2022
    !
    !****************************************************************

    USE tensor, ONLY: jxx, jxy, jxz, jyy, jyz, jzz, n_sym3x3

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: A(jxx:jzz)
    !# The \(3\times 3\) symmetric matrix, given as a 6-vector.
    DOUBLE PRECISION, INTENT(OUT):: det
    !! Determinant of the \(3\times 3\) symmetric matrix

    IF( SIZE(A) /= n_sym3x3 )THEN
      PRINT *, "** ERROR in determinant_sym3x3_grid in MODULE utility.", &
               " This subroutine needs an array with 6 components,",&
               " and a ", SIZE(A), "component array was given instead."
      STOP
    ENDIF

    det=   A(jxx)*A(jyy)*A(jzz) &
         + A(jxy)*A(jyz)*A(jxz) &
         + A(jxz)*A(jxy)*A(jyz) &
         - A(jxy)*A(jxy)*A(jzz) &
         - A(jxz)*A(jyy)*A(jxz) &
         - A(jxx)*A(jyz)*A(jyz)

  END SUBROUTINE determinant_sym3x3


  SUBROUTINE spacetime_vector_norm_sym4x4( g4, v, norm )

    !****************************************************************
    !
    !# Compute the spacetime squared norm of a vector, using the
    !  metric given as an array of 10 components
    !
    !  FT 07.02.2022
    !
    !****************************************************************

    USE tensor,    ONLY: itt, itx, ity, itz, ixx, ixy, ixz, iyy, iyz, izz, &
                         it, ix, iy, iz, n_sym4x4

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: g4(itt:izz)
    !# The \(4\times 4\) spacetime metric, given as a 10-vector.
    DOUBLE PRECISION, INTENT(IN):: v(it:iz)
    !# The \(4\)-vector whose norm has to be computed.
    DOUBLE PRECISION, INTENT(OUT):: norm
    !! Spacetime norm of the vector v.

    IF( SIZE(g4) /= n_sym4x4 )THEN
      PRINT *, "** ERROR in spacetime_vector_norm_sym4x4 in MODULE utility.", &
               " This subroutine needs a symmetric matrix with 10 components,",&
               " and a ", SIZE(g4), "component matrix was given instead."
      STOP
    ENDIF

    norm= g4(itt)*v(it)*v(it)     + two*g4(itx)*v(it)*v(ix) &
        + two*g4(ity)*v(it)*v(iy) + two*g4(itz)*v(it)*v(iz) &
        + g4(ixx)*v(ix)*v(ix)     + two*g4(ixy)*v(ix)*v(iy) &
        + two*g4(ixz)*v(ix)*v(iz) + g4(iyy)*v(iy)*v(iy)     &
        + two*g4(iyz)*v(iy)*v(iz) + g4(izz)*v(iz)*v(iz)

  END SUBROUTINE spacetime_vector_norm_sym4x4


  SUBROUTINE spatial_vector_norm_sym3x3( g3, v, norm )

    !****************************************************************
    !
    !# Compute the spatial squared norm of a vector, using the
    !  spatial metric given as an array of 6 components
    !
    !  FT 14.02.2022
    !
    !****************************************************************

    USE tensor,    ONLY: jxx, jxy, jxz, jyy, jyz, jzz, jx, jy, jz, n_sym3x3

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: g3(jxx:jzz)
    !# The \(3\times 3\) spacetime metric, given as a 6-vector.
    DOUBLE PRECISION, INTENT(IN):: v(jx:jz)
    !# The \(3\)-vector whose norm has to be computed.
    DOUBLE PRECISION, INTENT(OUT):: norm
    !! Spatial norm of the vector v.

    IF( SIZE(g3) /= n_sym3x3 )THEN
      PRINT *, "** ERROR in spatial_vector_norm_sym3x3 in MODULE utility.", &
               " This subroutine needs a symmetric matrix with 6 components,",&
               " and a ", SIZE(g3), "component matrix was given instead."
      STOP
    ENDIF

    norm= g3(jxx)*v(jx)*v(jx)     + two*g3(jxy)*v(jx)*v(jy) &
        + two*g3(jxz)*v(jx)*v(jz) + g3(jyy)*v(jy)*v(jy)     &
        + two*g3(jyz)*v(jy)*v(jz) + g3(jzz)*v(jz)*v(jz)

  END SUBROUTINE spatial_vector_norm_sym3x3


  PURE SUBROUTINE vec2mat_sym3x3( vec, mat )

    !********************************************
    !
    !# Write the components of symmetric \(3\times 3\)
    !  matrix given as a \(6\)-vector, into a
    !  \(3\times 3\) matrix
    !
    ! FT 25.04.2022
    !
    !*********************************************

    USE tensor,    ONLY: jxx, jxy, jxz, jyy, jyz, jzz, jx, jy, jz, n_sym3x3

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN)  :: vec(n_sym3x3)
    DOUBLE PRECISION, INTENT(OUT) :: mat(3,3)

    mat(1,1)= vec(jxx)
    mat(1,2)= vec(jxy)
    mat(1,3)= vec(jxz)

    mat(2,1)= mat(1,2)
    mat(2,2)= vec(jyy)
    mat(2,3)= vec(jyz)

    mat(3,1)= mat(1,3)
    mat(3,2)= mat(2,3)
    mat(3,3)= vec(jzz)

  END SUBROUTINE vec2mat_sym3x3


  PURE SUBROUTINE mat2vec_sym3x3( mat, vec )

    !************************************
    !                                   *
    ! transform symmetric 4x4-matrix    *
    ! into vector; SKR 30.11.2017       *
    !                                   *
    !************************************

    USE tensor,    ONLY: jxx, jxy, jxz, jyy, jyz, jzz, jx, jy, jz, n_sym3x3

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN)  :: mat(3,3)
    DOUBLE PRECISION, INTENT(OUT) :: vec(n_sym3x3)

    vec(jxx)= mat(1,1)
    vec(jxy)= mat(1,2)
    vec(jxz)= mat(1,3)
    vec(jyy)= mat(2,2)
    vec(jyz)= mat(2,3)
    vec(jzz)= mat(3,3)

  END SUBROUTINE mat2vec_sym3x3


  PURE SUBROUTINE spherical_from_cartesian( x, y, z, xo, yo, zo, r, theta, phi )

    !****************************************************************
    !
    !# Compute the spherical polar coordinates of a point \(p\)
    !  relative to a point \(O\), starting from the Cartesian
    !  coordinates of the points \(p\) and \(O\)
    !
    !  FT 14.01.2021
    !
    !****************************************************************

    USE constants, ONLY: pi

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: x
    !! \(x\) coordinate of the point \(p\)
    DOUBLE PRECISION, INTENT(IN):: y
    !! \(y\) coordinate of the point \(p\)
    DOUBLE PRECISION, INTENT(IN):: z
    !! \(z\) coordinate of the point \(p\)

    DOUBLE PRECISION, INTENT(IN):: xo
    !! \(x\) coordinate of the point \(O\)
    DOUBLE PRECISION, INTENT(IN):: yo
    !! \(y\) coordinate of the point \(O\)
    DOUBLE PRECISION, INTENT(IN):: zo
    !! \(z\) coordinate of the point \(O\)

    DOUBLE PRECISION, INTENT(OUT):: r
    !! \(r\) coordinate of the point \(p\), relative to \(O\)
    DOUBLE PRECISION, INTENT(OUT):: theta
    !! \(\theta\) coordinate (colatitude) of the point \(p\), relative to \(O\)
    DOUBLE PRECISION, INTENT(OUT):: phi
    !! \(\phi\) coordinate (azimuth) of the point \(p\), relative to \(O\)

    DOUBLE PRECISION:: xd, yd, zd

    xd= x - xo
    yd= y - yo
    zd= z - zo

    IF( xd > zero )THEN

      phi= ATAN( yd/xd )

    ELSEIF( xd < zero )THEN

      phi= ATAN( yd/xd ) + pi

    ELSE

      phi= pi/two

    ENDIF

    DO WHILE(phi < zero)

      phi= phi + two*pi

    ENDDO

    r= SQRT( xd**2 + yd**2 + zd**2 )

    theta= ACOS( zd/r )

  END SUBROUTINE spherical_from_cartesian


  PURE SUBROUTINE cartesian_from_spherical &
    ( r, theta, phi, xo, yo, zo, x, y, z, a_y, a_z )

    !****************************************************************
    !
    !# Compute the Cartesian coordinates of a points \(p\),
    !  starting from the spherical, or optionally elliptical, polar
    !  coordinates of the point \(p\) relative to a point \(O\)
    !
    !  FT 28.11.2022
    !
    !****************************************************************

    USE constants, ONLY: pi

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN):: r
    !! \(r\) coordinate of the point \(p\), relative to \(O\)
    DOUBLE PRECISION, INTENT(IN):: theta
    !! \(\theta\) coordinate (colatitude) of the point \(p\), relative to \(O\)
    DOUBLE PRECISION, INTENT(IN):: phi
    !! \(\phi\) coordinate (azimuth) of the point \(p\), relative to \(O\)

    DOUBLE PRECISION, INTENT(IN):: xo
    !! \(x\) coordinate of the point \(O\)
    DOUBLE PRECISION, INTENT(IN):: yo
    !! \(y\) coordinate of the point \(O\)
    DOUBLE PRECISION, INTENT(IN):: zo
    !! \(z\) coordinate of the point \(O\)

    DOUBLE PRECISION, INTENT(OUT):: x
    !! \(x\) coordinate of the point \(p\)
    DOUBLE PRECISION, INTENT(OUT):: y
    !! \(y\) coordinate of the point \(p\)
    DOUBLE PRECISION, INTENT(OUT):: z
    !! \(z\) coordinate of the point \(p\)

    DOUBLE PRECISION, INTENT(IN), OPTIONAL:: a_y
    !! Ratio between the y and x semiaxes of the ellipse
    DOUBLE PRECISION, INTENT(IN), OPTIONAL:: a_z
    !! Ratio between the z and x semiaxes of the ellipse

    DOUBLE PRECISION:: ay, az

    IF(PRESENT(a_y))THEN
      ay= a_y
    ELSE
      ay= one
    ENDIF
    IF(PRESENT(a_z))THEN
      az= a_z
    ELSE
      az= one
    ENDIF

    x= xo +    r*SIN(theta)*COS(phi)
    y= yo + ay*r*SIN(theta)*SIN(phi)
    z= zo + az*r*COS(theta)

  END SUBROUTINE cartesian_from_spherical


  PURE FUNCTION k_lorene2cu( gam )

    !****************************************************************
    !
    !# Compute the constant to convert the polytropic constant \(K\)
    !  from the units used in |lorene| for the single polytropic |eos|,
    !  to the units used in |sphincs|
    !
    !  FT xx.xx.2020
    !
    !****************************************************************

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN) :: gam
    !! Polytropic exponent \(\gamma\)
    DOUBLE PRECISION :: k_lorene2cu

    ! LORENE's EOS is in terms on number density n = rho/m_nucleon:
    ! P = K n^Gamma
    ! to convert to SI units:
    ! K_SI(n) = K_LORENE rho_nuc c^2 / n_nuc^gamma
    ! Converting this to be in terms of the mass density rho = n m_nucleon gets
    ! changes n_nuc to rho_nuc:
    ! K_SI(rho) = K_LORENE c^2 / rho_nuc^(gamma-1)
    ! In SI units P has units of M / (L T^2) and rho has units of M/L^3 thus
    ! K_SI has units of (L^3/M)^Gamma M/(L T^2).
    ! In Cactus units P and rho have the same units thus K_Cactus is unitless.
    ! Conversion between K_SI and K_Cactus thus amounts to dividing out the
    ! units of the SI quantity.

    k_lorene2cu= ( (MSun*g2kg)/((MSun_geo*km2m)**3*(1.66D+17)) )**(gam - one)


  END FUNCTION k_lorene2cu


  PURE FUNCTION k_lorene2cu_pwp( gamma0 )

    !****************************************************************
    !
    !# Compute the constant to convert the polytropic constant \(K\)
    !  from the units used in |lorene| for the piecewise polytropic
    !  |eos|, to the units used in |sphincs|
    !
    !  FT xx.xx.2020
    !
    !****************************************************************

    USE constants,  ONLY: m2cm

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN) :: gamma0
    !! Polytropic exponent \(\gamma_0\)
    DOUBLE PRECISION :: k_lorene2cu_pwp

    ! LORENE has K0 in units of (g cm^{-3})^{1-gamma0} for the piecewise
    ! polytropes. This factor writes it in SPHINCS units

    k_lorene2cu_pwp= ( MSun/((MSun_geo*km2m*m2cm)**3) )**(gamma0 - one)


  END FUNCTION k_lorene2cu_pwp


  FUNCTION shorten_eos_name( eos_long ) RESULT( eos_str )

    !*********************************************
    !                                            *
    ! Converts the longer names from LORENE      *
    ! to the 4-character strings needed by       *
    ! select_EOS_parameters                      *
    !                                            *
    ! FT 11.08.2021                              *
    !                                            *
    ! Last update FT 28.03.2023                  *
    !                                            *
    !*********************************************

    IMPLICIT NONE

    CHARACTER(4):: eos_str
    CHARACTER(LEN=*), INTENT(IN):: eos_long

    IF( INDEX(eos_long, 'SLy ') /= 0 .OR. INDEX(eos_long, 'sly ') /= 0 )THEN

      eos_str= 'SLy '

    ELSEIF( INDEX(eos_long, 'ALF2') /= 0 .OR. INDEX(eos_long, 'alf2') /= 0 )THEN

      eos_str= 'ALF2'

    ELSEIF( INDEX(eos_long, 'ALF4') /= 0 .OR. INDEX(eos_long, 'alf4') /= 0 )THEN

      eos_str= 'ALF4'

    ELSEIF( INDEX(eos_long, 'ENG ') /= 0 .OR. INDEX(eos_long, 'eng ') /= 0 )THEN

      eos_str= 'ENG '

    ELSEIF( INDEX(eos_long, 'H4  ') /= 0 .OR. INDEX(eos_long, 'h4  ') /= 0 )THEN

      eos_str= 'H4  '

    ELSEIF( INDEX(eos_long, 'MPA1') /= 0 .OR. INDEX(eos_long, 'mpa1') /= 0 )THEN

      eos_str= 'MPA1'

    ELSEIF( INDEX(eos_long, 'MS1 ') /= 0 .OR. INDEX(eos_long, 'ms1 ') /= 0 )THEN

      eos_str= 'MS1 '

    ELSEIF( INDEX(eos_long, 'MS1b') /= 0 .OR. INDEX(eos_long, 'ms1b') /= 0 )THEN

      eos_str= 'MS1b'

    ELSEIF( INDEX(eos_long, 'APR3') /= 0 .OR. INDEX(eos_long, 'apr3') /= 0 )THEN

      eos_str= 'AP3'

    ELSEIF( INDEX(eos_long, 'AP3 ') /= 0 .OR. INDEX(eos_long, 'ap3 ') /= 0 )THEN

      eos_str= 'AP3 '

    ELSEIF( INDEX(eos_long, 'APR4') /= 0 .OR. INDEX(eos_long, 'apr4') /= 0 )THEN

      eos_str= 'APR4'

    ELSEIF( INDEX(eos_long, 'AP4 ') /= 0 .OR. INDEX(eos_long, 'ap4 ') /= 0 )THEN

      eos_str= 'APR4'

    ELSEIF( INDEX(eos_long, 'WFF1') /= 0 .OR. INDEX(eos_long, 'wff1') /= 0 )THEN

      eos_str= 'WFF1'

    ELSEIF( INDEX(eos_long, 'WFF2') /= 0 .OR. INDEX(eos_long, 'wff2') /= 0 )THEN

      eos_str= 'WFF2'

    ELSEIF( INDEX(eos_long, 'HASO') /= 0 .OR. INDEX(eos_long, 'haso') /= 0 )THEN

      eos_str= 'haso'

    ELSEIF( eos_long == 'Polytropic EOS' .OR. INDEX(eos_long, 'gam') /= 0 )THEN

      eos_str= 'POLY'

    ELSE

      PRINT *, "** ERROR! Unknown EOS name: ", TRIM(eos_long)
      PRINT *, " * Please add the name to SUBROUTINE shorten_eos_name. "
      PRINT *, " * Stopping..."
      PRINT *
      STOP

    ENDIF

  END FUNCTION shorten_eos_name


END MODULE utility