00001 subroutine bh_boundary_nh_alph_test(dsou_surf)
00002 use phys_constant, only : long
00003 use grid_parameter, only : nrg, ntg, npg, rgin
00004 use trigonometry_grav_theta, only : hsinthg
00005 use trigonometry_grav_phi, only : hcosphig
00006 use def_binary_parameter, only : sepa
00007 implicit none
00008
00009
00010 real(long), pointer :: dsou_surf(:,:)
00011 real(long) :: st, cp, rad1, rad2, dr2dr1, bhmass
00012 real(long) :: psisurf, alphsurf, dpsidr1surf
00013 integer :: itg, ipg
00014
00015 bhmass = 2.0d0*rgin
00016 do ipg = 1, npg
00017 do itg = 1, ntg
00018 st = hsinthg(itg)
00019 cp = hcosphig(ipg)
00020 rad1 = rgin
00021 rad2 = sqrt(rad1**2 - 2.0d0*rad1*sepa*st*cp + sepa**2)
00022 dr2dr1 = (rad1 - sepa*st*cp)/rad2
00023
00024 psisurf = 1.0d0 + 0.5d0*bhmass/rad1 + 0.5d0*bhmass/rad2
00025 dpsidr1surf = -0.5d0*bhmass/rad1**2 - 0.5d0*bhmass/rad2**2*dr2dr1
00026 alphsurf = (1.0d0 - 0.5d0*bhmass/rad1 - 0.5d0*bhmass/rad2)/psisurf
00027
00028 dsou_surf(itg,ipg) = (0.5d0*bhmass/rad1**2 &
00029 + 0.5d0*bhmass/rad2**2*dr2dr1 &
00030 & - alphsurf*dpsidr1surf)/psisurf
00031 end do
00032 end do
00033
00034 end subroutine bh_boundary_nh_alph_test