00001 subroutine test_analytic_solution_bhex_mpt(impt)
00002   use phys_constant, only  :   long, pi
00003   use grid_parameter, only  :   nrg, ntg, npg, nrf, rgin
00004   use def_metric, only  :   bvxd, bvyd
00005   use coordinate_grav_r, only : rg
00006   use grid_points_binary_excision, only : rb
00007   use grid_parameter_mpt, only : grid_param_real_
00008   use grid_points_asymptotic_patch_mpt, only : ra_
00009   use make_array_3d
00010   use copy_array_4dto3d_mpt
00011   implicit none
00012   integer     ::   irg,itg,ipg, impt, impt_ce, impt_ex
00013   real(long)  ::   zfac, small = 1.0d-15, alps_tmp
00014   real(long)  ::   mass_ce, mass_ex, rgin_ce, rgin_ex, rad_1, rad_2
00015   real(long), pointer :: ra_1(:,:,:), ra_2(:,:,:)
00016 
00017   call alloc_array3d(ra_1,-2,nrg+2,0,ntg,0,npg)
00018   call alloc_array3d(ra_2,-2,nrg+2,0,ntg,0,npg)
00019 
00020   if (impt.eq.1) then ; impt_ce = impt ; impt_ex = 2 ; end if
00021   if (impt.eq.2) then ; impt_ce = impt ; impt_ex = 1 ; end if
00022   if (impt.eq.3) then ; impt_ce = 1    ; impt_ex = 2 ; end if
00023 
00024   rgin_ce = grid_param_real_(1,impt_ce)
00025   rgin_ex = grid_param_real_(1,impt_ex)
00026   mass_ce = 2.0d0*rgin_ce
00027   mass_ex = 2.0d0*rgin_ex
00028 
00029   if (impt.eq.3) then 
00030     call copy_array4dto3d_mpt(impt_ce,ra_,ra_1,-2,nrg+2,0,ntg,0,npg)
00031     call copy_array4dto3d_mpt(impt_ex,ra_,ra_2,-2,nrg+2,0,ntg,0,npg)
00032   end if
00033 
00034   do ipg = 0, npg
00035     do itg = 0, ntg
00036       do irg = 0, nrg
00037         if (impt.ne.3) then 
00038           rad_1 = rg(irg)
00039           rad_2 = rb(irg,itg,ipg)
00040         end if 
00041         if (impt.eq.3) then 
00042           rad_1 = ra_1(irg,itg,ipg) 
00043           rad_2 = ra_2(irg,itg,ipg)
00044         end if
00045         bvxd(irg,itg,ipg) = 1.0d0 + mass_ce/(2.0d0*rad_1) &
00046         &                         + mass_ex/(2.0d0*rad_2)
00047         alps_tmp          = 1.0d0 - mass_ce/(2.0d0*rad_1) &
00048         &                         - mass_ex/(2.0d0*rad_2)
00049         bvyd(irg,itg,ipg) = alps_tmp/bvxd(irg,itg,ipg)
00050       end do
00051     end do
00052   end do
00053 
00054   deallocate(ra_1)
00055   deallocate(ra_2)
00056 end subroutine test_analytic_solution_bhex_mpt