00001 subroutine bh_boundary_psi_test2_mpt(impt,char_bc,sou_surf,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   use grid_parameter_mpt,      only : grid_param_real_
00008   implicit none
00009   character(len=1), intent(in) :: char_bc
00010   real(long), pointer :: sou_surf(:,:), dsou_surf(:,:)
00011   real(long) :: st, cp, rad1, rad2, dr2dr1, bhmass
00012   integer    :: itg, ipg
00013   integer    :: impt, impt_ex
00014   real(long) :: bhmass_ex, rgin_ex
00015 
00016   if (impt.eq.1) impt_ex = 2
00017   if (impt.eq.2) impt_ex = 1
00018   rgin_ex = grid_param_real_(1,impt_ex)
00019   bhmass_ex = 2.0d0*rgin_ex*0.8d0
00020   bhmass = 2.0d0*rgin*0.8d0
00021 
00022   do ipg = 1, npg
00023     do itg = 1, ntg
00024       st = hsinthg(itg)
00025       cp = hcosphig(ipg)
00026       rad1 = rgin
00027       rad2 = sqrt(rad1**2 - 2.0d0*rad1*sepa*st*cp + sepa**2)
00028       dr2dr1 = (rad1 - sepa*st*cp)/rad2
00029       if (char_bc.eq.'d') then
00030         sou_surf(itg,ipg) = 1.0d0 + 0.5d0*bhmass/rad1 + 0.5d0*bhmass_ex/rad2
00031       end if
00032       if (char_bc.eq.'n') then
00033         dsou_surf(itg,ipg) = - 0.5d0*bhmass/rad1**2 &
00034         &                    - 0.5d0*bhmass_ex/rad2**2*dr2dr1
00035       end if
00036     end do
00037   end do
00038 
00039 end subroutine bh_boundary_psi_test2_mpt