place_surfaces Subroutine

subroutine place_surfaces(central_dens, center, radius, m_p, n_surfaces, surface_radii, last_r, get_dens)

Uses

  • proc~~place_surfaces~~UsesGraph proc~place_surfaces place_surfaces constants constants proc~place_surfaces->constants module~utility utility proc~place_surfaces->module~utility module~utility->constants matrix matrix module~utility->matrix

Place the surfaces, according to the baryon mass density of the star along the larger equatorial radius

FT 23.07.2021

Arguments

Type IntentOptional Attributes Name
double precision, intent(in) :: central_dens
double precision, intent(in) :: center
double precision, intent(in) :: radius
double precision, intent(in) :: m_p
integer, intent(in) :: n_surfaces
double precision, intent(inout), DIMENSION( n_surfaces ) :: surface_radii
double precision, intent(in) :: last_r
function get_dens(x, y, z) result(density)

Returns the baryon mass density at the desired point

Arguments
Type IntentOptional Attributes Name
double precision, intent(in) :: x

coordinate of the desired point
double precision, intent(in) :: y

coordinate of the desired point
double precision, intent(in) :: z

coordinate of the desired point
Return Value double precision

Called by

proc~~place_surfaces~~CalledByGraph proc~place_surfaces place_surfaces proc~place_particles_ellipsoidal_surfaces place_particles_ellipsoidal_surfaces proc~place_particles_ellipsoidal_surfaces->proc~place_surfaces interface~place_particles_ellipsoidal_surfaces place_particles_ellipsoidal_surfaces interface~place_particles_ellipsoidal_surfaces->proc~place_particles_ellipsoidal_surfaces

Contents

Variables

Source Code

place_surfaces

Variables

Type Visibility Attributes Name Initial
integer, private :: i
double precision, private :: lat
double precision, private :: long
integer, private :: phi
double precision, private :: rad
double precision, private :: rho_tmp
integer, private :: th

Source Code

  SUBROUTINE place_surfaces( central_dens, center, radius, m_p, n_surfaces, &
                             surface_radii, last_r, get_dens )

    !************************************************
    !
    !# Place the surfaces, according to
    !  the baryon mass density of the star
    !  along the larger equatorial radius
    !
    !  FT 23.07.2021
    !
    !************************************************

    USE constants,  ONLY: third, pi
    USE utility,    ONLY: zero, two, ten

    IMPLICIT NONE

    INTEGER,          INTENT( IN ):: n_surfaces
    DOUBLE PRECISION, INTENT( IN ):: central_dens, center, radius, m_p, last_r
    DOUBLE PRECISION, DIMENSION( n_surfaces ), INTENT( IN OUT ):: surface_radii
    INTERFACE
      FUNCTION get_dens( x, y, z ) RESULT( density )
        !! Returns the baryon mass density at the desired point
        DOUBLE PRECISION, INTENT(IN):: x
        !! \(x\) coordinate of the desired point
        DOUBLE PRECISION, INTENT(IN):: y
        !! \(y\) coordinate of the desired point
        DOUBLE PRECISION, INTENT(IN):: z
        !! \(z\) coordinate of the desired point
        DOUBLE PRECISION:: density
        !> Baryon mass density at \((x,y,z)\)
      END FUNCTION get_dens
    END INTERFACE

    INTEGER:: i, th, phi
    DOUBLE PRECISION:: rho_tmp, long, lat, rad

    surface_radii= zero

    IF( central_dens > zero )THEN
      surface_radii(1)= ( central_dens/m_p )**(-third)
    ELSE
      DO i= 1, 1000, 1
        IF( get_dens( center + radius*i/(ten*ten*ten), zero, zero ) > zero )THEN
          surface_radii(1)= &
      ( get_dens( center + radius*i/(ten*ten*ten), zero, zero )/m_p )**(-third)
          EXIT
        ENDIF
      ENDDO
    ENDIF

    DO itr= 2, n_surfaces, 1

      rho_tmp= get_dens( center + surface_radii( itr - 1 ), zero, zero )

      IF( rho_tmp <= 1.0D-13 )THEN

        rad_loop: DO i= 1, 1000, 1
          DO th= 0, 10, 1
            DO phi= 0, 20, 1

              long= phi/(two*ten)*(two*pi)
              lat = th/ten*pi
              rad = (center + surface_radii( itr - 1 ) + radius*i/(ten*ten*ten))

              IF( get_dens( rad*SIN(lat)*COS(long),&
                            rad*SIN(lat)*SIN(long), &
                            rad*COS(lat) ) > 1.0D-13 )THEN

                rho_tmp= get_dens( rad*SIN(lat)*COS(long),&
                                   rad*SIN(lat)*SIN(long), &
                                   rad*COS(lat) )

                EXIT rad_loop

              ENDIF

            ENDDO
          ENDDO
        ENDDO rad_loop

      ENDIF

      IF( rho_tmp == zero )THEN
        surface_radii= surface_radii*itr/n_surfaces
      ENDIF

      surface_radii( itr )= surface_radii( itr - 1 ) + ( rho_tmp/m_p )**(-third)

    ENDDO
    surface_radii= surface_radii*(radius*last_r/surface_radii(n_surfaces))

  END SUBROUTINE place_surfaces