module arrays ! This module defines global arrays of ! fixed size Nz integer, parameter :: Nz = 10000 real(8), dimension (1:Nz) :: A real(8), dimension (1:Nz) :: B real(8), dimension (1:Nz) :: alpha real(8), dimension (1:Nz) :: beta real(8), dimension (1:Nz) :: r real(8), dimension (1:Nz) :: KB end module program horizon use arrays implicit none ! This program finds the apparent horizon of the Schwarschild spacetime ! using a few sets of coordinates. ! The element is writen as ! ds^2 = -( alpha^2 - A^2*betadn^2 )dt +2*A*betadn*dr*dt +A^2*dr^2 ! + B^2*r^2*( dtheta^2 + sin^2(theta)*dphi^2 ). ! The expansion of the null geodesics is then calculated as ! ! / 1 beta \ d ( Br ) KB !theta = | --- - ------ | ---- - ---- ! \ A alpha / dr Br real(8) rmax ! Maximum radius real(8) M ! Mass of the black hole real(8) dr ! grid size real(8) one, two ! numbers integer coordinates ! Set the type of coordinate for Schwarzschild spacetime integer i,j,k ! Counters integer Nr ! number of points ! ==================== ! Numbers ! ==================== one = 1.0d0 two = 2.0d0 ! ==================== ! Set some parameters ! ==================== ! Maximum number of points Nr Nr =10000 ! Mass of the black hole M = 1.0d0 ! rmax rmax = 50.0d0 ! Metric coordinates ! 1 => Boyer Lidquist ! 2 => isotropic print*, "type of coordinates" read*, coordinates ! size of step dr = rmax/dble(Nr) ! Grid do i = 1, Nr r(i) = (0.1+ dble(i) )*dr end do ! Schwarzschild if(coordinates==1) then A = one/(one-two*M/r) B = one beta = 0 alpha = one-two*M/r ! Isotropic else if (coordinates==2) then A = (one + 0.5d0*M/r)**2 B = A beta = 0.0 alpha = (1.0d0-0.5*M/r)/(1.0d0+0.5*M/r) ! Kerr-Schild else if (coordinates==3) then A = sqrt(one + two*M/r) B = one beta = two*M/r * one/(one+two*M/r) alpha = one /sqrt(one+two*M/r) else print*, "Coordinates not available" end if ! Print Coordinates print*, "==================================" print*, "= Apparent horizon parameters =" if (coordinates==1)then print*, "= Schwarzschild coordinates: =" print*, "==================================" else if (coordinates==2) then print*, "= Isotropic coordinates: =" print*, "==================================" else if(coordinates==3) then print*, "= Kerr-Schild coordinates: =" print*, "==================================" end if ! Find the apparent horizon call horizon_finder(Nr,dr) ! End end program !========================================= !========================================= subroutine horizon_finder(Nr,dr) use arrays implicit none integer i,Nr integer flag real(8) dr real(8) theta real(8) theta_old real(8) r_ah,ah_schw,ah_area,ah_mass,ah_A,ah_B,ah_alpha real(8) zero,half,one,two,four,smallpi,idr real(8), dimension (1:Nr) :: DB,betadn ! Numbers zero = 0.0d0 one = 1.0d0 two = 2.0d0 four = 4.0d0 half = 0.5d0 smallpi = acos(-one) ! Shift contra covariant betadn = beta*A**2 ! Calculate the derivative of the metric coefficient B do i =2,Nr-1 DB(i) = 0.5d0*(B(i+1)-B(i-1))/dr end do ! Horizon position. r_ah = zero ! Metric variables at horizon. ah_alpha = one ah_A = one ah_B = one ! Area, mass and Schwarzschild radius of horizons. ah_area = zero ah_mass = zero ah_schw = zero ! **************************** ! *** HORIZON EQUATION *** ! **************************** ! Initialize flag for horizon. flag = 0 ! Initialize trial horizon position. i = Nr-1 ! Start moving inward from boundary. do while ((flag.eq.0).and.(i.ge.1)) ! Find expansion. theta = - KB(i)/(r(i)*B(i)) theta = theta+ (one/A(i) - beta(i)/alpha(i))*(B(i)+ r(i)* DB(i)) ! If the expansion changes sign then find the horizon. if (theta.lt.zero) then ! Set horizon flag. flag = 1 ! Calculate the radius of the horizon. r_ah = (r(i+1)-r(i))/(theta_old-theta)*(-theta_old) + r(i+1) ! Calculate the metric function B at the horizon ! (required to calculate the radius in Schwarzschild coordinates). ah_B = (B(i+1)-B(i))/(r(i+1)-r(i))*(r_ah-r(i+1)) + B(i+1) ! Calculate the metric function A at the horizon. ah_A = (A(i+1)-A(i))/(r(i+1)-r(i))*(r_ah-r(i+1)) + A(i+1) ! Calculate the lapse at the horizon. ah_alpha = (alpha(i+1)-alpha(i))/(r(i+1)-r(i))*(r_ah-r(i+1)) + alpha(i+1) ! Calculate the radius of the horizon in Schwarzschild coordinates. ah_schw = r_ah*ah_B ! Calculate the area of the horizon. ah_area = 4.0d0*smallpi*(r_ah*ah_B)**2 ! Calculate the mass of the horizon. ah_mass = sqrt(ah_area/16.0D0/smallpi) endif ! Save value of theta theta_old = theta ! Move one grid point it. i = i - 1 end do ! ***************************** ! *** PRINT HORIZON DATA *** ! ***************************** ! Print parameters of the apparent horizon. print*, "| AH position:",real(r_ah) print*, "| AH area: ",real(ah_area) print*, "| AH mass: ",real(ah_mass) print*, "==================================" end subroutine