00001
00002 include '../Include_file/include_modulefiles_plot_BNS_CF_3mpt.f90'
00003 include '../Include_file/include_PEOS_modulefile.f90'
00004 include '../Include_file/include_modulefiles_analysis_plot_BNS_CF_3mpt.f90'
00005 include '../Include_file/include_interface_modulefiles_plot_BNS_CF_3mpt.f90'
00006 include '../Include_file/include_interface_modulefiles_analysis_plot_BNS_CF_3mpt.f90'
00007 include '../Include_file/include_subroutines_plot_BNS_CF_3mpt.f90'
00008 include '../Include_file/include_subroutines_analysis_plot_BNS_CF_3mpt.f90'
00009 include '../Include_file/include_PEOS_subroutines.f90'
00010 include '../Include_file/include_functions.f90'
00011
00012
00013
00014
00015 PROGRAM coc2cac
00016
00017 use grid_parameter_binary_excision
00018 use phys_constant
00019 use grid_parameter
00020 use interface_modules_cartesian
00021 use interface_calc_gradvep_export
00022 use trigonometry_grav_phi
00023 use def_binary_parameter, only : dis
00024 use def_matter_parameter_mpt
00025 use interface_IO_input_CF_grav_export
00026 use interface_IO_input_CF_flsp_export
00027 use interface_IO_input_CF_flco_export
00028 use interface_IO_input_CF_flir_export
00029 use interface_IO_input_CF_surf_export
00030 use interface_excurve_CF_gridpoint_export
00031 use interface_interpo_gr2fl_metric_CF_export
00032 implicit none
00033 integer :: impt, impt_ex, ico, irr, isp, ipo, npoints, iofile
00034 real(8) :: confpow, psifcacp, wxspfca, wyspfca, wzspfca
00035 character(30) :: char1
00036 character*400 :: dir_path
00037 character(len=1) :: np(5) = (/'1', '2','3', '4', '5'/)
00038 real(8) :: rr3, rrcm, xc,yc,zc, xc_p1, yc_p1, zc_p1, xc_p2, yc_p2, zc_p2, dis_cm
00039 real(8) :: xc_p3, yc_p3, zc_p3
00040 real(8) :: xcac, ycac, zcac
00041 real(8) :: xcoc, ycoc, zcoc
00042 real(8) :: emdca, vepca, psica, alphca, bvxdca, bvydca, bvzdca, psi4ca, psif4ca
00043 real(8) :: hca, preca, rhoca, eneca, epsca
00044 real(8) :: vepxfca, vepyfca, vepzfca, vxu, vyu, vzu, lam_p1, lam_p2
00045 real(8) :: bxcor, bycor, bzcor, bvxdfca, bvydfca, bvzdfca, psifca, alphfca
00046 real(8) :: gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, kzz
00047 real(8) :: axx, axy, axz, ayy, ayz, azz
00048 real(8) :: ome_p1, ber_p1, radi_p1, r_surf_p1
00049 real(8) :: ome_p2, ber_p2, radi_p2, r_surf_p2
00050
00051 real(8), pointer :: emd_p1(:,:,:), vep_p1(:,:,:), rs_p1(:,:)
00052 real(8), pointer :: vepxf_p1(:,:,:), vepyf_p1(:,:,:), vepzf_p1(:,:,:)
00053 real(8), pointer :: wxspf_p1(:,:,:), wyspf_p1(:,:,:), wzspf_p1(:,:,:)
00054 real(8), pointer :: psif_p1(:,:,:), alphf_p1(:,:,:), bvxdf_p1(:,:,:), bvydf_p1(:,:,:), bvzdf_p1(:,:,:)
00055 real(8), pointer :: psi_p1(:,:,:), alph_p1(:,:,:), bvxd_p1(:,:,:), bvyd_p1(:,:,:), bvzd_p1(:,:,:)
00056 real(8), pointer :: axx_p1(:,:,:), axy_p1(:,:,:) , axz_p1(:,:,:) , ayy_p1(:,:,:) , ayz_p1(:,:,:), azz_p1(:,:,:)
00057
00058 real(8), pointer :: emd_p2(:,:,:), vep_p2(:,:,:), rs_p2(:,:)
00059 real(8), pointer :: vepxf_p2(:,:,:), vepyf_p2(:,:,:), vepzf_p2(:,:,:)
00060 real(8), pointer :: wxspf_p2(:,:,:), wyspf_p2(:,:,:), wzspf_p2(:,:,:)
00061 real(8), pointer :: psif_p2(:,:,:), alphf_p2(:,:,:), bvxdf_p2(:,:,:), bvydf_p2(:,:,:), bvzdf_p2(:,:,:)
00062 real(8), pointer :: psi_p2(:,:,:), alph_p2(:,:,:), bvxd_p2(:,:,:), bvyd_p2(:,:,:), bvzd_p2(:,:,:)
00063 real(8), pointer :: axx_p2(:,:,:), axy_p2(:,:,:) , axz_p2(:,:,:) , ayy_p2(:,:,:) , ayz_p2(:,:,:), azz_p2(:,:,:)
00064
00065 real(8), pointer :: psi_p3(:,:,:), alph_p3(:,:,:), bvxd_p3(:,:,:), bvyd_p3(:,:,:), bvzd_p3(:,:,:)
00066 real(8), pointer :: axx_p3(:,:,:), axy_p3(:,:,:) , axz_p3(:,:,:) , ayy_p3(:,:,:) , ayz_p3(:,:,:), azz_p3(:,:,:)
00067
00068 integer :: kij_parity
00069 confpow = -6.0d0
00070
00071 gxx=0.0d0; gxy=0.0d0; gxz=0.0d0; gyy=0.0d0; gyz=0.0d0; gzz=0.0d0
00072 kxx=0.0d0; kxy=0.0d0; kxz=0.0d0; kyy=0.0d0; kyz=0.0d0; kzz=0.0d0
00073 axx=0.0d0; axy=0.0d0; axz=0.0d0; ayy=0.0d0; ayz=0.0d0; azz=0.0d0
00074
00075
00076
00077
00078 dir_path='.'
00079
00080
00081 call allocate_grid_parameter_mpt
00082 call allocate_grid_parameter_binary_excision_mpt
00083 call allocate_def_matter_parameter_mpt
00084 do impt = 1, nmpt
00085 call read_parameter_mpt_cactus(impt,dir_path)
00086 if (impt .le. 2) call read_surf_parameter_mpt_cactus(impt,dir_path)
00087 call copy_grid_parameter_to_mpt(impt)
00088 call read_parameter_binary_excision_mpt_cactus(impt,dir_path)
00089 call copy_grid_parameter_binary_excision_to_mpt(impt)
00090 if (impt .le. 2) call peos_initialize_mpt_cactus(impt,dir_path)
00091 call copy_def_peos_parameter_to_mpt(impt)
00092 end do
00093
00094 call set_allocate_size_mpt
00095 call allocate_trig_grav_mphi
00096 call allocate_trigonometry_grav_phi_mpt
00097
00098 do impt = 1, nmpt
00099 call copy_grid_parameter_from_mpt(impt)
00100 call copy_grid_parameter_binary_excision_from_mpt(impt)
00101 call copy_def_peos_parameter_from_mpt(impt)
00102 call coordinate_patch_kit_grav_grid_coc2cac_mpt(3)
00103 call calc_parameter_binary_excision
00104 call copy_grid_parameter_to_mpt(impt)
00105 call copy_grid_parameter_binary_excision_to_mpt(impt)
00106 call copy_coordinate_grav_extended_to_mpt(impt)
00107 call copy_coordinate_grav_phi_to_mpt(impt)
00108 call copy_coordinate_grav_r_to_mpt(impt)
00109 call copy_coordinate_grav_theta_to_mpt(impt)
00110 call copy_def_binary_parameter_to_mpt(impt)
00111 call copy_trigonometry_grav_phi_to_mpt(impt)
00112 call copy_trigonometry_grav_theta_to_mpt(impt)
00113 end do
00114
00115 do impt = 1,nmpt
00116
00117 call copy_grid_parameter_from_mpt(impt)
00118 call copy_grid_parameter_binary_excision_from_mpt(impt)
00119 call copy_coordinate_grav_extended_from_mpt(impt)
00120 call copy_coordinate_grav_phi_from_mpt(impt)
00121 call copy_coordinate_grav_r_from_mpt(impt)
00122 call copy_coordinate_grav_theta_from_mpt(impt)
00123 call copy_def_binary_parameter_from_mpt(impt)
00124 call copy_trigonometry_grav_phi_from_mpt(impt)
00125 call copy_trigonometry_grav_theta_from_mpt(impt)
00126
00127
00128 if (impt==1) then
00129 if(NS_shape.eq.'SP') then
00130 ico=0; irr=0; isp=1
00131 write(6,*) "************ Spinning configuration ************"
00132 end if
00133 if(NS_shape.eq.'IR') then
00134 ico=0; irr=1; isp=0
00135 write(6,*) "************ Irrotational configuration ************"
00136 end if
00137 if(NS_shape.eq.'CO') then
00138 ico=1; irr=0; isp=0
00139 write(6,*) "************ Corotating configuration ************"
00140 end if
00141
00142 rr3 = 0.7d0*(rgout - rgmid)
00143 dis_cm = dis
00144 r_surf_p1 = r_surf
00145 write(6,'(a13,i1,a8)') "Reading COCP-", impt, " data..."
00146
00147 allocate ( emd_p1(0:nrf,0:ntf,0:npf))
00148 allocate ( vep_p1(0:nrf,0:ntf,0:npf))
00149 allocate (wxspf_p1(0:nrf,0:ntf,0:npf))
00150 allocate (wyspf_p1(0:nrf,0:ntf,0:npf))
00151 allocate (wzspf_p1(0:nrf,0:ntf,0:npf))
00152 allocate (vepxf_p1(0:nrf,0:ntf,0:npf))
00153 allocate (vepyf_p1(0:nrf,0:ntf,0:npf))
00154 allocate (vepzf_p1(0:nrf,0:ntf,0:npf))
00155 allocate ( psif_p1(0:nrf,0:ntf,0:npf))
00156 allocate (alphf_p1(0:nrf,0:ntf,0:npf))
00157 allocate (bvxdf_p1(0:nrf,0:ntf,0:npf))
00158 allocate (bvydf_p1(0:nrf,0:ntf,0:npf))
00159 allocate (bvzdf_p1(0:nrf,0:ntf,0:npf))
00160 allocate ( rs_p1(0:ntf,0:npf))
00161 allocate ( psi_p1(0:nrg,0:ntg,0:npg))
00162 allocate ( alph_p1(0:nrg,0:ntg,0:npg))
00163 allocate ( bvxd_p1(0:nrg,0:ntg,0:npg))
00164 allocate ( bvyd_p1(0:nrg,0:ntg,0:npg))
00165 allocate ( bvzd_p1(0:nrg,0:ntg,0:npg))
00166 allocate ( axx_p1(0:nrg,0:ntg,0:npg))
00167 allocate ( axy_p1(0:nrg,0:ntg,0:npg))
00168 allocate ( axz_p1(0:nrg,0:ntg,0:npg))
00169 allocate ( ayy_p1(0:nrg,0:ntg,0:npg))
00170 allocate ( ayz_p1(0:nrg,0:ntg,0:npg))
00171 allocate ( azz_p1(0:nrg,0:ntg,0:npg))
00172 emd_p1=0.0d0; vep_p1 =0.0d0; rs_p1 =0.0d0; wxspf_p1=0.0d0; wyspf_p1=0.0d0; wzspf_p1=0.0d0
00173 psi_p1=0.0d0; alph_p1=0.0d0; bvxd_p1=0.0d0; bvyd_p1=0.0d0; bvzd_p1=0.0d0
00174 axx_p1=0.0d0; axy_p1 =0.0d0; axz_p1 =0.0d0; ayy_p1=0.0d0; ayz_p1=0.0d0; azz_p1=0.0d0
00175
00176
00177
00178
00179
00180 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt1.las",psi_p1,alph_p1,bvxd_p1,bvyd_p1,bvzd_p1)
00181
00182 if(NS_shape.eq.'CO') call IO_input_CF_flco_export(trim(dir_path)//"/bnsflu_3D_mpt1.las",emd_p1,ome_p1,ber_p1,radi_p1)
00183 if(NS_shape.eq.'IR') call IO_input_CF_flir_export(trim(dir_path)//"/bnsflu_3D_mpt1.las",emd_p1,vep_p1,ome_p1,ber_p1,radi_p1)
00184 if(NS_shape.eq.'SP') call IO_input_CF_flsp_export(trim(dir_path)//"/bnsflu_3D_mpt1.las",emd_p1,vep_p1,wxspf_p1,wyspf_p1, &
00185 & wzspf_p1,ome_p1,ber_p1,radi_p1)
00186
00187 call IO_input_CF_surf_export(trim(dir_path)//"/bnssur_3D_mpt1.las",rs_p1)
00188
00189 call excurve_CF_gridpoint_export(alph_p1,bvxd_p1,bvyd_p1,bvzd_p1, &
00190 & axx_p1, axy_p1, axz_p1, ayy_p1, ayz_p1, azz_p1)
00191
00192 call interpo_gr2fl_metric_CF_export(alph_p1, psi_p1, bvxd_p1, bvyd_p1, bvzd_p1, &
00193 & alphf_p1, psif_p1, bvxdf_p1, bvydf_p1, bvzdf_p1, rs_p1)
00194
00195 if((NS_shape.eq.'IR').or.(NS_shape.eq.'SP')) call calc_gradvep_export(vep_p1,vepxf_p1,vepyf_p1,vepzf_p1,rs_p1)
00196 end if
00197 if (impt==2) then
00198 if(NS_shape.eq.'SP') then
00199 write(6,*) "************ Spinning configuration ************"
00200 if (isp.ne.1) stop
00201 end if
00202 if(NS_shape.eq.'IR') then
00203 write(6,*) "************ Irrotational configuration ************"
00204 if(irr.ne.1) stop
00205 endif
00206 if(NS_shape.eq.'CO') then
00207 write(6,*) "************ Corotating configuration ************"
00208 if(ico.ne.1) stop
00209 end if
00210 r_surf_p2 = r_surf
00211 write(6,'(a13,i1,a8)') "Reading COCP-", impt, " data..."
00212 allocate ( emd_p2(0:nrf,0:ntf,0:npf))
00213 allocate ( vep_p2(0:nrf,0:ntf,0:npf))
00214 allocate (wxspf_p2(0:nrf,0:ntf,0:npf))
00215 allocate (wyspf_p2(0:nrf,0:ntf,0:npf))
00216 allocate (wzspf_p2(0:nrf,0:ntf,0:npf))
00217 allocate (vepxf_p2(0:nrf,0:ntf,0:npf))
00218 allocate (vepyf_p2(0:nrf,0:ntf,0:npf))
00219 allocate (vepzf_p2(0:nrf,0:ntf,0:npf))
00220 allocate ( psif_p2(0:nrf,0:ntf,0:npf))
00221 allocate (alphf_p2(0:nrf,0:ntf,0:npf))
00222 allocate (bvxdf_p2(0:nrf,0:ntf,0:npf))
00223 allocate (bvydf_p2(0:nrf,0:ntf,0:npf))
00224 allocate (bvzdf_p2(0:nrf,0:ntf,0:npf))
00225 allocate ( rs_p2(0:ntf,0:npf))
00226 allocate ( psi_p2(0:nrg,0:ntg,0:npg))
00227 allocate ( alph_p2(0:nrg,0:ntg,0:npg))
00228 allocate ( bvxd_p2(0:nrg,0:ntg,0:npg))
00229 allocate ( bvyd_p2(0:nrg,0:ntg,0:npg))
00230 allocate ( bvzd_p2(0:nrg,0:ntg,0:npg))
00231 allocate ( axx_p2(0:nrg,0:ntg,0:npg))
00232 allocate ( axy_p2(0:nrg,0:ntg,0:npg))
00233 allocate ( axz_p2(0:nrg,0:ntg,0:npg))
00234 allocate ( ayy_p2(0:nrg,0:ntg,0:npg))
00235 allocate ( ayz_p2(0:nrg,0:ntg,0:npg))
00236 allocate ( azz_p2(0:nrg,0:ntg,0:npg))
00237 emd_p2=0.0d0; vep_p2=0.0d0; rs_p2=0.0d0; wxspf_p2=0.0d0; wyspf_p2=0.0d0; wzspf_p2=0.0d0
00238 psi_p2=0.0d0; alph_p2=0.0d0; bvxd_p2=0.0d0; bvyd_p2=0.0d0; bvzd_p2=0.0d0
00239 axx_p2=0.0d0; axy_p2=0.0d0; axz_p2=0.0d0; ayy_p2=0.0d0; ayz_p2=0.0d0; azz_p2=0.0d0
00240
00241
00242
00243
00244
00245 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt2.las",psi_p2,alph_p2,bvxd_p2,bvyd_p2,bvzd_p2)
00246
00247 if(NS_shape.eq.'CO') call IO_input_CF_flco_export(trim(dir_path)//"/bnsflu_3D_mpt2.las",emd_p2,ome_p2,ber_p2,radi_p2)
00248 if(NS_shape.eq.'IR') call IO_input_CF_flir_export(trim(dir_path)//"/bnsflu_3D_mpt2.las",emd_p2,vep_p2,ome_p2,ber_p2,radi_p2)
00249 if(NS_shape.eq.'SP') call IO_input_CF_flsp_export(trim(dir_path)//"/bnsflu_3D_mpt2.las",emd_p2,vep_p2,wxspf_p2,wyspf_p2, &
00250 & wzspf_p2,ome_p2,ber_p2,radi_p2)
00251
00252 call IO_input_CF_surf_export(trim(dir_path)//"/bnssur_3D_mpt2.las",rs_p2)
00253
00254 call excurve_CF_gridpoint_export(alph_p2,bvxd_p2,bvyd_p2,bvzd_p2, &
00255 & axx_p2, axy_p2, axz_p2, ayy_p2, ayz_p2, azz_p2)
00256
00257 call interpo_gr2fl_metric_CF_export(alph_p2, psi_p2, bvxd_p2, bvyd_p2, bvzd_p2, &
00258 & alphf_p2, psif_p2, bvxdf_p2, bvydf_p2, bvzdf_p2, rs_p2)
00259
00260 if((NS_shape.eq.'IR').or.(NS_shape.eq.'SP')) call calc_gradvep_export(vep_p2,vepxf_p2,vepyf_p2,vepzf_p2,rs_p2)
00261
00262 vepxf_p2 = - vepxf_p2
00263 vepyf_p2 = - vepyf_p2
00264 bvxdf_p2 = - bvxdf_p2
00265 bvydf_p2 = - bvydf_p2
00266 bvxd_p2 = - bvxd_p2
00267 bvyd_p2 = - bvyd_p2
00268 end if
00269 if (impt==3) then
00270 write(6,'(a13,i1,a8)') "Reading ARCP-", impt, " data..."
00271 allocate ( psi_p3(0:nrg,0:ntg,0:npg))
00272 allocate (alph_p3(0:nrg,0:ntg,0:npg))
00273 allocate (bvxd_p3(0:nrg,0:ntg,0:npg))
00274 allocate (bvyd_p3(0:nrg,0:ntg,0:npg))
00275 allocate (bvzd_p3(0:nrg,0:ntg,0:npg))
00276 allocate ( axx_p3(0:nrg,0:ntg,0:npg))
00277 allocate ( axy_p3(0:nrg,0:ntg,0:npg))
00278 allocate ( axz_p3(0:nrg,0:ntg,0:npg))
00279 allocate ( ayy_p3(0:nrg,0:ntg,0:npg))
00280 allocate ( ayz_p3(0:nrg,0:ntg,0:npg))
00281 allocate ( azz_p3(0:nrg,0:ntg,0:npg))
00282 psi_p3=0.0d0; alph_p3=0.0d0; bvxd_p3=0.0d0; bvyd_p3=0.0d0; bvzd_p3=0.0d0
00283 axx_p3=0.0d0; axy_p3=0.0d0; axz_p3=0.0d0; ayy_p3=0.0d0; ayz_p3=0.0d0; azz_p3=0.0d0
00284
00285 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt3.las",psi_p3,alph_p3,bvxd_p3,bvyd_p3,bvzd_p3)
00286
00287 call excurve_CF_gridpoint_export(alph_p3,bvxd_p3,bvyd_p3,bvzd_p3, &
00288 & axx_p3, axy_p3, axz_p3, ayy_p3, ayz_p3, azz_p3)
00289 end if
00290 end do
00291 write(6,'(2e20.12)') emd_p1(0,0,0), emd_p1(58,0,0)
00292 write(6,'(3e20.12)') ome_p1, ber_p1, radi_p1
00293 write(6,'(e20.12)') dis_cm
00294
00295 open(850,file='rho_ini0.txt',status='old')
00296 npoints=0
00297 do
00298 read(850,*,iostat=iofile) xcac
00299 if (iofile > 0) then
00300 write(6,*) "Problem reading file rho_ini0.txt...exiting"
00301 stop
00302 else if (iofile < 0) THEN
00303 write(6,*) "Total number of points is =", npoints
00304 exit
00305 else
00306 npoints = npoints + 1
00307 end if
00308 end do
00309 close(850)
00310
00311 open(860,file='rho_ini0.txt',status='old')
00312 write(6,'(21a16)') "# xcac", "psica", "alphca", "bvxdca", "bvydca", "bvzdca", "preca", "rhoca", "epsca", &
00313 & "vxu", "vyu", "vzu", "vepxfca", "vepyfca", "vepzfca", "kxx", "kxy", "kxz", "kyy", "kyz", "kzz"
00314
00315 do ipo=1,npoints
00316 read(860,*,iostat=iofile) xcac
00317 ycac=0.0d0
00318 zcac=0.0d0
00319
00320
00321
00322
00323
00324
00325 xcoc = xcac/(radi_p1)
00326 ycoc = ycac/(radi_p1)
00327 zcoc = zcac/(radi_p1)
00328
00329
00330 rrcm = sqrt(xcoc**2 + ycoc**2 + zcoc**2)
00331
00332 if (rrcm > rr3) then
00333
00334 call copy_grid_parameter_from_mpt(3)
00335 call copy_grid_parameter_binary_excision_from_mpt(3)
00336 call copy_coordinate_grav_extended_from_mpt(3)
00337 call copy_coordinate_grav_phi_from_mpt(3)
00338 call copy_coordinate_grav_r_from_mpt(3)
00339 call copy_coordinate_grav_theta_from_mpt(3)
00340 call copy_def_binary_parameter_from_mpt(3)
00341 call copy_trigonometry_grav_phi_from_mpt(3)
00342 call copy_trigonometry_grav_theta_from_mpt(3)
00343
00344 kij_parity = 1
00345 xc_p3 = xcoc
00346 yc_p3 = ycoc
00347 zc_p3 = zcoc
00348
00349
00350 call interpo_gr2cgr_4th(psi_p3 , psica , xc_p3, yc_p3, zc_p3)
00351 call interpo_gr2cgr_4th(alph_p3, alphca, xc_p3, yc_p3, zc_p3)
00352 call interpo_gr2cgr_4th(bvxd_p3, bvxdca, xc_p3, yc_p3, zc_p3)
00353 call interpo_gr2cgr_4th(bvyd_p3, bvydca, xc_p3, yc_p3, zc_p3)
00354 call interpo_gr2cgr_4th(bvzd_p3, bvzdca, xc_p3, yc_p3, zc_p3)
00355 call interpo_gr2cgr_4th(axx_p3 , axx , xc_p3, yc_p3, zc_p3)
00356 call interpo_gr2cgr_4th(axy_p3 , axy , xc_p3, yc_p3, zc_p3)
00357 call interpo_gr2cgr_4th(axz_p3 , axz , xc_p3, yc_p3, zc_p3)
00358 call interpo_gr2cgr_4th(ayy_p3 , ayy , xc_p3, yc_p3, zc_p3)
00359 call interpo_gr2cgr_4th(ayz_p3 , ayz , xc_p3, yc_p3, zc_p3)
00360 call interpo_gr2cgr_4th(azz_p3 , azz , xc_p3, yc_p3, zc_p3)
00361 psi4ca = psica**4
00362 vxu = 0.0d0
00363 vyu = 0.0d0
00364 vzu = 0.0d0
00365 emdca = 0.0d0
00366 else
00367 if (xcoc<=0.0d0) then
00368
00369 call copy_grid_parameter_from_mpt(1)
00370 call copy_grid_parameter_binary_excision_from_mpt(1)
00371 call copy_coordinate_grav_extended_from_mpt(1)
00372 call copy_coordinate_grav_phi_from_mpt(1)
00373 call copy_coordinate_grav_r_from_mpt(1)
00374 call copy_coordinate_grav_theta_from_mpt(1)
00375 call copy_def_binary_parameter_from_mpt(1)
00376 call copy_trigonometry_grav_phi_from_mpt(1)
00377 call copy_trigonometry_grav_theta_from_mpt(1)
00378
00379 kij_parity = 1
00380 xc_p1 = xcoc + dis_cm
00381 yc_p1 = ycoc
00382 zc_p1 = zcoc
00383
00384 call interpo_gr2cgr_4th(psi_p1 , psica , xc_p1, yc_p1, zc_p1)
00385 call interpo_gr2cgr_4th(alph_p1, alphca, xc_p1, yc_p1, zc_p1)
00386 call interpo_gr2cgr_4th(bvxd_p1, bvxdca, xc_p1, yc_p1, zc_p1)
00387 call interpo_gr2cgr_4th(bvyd_p1, bvydca, xc_p1, yc_p1, zc_p1)
00388 call interpo_gr2cgr_4th(bvzd_p1, bvzdca, xc_p1, yc_p1, zc_p1)
00389 call interpo_gr2cgr_4th(axx_p1 , axx , xc_p1, yc_p1, zc_p1)
00390 call interpo_gr2cgr_4th(axy_p1 , axy , xc_p1, yc_p1, zc_p1)
00391 call interpo_gr2cgr_4th(axz_p1 , axz , xc_p1, yc_p1, zc_p1)
00392 call interpo_gr2cgr_4th(ayy_p1 , ayy , xc_p1, yc_p1, zc_p1)
00393 call interpo_gr2cgr_4th(ayz_p1 , ayz , xc_p1, yc_p1, zc_p1)
00394 call interpo_gr2cgr_4th(azz_p1 , azz , xc_p1, yc_p1, zc_p1)
00395 call interpo_fl2cgr_4th_export(emd_p1 , emdca , xc_p1, yc_p1, zc_p1, rs_p1)
00396
00397 call interpo_fl2cgr_4th_export(vepxf_p1, vepxfca , xc_p1, yc_p1, zc_p1, rs_p1)
00398 call interpo_fl2cgr_4th_export(vepyf_p1, vepyfca , xc_p1, yc_p1, zc_p1, rs_p1)
00399 call interpo_fl2cgr_4th_export(vepzf_p1, vepzfca , xc_p1, yc_p1, zc_p1, rs_p1)
00400
00401 call interpo_fl2cgr_4th_export(wxspf_p1, wxspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00402 call interpo_fl2cgr_4th_export(wyspf_p1, wyspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00403 call interpo_fl2cgr_4th_export(wzspf_p1, wzspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00404
00405 call interpo_fl2cgr_4th_export(psif_p1 , psifca , xc_p1, yc_p1, zc_p1, rs_p1)
00406 call interpo_fl2cgr_4th_export(alphf_p1, alphfca , xc_p1, yc_p1, zc_p1, rs_p1)
00407 call interpo_fl2cgr_4th_export(bvxdf_p1, bvxdfca , xc_p1, yc_p1, zc_p1, rs_p1)
00408 call interpo_fl2cgr_4th_export(bvydf_p1, bvydfca , xc_p1, yc_p1, zc_p1, rs_p1)
00409 call interpo_fl2cgr_4th_export(bvzdf_p1, bvzdfca , xc_p1, yc_p1, zc_p1, rs_p1)
00410 bxcor = bvxdfca + ome_p1*(-ycoc)
00411 bycor = bvydfca + ome_p1*(xcoc)
00412 bzcor = bvzdfca
00413 lam_p1 = ber_p1 + bxcor*vepxfca + bycor*vepyfca + bzcor*vepzfca
00414 psi4ca = psica**4
00415 psif4ca = psifca**4
00416 psifcacp= psifca**confpow
00417
00418 if (dabs(emdca) > 1.0d-14) then
00419 vxu = ico*( ome_p1*(-ycoc) ) + &
00420 & irr*( alphfca*vepxfca/psif4ca/lam_p1 ) + &
00421 & isp*( alphfca*(vepxfca/psif4ca + psifcacp*wxspfca)/lam_p1 )
00422
00423 vyu = ico*( ome_p1*( xcoc) ) + &
00424 & irr*( alphfca*vepyfca/psif4ca/lam_p1 ) + &
00425 & isp*( alphfca*(vepyfca/psif4ca + psifcacp*wyspfca)/lam_p1 )
00426
00427 vzu = ico*( 0.0d0 ) + &
00428 & irr*( alphfca*vepzfca/psif4ca/lam_p1 ) + &
00429 & isp*( alphfca*(vepzfca/psif4ca + psifcacp*wzspfca)/lam_p1 )
00430 else
00431 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00432 end if
00433 else
00434
00435 call copy_grid_parameter_from_mpt(2)
00436 call copy_grid_parameter_binary_excision_from_mpt(2)
00437 call copy_coordinate_grav_extended_from_mpt(2)
00438 call copy_coordinate_grav_phi_from_mpt(2)
00439 call copy_coordinate_grav_r_from_mpt(2)
00440 call copy_coordinate_grav_theta_from_mpt(2)
00441 call copy_def_binary_parameter_from_mpt(2)
00442 call copy_trigonometry_grav_phi_from_mpt(2)
00443 call copy_trigonometry_grav_theta_from_mpt(2)
00444
00445 kij_parity = -1
00446 xc_p2 = -(xcoc - dis_cm)
00447 yc_p2 = -ycoc
00448 zc_p2 = zcoc
00449
00450 call interpo_gr2cgr_4th(psi_p2 , psica , xc_p2, yc_p2, zc_p2)
00451 call interpo_gr2cgr_4th(alph_p2, alphca, xc_p2, yc_p2, zc_p2)
00452 call interpo_gr2cgr_4th(bvxd_p2, bvxdca, xc_p2, yc_p2, zc_p2)
00453 call interpo_gr2cgr_4th(bvyd_p2, bvydca, xc_p2, yc_p2, zc_p2)
00454 call interpo_gr2cgr_4th(bvzd_p2, bvzdca, xc_p2, yc_p2, zc_p2)
00455 call interpo_gr2cgr_4th(axx_p2 , axx , xc_p2, yc_p2, zc_p2)
00456 call interpo_gr2cgr_4th(axy_p2 , axy , xc_p2, yc_p2, zc_p2)
00457 call interpo_gr2cgr_4th(axz_p2 , axz , xc_p2, yc_p2, zc_p2)
00458 call interpo_gr2cgr_4th(ayy_p2 , ayy , xc_p2, yc_p2, zc_p2)
00459 call interpo_gr2cgr_4th(ayz_p2 , ayz , xc_p2, yc_p2, zc_p2)
00460 call interpo_gr2cgr_4th(azz_p2 , azz , xc_p2, yc_p2, zc_p2)
00461 call interpo_fl2cgr_4th_export(emd_p2 , emdca , xc_p2, yc_p2, zc_p2, rs_p2)
00462
00463 call interpo_fl2cgr_4th_export(vepxf_p2, vepxfca , xc_p2, yc_p2, zc_p2, rs_p2)
00464 call interpo_fl2cgr_4th_export(vepyf_p2, vepyfca , xc_p2, yc_p2, zc_p2, rs_p2)
00465 call interpo_fl2cgr_4th_export(vepzf_p2, vepzfca , xc_p2, yc_p2, zc_p2, rs_p2)
00466
00467 call interpo_fl2cgr_4th_export(wxspf_p2, wxspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00468 call interpo_fl2cgr_4th_export(wyspf_p2, wyspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00469 call interpo_fl2cgr_4th_export(wzspf_p2, wzspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00470
00471 call interpo_fl2cgr_4th_export(psif_p2 , psifca , xc_p2, yc_p2, zc_p2, rs_p2)
00472 call interpo_fl2cgr_4th_export(alphf_p2, alphfca , xc_p2, yc_p2, zc_p2, rs_p2)
00473 call interpo_fl2cgr_4th_export(bvxdf_p2, bvxdfca , xc_p2, yc_p2, zc_p2, rs_p2)
00474 call interpo_fl2cgr_4th_export(bvydf_p2, bvydfca , xc_p2, yc_p2, zc_p2, rs_p2)
00475 call interpo_fl2cgr_4th_export(bvzdf_p2, bvzdfca , xc_p2, yc_p2, zc_p2, rs_p2)
00476 bxcor = bvxdfca + ome_p2*(-ycoc)
00477 bycor = bvydfca + ome_p2*(xcoc)
00478 bzcor = bvzdfca
00479 lam_p2 = ber_p2 + bxcor*vepxfca + bycor*vepyfca + bzcor*vepzfca
00480 psi4ca = psica**4
00481 psif4ca = psifca**4
00482 psifcacp= psifca**confpow
00483
00484 if (dabs(emdca) > 1.0d-14) then
00485 vxu = ico*( ome_p2*(-ycoc) ) + &
00486 & irr*( alphfca*vepxfca/psif4ca/lam_p2 ) + &
00487 & isp*( alphfca*(vepxfca/psif4ca + psifcacp*wxspfca)/lam_p2 )
00488
00489 vyu = ico*( ome_p2*( xcoc) ) + &
00490 & irr*( alphfca*vepyfca/psif4ca/lam_p2 ) + &
00491 & isp*( alphfca*(vepyfca/psif4ca + psifcacp*wyspfca)/lam_p2 )
00492
00493 vzu = ico*( 0.0d0 ) + &
00494 & irr*( alphfca*vepzfca/psif4ca/lam_p2 ) + &
00495 & isp*( alphfca*(vepzfca/psif4ca + psifcacp*wzspfca)/lam_p2 )
00496 else
00497 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00498 end if
00499 end if
00500 end if
00501
00502 gxx = psi4ca
00503 gxy = 0.0d0
00504 gxz = 0.0d0
00505 gyy = psi4ca
00506 gyz = 0.0d0
00507 gzz = psi4ca
00508
00509 kxx = psi4ca*axx/(radi_p1)
00510 kxy = psi4ca*axy/(radi_p1)
00511 kxz = kij_parity*psi4ca*axz/(radi_p1)
00512 kyy = psi4ca*ayy/(radi_p1)
00513 kyz = kij_parity*psi4ca*ayz/(radi_p1)
00514 kzz = psi4ca*azz/(radi_p1)
00515
00516 call peos_q2hprho(emdca, hca, preca, rhoca, eneca)
00517
00518 epsca = eneca/rhoca - 1.0d0
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546 write(6,'(21e16.8)') xcac, psica, alphca, bvxdca, bvydca, bvzdca, preca, rhoca, &
00547 & epsca, vxu, vyu, vzu, vepxfca, vepyfca, vepzfca, kxx, kxy, kxz, kyy, kyz, kzz
00548 end do
00549 close(860)
00550
00551 write(6,'(a16)') "Deallocating...."
00552 deallocate( emd_p1); deallocate( emd_p2);
00553 deallocate( vep_p1); deallocate( vep_p2);
00554 deallocate(wxspf_p1); deallocate(wxspf_p2);
00555 deallocate(wyspf_p1); deallocate(wyspf_p2);
00556 deallocate(wzspf_p1); deallocate(wzspf_p2);
00557 deallocate(vepxf_p1); deallocate(vepxf_p2);
00558 deallocate(vepyf_p1); deallocate(vepyf_p2);
00559 deallocate(vepzf_p1); deallocate(vepzf_p2);
00560 deallocate( psif_p1); deallocate( psif_p2);
00561 deallocate(alphf_p1); deallocate(alphf_p2);
00562 deallocate(bvxdf_p1); deallocate(bvxdf_p2);
00563 deallocate(bvydf_p1); deallocate(bvydf_p2);
00564 deallocate(bvzdf_p1); deallocate(bvzdf_p2);
00565 deallocate( rs_p1); deallocate( rs_p2);
00566 deallocate( psi_p1); deallocate( psi_p2); deallocate(psi_p3)
00567 deallocate( alph_p1); deallocate( alph_p2); deallocate(alph_p3)
00568 deallocate( bvxd_p1); deallocate( bvxd_p2); deallocate(bvxd_p3)
00569 deallocate( bvyd_p1); deallocate( bvyd_p2); deallocate(bvyd_p3)
00570 deallocate( bvzd_p1); deallocate( bvzd_p2); deallocate(bvzd_p3)
00571 deallocate( axx_p1); deallocate( axx_p2); deallocate(axx_p3)
00572 deallocate( axy_p1); deallocate( axy_p2); deallocate(axy_p3)
00573 deallocate( axz_p1); deallocate( axz_p2); deallocate(axz_p3)
00574 deallocate( ayy_p1); deallocate( ayy_p2); deallocate(ayy_p3)
00575 deallocate( ayz_p1); deallocate( ayz_p2); deallocate(ayz_p3)
00576 deallocate( azz_p1); deallocate( azz_p2); deallocate(azz_p3)
00577
00578
00579 END PROGRAM coc2cac
