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
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 write(6,'(a56)', ADVANCE = "NO") "Give cartesian coordinates (x,y,z) separated by a space:"
00296 read(5,*) xcac,ycac,zcac
00297 write(6,'(a23,3e20.12)') "Point given wrt CACTUS:", xcac,ycac,zcac
00298 write(6,'(a38,2e20.12)') "Cocal radius scale in COCP-1, COCP-2 :", radi_p1, radi_p2
00299 write(6,'(a38,2e20.12)') "Cocal surface scale in COCP-1, COCP-2:", r_surf_p1, r_surf_p2
00300 xcoc = xcac/(radi_p1)
00301 ycoc = ycac/(radi_p1)
00302 zcoc = zcac/(radi_p1)
00303 write(6,'(a23,3e20.12)') "Point given wrt COCAL:", xcoc,ycoc,zcoc
00304
00305 rrcm = sqrt(xcoc**2 + ycoc**2 + zcoc**2)
00306 write(6,*) rrcm, rr3
00307 if (rrcm > rr3) then
00308
00309 call copy_grid_parameter_from_mpt(3)
00310 call copy_grid_parameter_binary_excision_from_mpt(3)
00311 call copy_coordinate_grav_extended_from_mpt(3)
00312 call copy_coordinate_grav_phi_from_mpt(3)
00313 call copy_coordinate_grav_r_from_mpt(3)
00314 call copy_coordinate_grav_theta_from_mpt(3)
00315 call copy_def_binary_parameter_from_mpt(3)
00316 call copy_trigonometry_grav_phi_from_mpt(3)
00317 call copy_trigonometry_grav_theta_from_mpt(3)
00318
00319 kij_parity = 1
00320 xc_p3 = xcoc
00321 yc_p3 = ycoc
00322 zc_p3 = zcoc
00323 write(6,'(a23,3e20.12)') "Point given wrt COCP-3:", xc_p3,yc_p3,zc_p3
00324
00325 call interpo_gr2cgr_4th(psi_p3 , psica , xc_p3, yc_p3, zc_p3)
00326 call interpo_gr2cgr_4th(alph_p3, alphca, xc_p3, yc_p3, zc_p3)
00327 call interpo_gr2cgr_4th(bvxd_p3, bvxdca, xc_p3, yc_p3, zc_p3)
00328 call interpo_gr2cgr_4th(bvyd_p3, bvydca, xc_p3, yc_p3, zc_p3)
00329 call interpo_gr2cgr_4th(bvzd_p3, bvzdca, xc_p3, yc_p3, zc_p3)
00330 call interpo_gr2cgr_4th(axx_p3 , axx , xc_p3, yc_p3, zc_p3)
00331 call interpo_gr2cgr_4th(axy_p3 , axy , xc_p3, yc_p3, zc_p3)
00332 call interpo_gr2cgr_4th(axz_p3 , axz , xc_p3, yc_p3, zc_p3)
00333 call interpo_gr2cgr_4th(ayy_p3 , ayy , xc_p3, yc_p3, zc_p3)
00334 call interpo_gr2cgr_4th(ayz_p3 , ayz , xc_p3, yc_p3, zc_p3)
00335 call interpo_gr2cgr_4th(azz_p3 , azz , xc_p3, yc_p3, zc_p3)
00336 psi4ca = psica**4
00337 vxu = 0.0d0
00338 vyu = 0.0d0
00339 vzu = 0.0d0
00340 emdca = 0.0d0
00341 else
00342 if (xcoc<=0.0d0) then
00343
00344 call copy_grid_parameter_from_mpt(1)
00345 call copy_grid_parameter_binary_excision_from_mpt(1)
00346 call copy_coordinate_grav_extended_from_mpt(1)
00347 call copy_coordinate_grav_phi_from_mpt(1)
00348 call copy_coordinate_grav_r_from_mpt(1)
00349 call copy_coordinate_grav_theta_from_mpt(1)
00350 call copy_def_binary_parameter_from_mpt(1)
00351 call copy_trigonometry_grav_phi_from_mpt(1)
00352 call copy_trigonometry_grav_theta_from_mpt(1)
00353
00354 kij_parity = 1
00355 xc_p1 = xcoc + dis_cm
00356 yc_p1 = ycoc
00357 zc_p1 = zcoc
00358 write(6,'(a23,3e20.12)') "Point given wrt COCP-1:", xc_p1,yc_p1,zc_p1
00359 call interpo_gr2cgr_4th(psi_p1 , psica , xc_p1, yc_p1, zc_p1)
00360 call interpo_gr2cgr_4th(alph_p1, alphca, xc_p1, yc_p1, zc_p1)
00361 call interpo_gr2cgr_4th(bvxd_p1, bvxdca, xc_p1, yc_p1, zc_p1)
00362 call interpo_gr2cgr_4th(bvyd_p1, bvydca, xc_p1, yc_p1, zc_p1)
00363 call interpo_gr2cgr_4th(bvzd_p1, bvzdca, xc_p1, yc_p1, zc_p1)
00364 call interpo_gr2cgr_4th(axx_p1 , axx , xc_p1, yc_p1, zc_p1)
00365 call interpo_gr2cgr_4th(axy_p1 , axy , xc_p1, yc_p1, zc_p1)
00366 call interpo_gr2cgr_4th(axz_p1 , axz , xc_p1, yc_p1, zc_p1)
00367 call interpo_gr2cgr_4th(ayy_p1 , ayy , xc_p1, yc_p1, zc_p1)
00368 call interpo_gr2cgr_4th(ayz_p1 , ayz , xc_p1, yc_p1, zc_p1)
00369 call interpo_gr2cgr_4th(azz_p1 , azz , xc_p1, yc_p1, zc_p1)
00370 call interpo_fl2cgr_4th_export(emd_p1 , emdca , xc_p1, yc_p1, zc_p1, rs_p1)
00371
00372 call interpo_fl2cgr_4th_export(vepxf_p1, vepxfca , xc_p1, yc_p1, zc_p1, rs_p1)
00373 call interpo_fl2cgr_4th_export(vepyf_p1, vepyfca , xc_p1, yc_p1, zc_p1, rs_p1)
00374 call interpo_fl2cgr_4th_export(vepzf_p1, vepzfca , xc_p1, yc_p1, zc_p1, rs_p1)
00375
00376 call interpo_fl2cgr_4th_export(wxspf_p1, wxspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00377 call interpo_fl2cgr_4th_export(wyspf_p1, wyspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00378 call interpo_fl2cgr_4th_export(wzspf_p1, wzspfca , xc_p1, yc_p1, zc_p1, rs_p1)
00379
00380 call interpo_fl2cgr_4th_export(psif_p1 , psifca , xc_p1, yc_p1, zc_p1, rs_p1)
00381 call interpo_fl2cgr_4th_export(alphf_p1, alphfca , xc_p1, yc_p1, zc_p1, rs_p1)
00382 call interpo_fl2cgr_4th_export(bvxdf_p1, bvxdfca , xc_p1, yc_p1, zc_p1, rs_p1)
00383 call interpo_fl2cgr_4th_export(bvydf_p1, bvydfca , xc_p1, yc_p1, zc_p1, rs_p1)
00384 call interpo_fl2cgr_4th_export(bvzdf_p1, bvzdfca , xc_p1, yc_p1, zc_p1, rs_p1)
00385 bxcor = bvxdfca + ome_p1*(-ycoc)
00386 bycor = bvydfca + ome_p1*(xcoc)
00387 bzcor = bvzdfca
00388 lam_p1 = ber_p1 + bxcor*vepxfca + bycor*vepyfca + bzcor*vepzfca
00389 psi4ca = psica**4
00390 psif4ca = psifca**4
00391 psifcacp= psifca**confpow
00392
00393 if (dabs(emdca) > 1.0d-14) then
00394 vxu = ico*( ome_p1*(-ycoc) ) + &
00395 & irr*( alphfca*vepxfca/psif4ca/lam_p1 ) + &
00396 & isp*( alphfca*(vepxfca/psif4ca + psifcacp*wxspfca)/lam_p1 )
00397
00398 vyu = ico*( ome_p1*( xcoc) ) + &
00399 & irr*( alphfca*vepyfca/psif4ca/lam_p1 ) + &
00400 & isp*( alphfca*(vepyfca/psif4ca + psifcacp*wyspfca)/lam_p1 )
00401
00402 vzu = ico*( 0.0d0 ) + &
00403 & irr*( alphfca*vepzfca/psif4ca/lam_p1 ) + &
00404 & isp*( alphfca*(vepzfca/psif4ca + psifcacp*wzspfca)/lam_p1 )
00405 else
00406 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00407 end if
00408 else
00409
00410 call copy_grid_parameter_from_mpt(2)
00411 call copy_grid_parameter_binary_excision_from_mpt(2)
00412 call copy_coordinate_grav_extended_from_mpt(2)
00413 call copy_coordinate_grav_phi_from_mpt(2)
00414 call copy_coordinate_grav_r_from_mpt(2)
00415 call copy_coordinate_grav_theta_from_mpt(2)
00416 call copy_def_binary_parameter_from_mpt(2)
00417 call copy_trigonometry_grav_phi_from_mpt(2)
00418 call copy_trigonometry_grav_theta_from_mpt(2)
00419
00420 kij_parity = -1
00421 xc_p2 = -(xcoc - dis_cm)
00422 yc_p2 = -ycoc
00423 zc_p2 = zcoc
00424 write(6,'(a23,3e20.12)') "Point given wrt COCP-2:", xc_p2,yc_p2,zc_p2
00425 call interpo_gr2cgr_4th(psi_p2 , psica , xc_p2, yc_p2, zc_p2)
00426 call interpo_gr2cgr_4th(alph_p2, alphca, xc_p2, yc_p2, zc_p2)
00427 call interpo_gr2cgr_4th(bvxd_p2, bvxdca, xc_p2, yc_p2, zc_p2)
00428 call interpo_gr2cgr_4th(bvyd_p2, bvydca, xc_p2, yc_p2, zc_p2)
00429 call interpo_gr2cgr_4th(bvzd_p2, bvzdca, xc_p2, yc_p2, zc_p2)
00430 call interpo_gr2cgr_4th(axx_p2 , axx , xc_p2, yc_p2, zc_p2)
00431 call interpo_gr2cgr_4th(axy_p2 , axy , xc_p2, yc_p2, zc_p2)
00432 call interpo_gr2cgr_4th(axz_p2 , axz , xc_p2, yc_p2, zc_p2)
00433 call interpo_gr2cgr_4th(ayy_p2 , ayy , xc_p2, yc_p2, zc_p2)
00434 call interpo_gr2cgr_4th(ayz_p2 , ayz , xc_p2, yc_p2, zc_p2)
00435 call interpo_gr2cgr_4th(azz_p2 , azz , xc_p2, yc_p2, zc_p2)
00436 call interpo_fl2cgr_4th_export(emd_p2 , emdca , xc_p2, yc_p2, zc_p2, rs_p2)
00437
00438 call interpo_fl2cgr_4th_export(vepxf_p2, vepxfca , xc_p2, yc_p2, zc_p2, rs_p2)
00439 call interpo_fl2cgr_4th_export(vepyf_p2, vepyfca , xc_p2, yc_p2, zc_p2, rs_p2)
00440 call interpo_fl2cgr_4th_export(vepzf_p2, vepzfca , xc_p2, yc_p2, zc_p2, rs_p2)
00441
00442 call interpo_fl2cgr_4th_export(wxspf_p2, wxspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00443 call interpo_fl2cgr_4th_export(wyspf_p2, wyspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00444 call interpo_fl2cgr_4th_export(wzspf_p2, wzspfca , xc_p2, yc_p2, zc_p2, rs_p2)
00445
00446 call interpo_fl2cgr_4th_export(psif_p2 , psifca , xc_p2, yc_p2, zc_p2, rs_p2)
00447 call interpo_fl2cgr_4th_export(alphf_p2, alphfca , xc_p2, yc_p2, zc_p2, rs_p2)
00448 call interpo_fl2cgr_4th_export(bvxdf_p2, bvxdfca , xc_p2, yc_p2, zc_p2, rs_p2)
00449 call interpo_fl2cgr_4th_export(bvydf_p2, bvydfca , xc_p2, yc_p2, zc_p2, rs_p2)
00450 call interpo_fl2cgr_4th_export(bvzdf_p2, bvzdfca , xc_p2, yc_p2, zc_p2, rs_p2)
00451 bxcor = bvxdfca + ome_p2*(-ycoc)
00452 bycor = bvydfca + ome_p2*(xcoc)
00453 bzcor = bvzdfca
00454 lam_p2 = ber_p2 + bxcor*vepxfca + bycor*vepyfca + bzcor*vepzfca
00455 psi4ca = psica**4
00456 psif4ca = psifca**4
00457 psifcacp= psifca**confpow
00458
00459 if (dabs(emdca) > 1.0d-14) then
00460 vxu = ico*( ome_p2*(-ycoc) ) + &
00461 & irr*( alphfca*vepxfca/psif4ca/lam_p2 ) + &
00462 & isp*( alphfca*(vepxfca/psif4ca + psifcacp*wxspfca)/lam_p2 )
00463
00464 vyu = ico*( ome_p2*( xcoc) ) + &
00465 & irr*( alphfca*vepyfca/psif4ca/lam_p2 ) + &
00466 & isp*( alphfca*(vepyfca/psif4ca + psifcacp*wyspfca)/lam_p2 )
00467
00468 vzu = ico*( 0.0d0 ) + &
00469 & irr*( alphfca*vepzfca/psif4ca/lam_p2 ) + &
00470 & isp*( alphfca*(vepzfca/psif4ca + psifcacp*wzspfca)/lam_p2 )
00471 else
00472 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00473 end if
00474 end if
00475 end if
00476
00477 gxx = psi4ca
00478 gxy = 0.0d0
00479 gxz = 0.0d0
00480 gyy = psi4ca
00481 gyz = 0.0d0
00482 gzz = psi4ca
00483
00484 kxx = psi4ca*axx/(radi_p1)
00485 kxy = psi4ca*axy/(radi_p1)
00486 kxz = kij_parity*psi4ca*axz/(radi_p1)
00487 kyy = psi4ca*ayy/(radi_p1)
00488 kyz = kij_parity*psi4ca*ayz/(radi_p1)
00489 kzz = psi4ca*azz/(radi_p1)
00490
00491 call peos_q2hprho(emdca, hca, preca, rhoca, eneca)
00492
00493 epsca = eneca/rhoca - 1.0d0
00494
00495 write(6,'(a6,e20.12)') "psi =", psica
00496 write(6,'(a6,e20.12)') "alph =", alphca
00497 write(6,'(a6,e20.12)') "bvxd =", bvxdca
00498 write(6,'(a6,e20.12)') "bvyd =", bvydca
00499 write(6,'(a6,e20.12)') "bvzd =", bvzdca
00500 write(6,'(a6,e20.12)') "Radi =", r_surf_p1*radi_p1
00501 write(6,'(a6,e20.12)') "Omeg =", ome_p1/radi_p1
00502 write(6,'(a6,e20.12)') "emd =", emdca
00503 write(6,'(a6,e20.12)') "h =", hca
00504 write(6,'(a6,e20.12)') "pre =", preca
00505 write(6,'(a6,e20.12)') "rho =", rhoca
00506 write(6,'(a6,e20.12)') "ene =", eneca
00507 write(6,'(a6,e20.12)') "eps =", epsca
00508 write(6,'(a6,e20.12)') "vepx =", vepxfca
00509 write(6,'(a6,e20.12)') "vepy =", vepyfca
00510 write(6,'(a6,e20.12)') "vepz =", vepzfca
00511
00512 write(6,'(a18)') "Kij at gridpoints:"
00513 write(6,'(3e20.12)') kxx, kxy, kxz
00514 write(6,'(3e20.12)') kxy, kyy, kyz
00515 write(6,'(3e20.12)') kxz, kyz, kzz
00516
00517 write(6,'(a13)') "v^i eulerian:"
00518 write(6,'(a6,e20.12)') "vxu =", vxu
00519 write(6,'(a6,e20.12)') "vyu =", vyu
00520 write(6,'(a6,e20.12)') "vzu =", vzu
00521
00522 write(6,'(a16)') "Deallocating...."
00523 deallocate( emd_p1); deallocate( emd_p2);
00524 deallocate( vep_p1); deallocate( vep_p2);
00525 deallocate(wxspf_p1); deallocate(wxspf_p2);
00526 deallocate(wyspf_p1); deallocate(wyspf_p2);
00527 deallocate(wzspf_p1); deallocate(wzspf_p2);
00528 deallocate(vepxf_p1); deallocate(vepxf_p2);
00529 deallocate(vepyf_p1); deallocate(vepyf_p2);
00530 deallocate(vepzf_p1); deallocate(vepzf_p2);
00531 deallocate( psif_p1); deallocate( psif_p2);
00532 deallocate(alphf_p1); deallocate(alphf_p2);
00533 deallocate(bvxdf_p1); deallocate(bvxdf_p2);
00534 deallocate(bvydf_p1); deallocate(bvydf_p2);
00535 deallocate(bvzdf_p1); deallocate(bvzdf_p2);
00536 deallocate( rs_p1); deallocate( rs_p2);
00537 deallocate( psi_p1); deallocate( psi_p2); deallocate(psi_p3)
00538 deallocate( alph_p1); deallocate( alph_p2); deallocate(alph_p3)
00539 deallocate( bvxd_p1); deallocate( bvxd_p2); deallocate(bvxd_p3)
00540 deallocate( bvyd_p1); deallocate( bvyd_p2); deallocate(bvyd_p3)
00541 deallocate( bvzd_p1); deallocate( bvzd_p2); deallocate(bvzd_p3)
00542 deallocate( axx_p1); deallocate( axx_p2); deallocate(axx_p3)
00543 deallocate( axy_p1); deallocate( axy_p2); deallocate(axy_p3)
00544 deallocate( axz_p1); deallocate( axz_p2); deallocate(axz_p3)
00545 deallocate( ayy_p1); deallocate( ayy_p2); deallocate(ayy_p3)
00546 deallocate( ayz_p1); deallocate( ayz_p2); deallocate(ayz_p3)
00547 deallocate( azz_p1); deallocate( azz_p2); deallocate(azz_p3)
00548
00549
00550 END PROGRAM coc2cac
