c Mahc_Init_Data: Various Initial Data Routines for Mahc_Evolve c Copyright (C) 2001 Mark Miller /*@@ @file boost_isons_pt.F @date January 2000 @author Mark Miller @desc Calculates the boosted NS solution pointwise from static TOV. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" /*@@ @routine boost_isons_pt @date January 2000 @author Mark Miller @desc This routine serves as a wrapper to BOOST_ISONS_PT_NOK, and adds in the Lie derivative of 3-metric wrt shift term to the extrinsic curvature. This routine can now be treated pointwise. @enddesc @calls @calledby @history @endhistory @@*/ c #define TEST subroutine boost_isons_pt(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime, . eos_gamma,eos_k,press,rho,eps,gxx,gxy,gxz,gyy,gyz,gzz,alp, . betax,betay,betaz,velx,vely,velz,kxx,kxy,kxz,kyy,kyz,kzz,dx, . mahc_init_eos_handle) implicit none integer ntov,imaxtov,mahc_init_eos_handle CCTK_REAL presstov(ntov),rhotov(ntov),epstov(ntov),mtov(ntov) CCTK_REAL phitov(ntov),r_schwartov(ntov),r_isotov(ntov) CCTK_REAL boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime CCTK_REAL press,rho,eps,eos_gamma,eos_k,gxx,gxy,gxz,gyy,gyz,gzz,alp CCTK_REAL betax,betay,betaz,velx,vely,velz CCTK_REAL kxx,kxy,kxz,kyy,kyz,kzz,dx,rmaxtov c LOCAL VARS CCTK_REAL gxx_px,gxy_px,gxz_px,gyy_px,gyz_px,gzz_px CCTK_REAL gxx_py,gxy_py,gxz_py,gyy_py,gyz_py,gzz_py CCTK_REAL gxx_pz,gxy_pz,gxz_pz,gyy_pz,gyz_pz,gzz_pz CCTK_REAL gxx_mx,gxy_mx,gxz_mx,gyy_mx,gyz_mx,gzz_mx CCTK_REAL gxx_my,gxy_my,gxz_my,gyy_my,gyz_my,gzz_my CCTK_REAL gxx_mz,gxy_mz,gxz_mz,gyy_mz,gyz_mz,gzz_mz CCTK_REAL betax_px,betay_px,betaz_px,betax_mx,betay_mx,betaz_mx CCTK_REAL betax_py,betay_py,betaz_py,betax_my,betay_my,betaz_my CCTK_REAL betax_pz,betay_pz,betaz_pz,betax_mz,betay_mz,betaz_mz CCTK_REAL press_temp,rho_temp,eps_temp,alp_temp,velx_temp,vely_temp CCTK_REAL velz_temp,kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp CCTK_REAL kzz_temp,h c metric derivatives CCTK_REAL dx_gxx,dx_gxy,dx_gxz,dx_gyy,dx_gyz,dx_gzz CCTK_REAL dy_gxx,dy_gxy,dy_gxz,dy_gyy,dy_gyz,dy_gzz CCTK_REAL dz_gxx,dz_gxy,dz_gxz,dz_gyy,dz_gyz,dz_gzz CCTK_REAL dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz CCTK_REAL dz_betax,dz_betay,dz_betaz CCTK_REAL liexx,liexy,liexz,lieyy,lieyz,liezz c get the basics... call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime, . eos_gamma,eos_k,press,rho,eps,gxx,gxy,gxz,gyy,gyz,gzz,alp, . betax,betay,betaz,velx,vely,velz,kxx,kxy,kxz,kyy,kyz,kzz, . mahc_init_eos_handle) c now, we need to add in the L_g beta term. So, get the c non-local stuff: h = dx/10.0d0 call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime+h,y_prime,z_prime, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_px,gxy_px,gxz_px,gyy_px,gyz_px,gzz_px,alp_temp, . betax_px,betay_px,betaz_px,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime-h,y_prime,z_prime, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_mx,gxy_mx,gxz_mx,gyy_mx,gyz_mx,gzz_mx,alp_temp, . betax_mx,betay_mx,betaz_mx,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime+h,z_prime, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_py,gxy_py,gxz_py,gyy_py,gyz_py,gzz_py,alp_temp, . betax_py,betay_py,betaz_py,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime-h,z_prime, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_my,gxy_my,gxz_my,gyy_my,gyz_my,gzz_my,alp_temp, . betax_my,betay_my,betaz_my,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime+h, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_pz,gxy_pz,gxz_pz,gyy_pz,gyz_pz,gzz_pz,alp_temp, . betax_pz,betay_pz,betaz_pz,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) call boost_isons_pt_noK(presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime-h, . eos_gamma,eos_k,press_temp,rho_temp,eps_temp, . gxx_mz,gxy_mz,gxz_mz,gyy_mz,gyz_mz,gzz_mz,alp_temp, . betax_mz,betay_mz,betaz_mz,velx_temp,vely_temp,velz_temp, . kxx_temp,kxy_temp,kxz_temp,kyy_temp,kyz_temp,kzz_temp, . mahc_init_eos_handle) c calculate the derivative terms needed. dx_gxx = (gxx_px - gxx_mx)/(2.0D0*h) dx_gxy = (gxy_px - gxy_mx)/(2.0D0*h) dx_gxz = (gxz_px - gxz_mx)/(2.0D0*h) dx_gyy = (gyy_px - gyy_mx)/(2.0D0*h) dx_gyz = (gyz_px - gyz_mx)/(2.0D0*h) dx_gzz = (gzz_px - gzz_mx)/(2.0D0*h) dy_gxx = (gxx_py - gxx_my)/(2.0D0*h) dy_gxy = (gxy_py - gxy_my)/(2.0D0*h) dy_gxz = (gxz_py - gxz_my)/(2.0D0*h) dy_gyy = (gyy_py - gyy_my)/(2.0D0*h) dy_gyz = (gyz_py - gyz_my)/(2.0D0*h) dy_gzz = (gzz_py - gzz_my)/(2.0D0*h) dz_gxx = (gxx_pz - gxx_mz)/(2.0D0*h) dz_gxy = (gxy_pz - gxy_mz)/(2.0D0*h) dz_gxz = (gxz_pz - gxz_mz)/(2.0D0*h) dz_gyy = (gyy_pz - gyy_mz)/(2.0D0*h) dz_gyz = (gyz_pz - gyz_mz)/(2.0D0*h) dz_gzz = (gzz_pz - gzz_mz)/(2.0D0*h) dx_betax = (betax_px - betax_mx)/(2.0D0*h) dx_betay = (betay_px - betay_mx)/(2.0D0*h) dx_betaz = (betaz_px - betaz_mx)/(2.0D0*h) dy_betax = (betax_py - betax_my)/(2.0D0*h) dy_betay = (betay_py - betay_my)/(2.0D0*h) dy_betaz = (betaz_py - betaz_my)/(2.0D0*h) dz_betax = (betax_pz - betax_mz)/(2.0D0*h) dz_betay = (betay_pz - betay_mz)/(2.0D0*h) dz_betaz = (betaz_pz - betaz_mz)/(2.0D0*h) liexx = . betax*dx_gxx + betay*dy_gxx + betaz*dz_gxx + . gxx*dx_betax + gxy*dx_betay + gxz*dx_betaz + . gxx*dx_betax + gxy*dx_betay + gxz*dx_betaz liexy = . betax*dx_gxy + betay*dy_gxy + betaz*dz_gxy + . gxy*dx_betax + gyy*dx_betay + gyz*dx_betaz + . gxx*dy_betax + gxy*dy_betay + gxz*dy_betaz liexz = . betax*dx_gxz + betay*dy_gxz + betaz*dz_gxz + . gxz*dx_betax + gyz*dx_betay + gzz*dx_betaz + . gxx*dz_betax + gxy*dz_betay + gxz*dz_betaz lieyy = . betax*dx_gyy + betay*dy_gyy + betaz*dz_gyy + . gxy*dy_betax + gyy*dy_betay + gyz*dy_betaz + . gxy*dy_betax + gyy*dy_betay + gyz*dy_betaz lieyz = . betax*dx_gyz + betay*dy_gyz + betaz*dz_gyz + . gxz*dy_betax + gyz*dy_betay + gzz*dy_betaz + . gxy*dz_betax + gyy*dz_betay + gyz*dz_betaz liezz = . betax*dx_gzz + betay*dy_gzz + betaz*dz_gzz + . gxz*dz_betax + gyz*dz_betay + gzz*dz_betaz + . gxz*dz_betax + gyz*dz_betay + gzz*dz_betaz kxx = kxx + liexx/(2.0d0 * alp) kxy = kxy + liexy/(2.0d0 * alp) kxz = kxz + liexz/(2.0d0 * alp) kyy = kyy + lieyy/(2.0d0 * alp) kyz = kyz + lieyz/(2.0d0 * alp) kzz = kzz + liezz/(2.0d0 * alp) return end /*@@ @routine boost_isons_pt @date January 2000 @author Mark Miller @desc This routine takes the tov solution, boost parameters, and input primed coordinates (moving frame coordinates), and returns the boosted TOV solution at that point. (Warning: the extrinsic curvature returned is only the term with the time derivative of the 3-metric. To properly reconstruct the extrinsic curvature, you will have to add the Lie derivative of 3-metric wrt shift term.) @enddesc @calls @calledby @history @endhistory @@*/ subroutine boost_isons_pt_noK(presstov,rhotov,epstov, . mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov, . boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime, . eos_gamma,eos_k,press,rho,eps,gxx,gxy,gxz,gyy,gyz,gzz,alp, . betax,betay,betaz,velx,vely,velz,kxx,kxy,kxz,kyy,kyz,kzz, . mahc_init_eos_handle) implicit none integer ntov,imaxtov,mahc_init_eos_handle CCTK_REAL presstov(ntov),rhotov(ntov),epstov(ntov),mtov(ntov) CCTK_REAL phitov(ntov),r_schwartov(ntov),r_isotov(ntov) CCTK_REAL boost_x,boost_y,boost_z,t_prime,x_prime,y_prime,z_prime CCTK_REAL press,rho,eps,eos_gamma,eos_k,gxx,gxy,gxz,gyy,gyz,gzz,alp CCTK_REAL betax,betay,betaz,velx,vely,velz CCTK_REAL kxx,kxy,kxz,kyy,kyz,kzz,rmaxtov c LOCAL VARS CCTK_REAL fourmetric_up(4,4),lam_p_up(4,4),lam_up_p(4,4),boost_gamma CCTK_REAL one,test,coord(4),coord_prime(4),zero,two CCTK_REAL fourmetric(4,4) integer i,j,k,l,m,ilow,ihigh c unprimed metric components CCTK_REAL gxx_up,gxy_up,gxz_up,gyy_up,gyz_up,gzz_up CCTK_REAL dxgxx_up,dxgxy_up,dxgxz_up,dxgyy_up,dxgyz_up,dxgzz_up CCTK_REAL dygxx_up,dygxy_up,dygxz_up,dygyy_up,dygyz_up,dygzz_up CCTK_REAL dzgxx_up,dzgxy_up,dzgxz_up,dzgyy_up,dzgyz_up,dzgzz_up CCTK_REAL alp_up,ax_up,ay_up,az_up CCTK_REAL betalowx,betalowy,betalowz,beta2 CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz CCTK_REAL fourvelocity_up(4),fourvelocity(4),w,v2 CCTK_REAL d_alp(3),d_g(3,3,3),boost_vec(3),qdot(3,3),temp CCTK_REAL tiny one = 1.0d0 zero = 0.0d0 two = 2.0d0 tiny = 1.0d-30 boost_gamma = one/sqrt(one - boost_x**2 - boost_y**2 - boost_z**2) c construct coordinate transformation matrices lam_p_up(1,1) = boost_gamma lam_p_up(1,2) = boost_x * boost_gamma lam_p_up(1,3) = boost_y * boost_gamma lam_p_up(1,4) = boost_z * boost_gamma lam_p_up(2,1) = lam_p_up(1,2) lam_p_up(3,1) = lam_p_up(1,3) lam_p_up(4,1) = lam_p_up(1,4) lam_p_up(2,2) = one + (boost_gamma - one) * boost_x**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(2,3) = (boost_gamma - one) * boost_x * boost_y / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(2,4) = (boost_gamma - one) * boost_x * boost_z / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(3,3) = one + (boost_gamma - one) * boost_y**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(3,4) = (boost_gamma - one) * boost_y * boost_z / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(4,4) = one + (boost_gamma - one) * boost_z**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_p_up(3,2) = lam_p_up(2,3) lam_p_up(4,2) = lam_p_up(2,4) lam_p_up(4,3) = lam_p_up(3,4) lam_up_p(1,1) = boost_gamma lam_up_p(1,2) = -boost_x * boost_gamma lam_up_p(1,3) = -boost_y * boost_gamma lam_up_p(1,4) = -boost_z * boost_gamma lam_up_p(2,1) = lam_up_p(1,2) lam_up_p(3,1) = lam_up_p(1,3) lam_up_p(4,1) = lam_up_p(1,4) lam_up_p(2,2) = one + (boost_gamma - one) * boost_x**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(2,3) = (boost_gamma - one) * boost_x * boost_y / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(2,4) = (boost_gamma - one) * boost_x * boost_z / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(3,3) = one + (boost_gamma - one) * boost_y**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(3,4) = (boost_gamma - one) * boost_y * boost_z / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(4,4) = one + (boost_gamma - one) * boost_z**2 / . (boost_x**2 + boost_y**2 + boost_z**2 + tiny) lam_up_p(3,2) = lam_up_p(2,3) lam_up_p(4,2) = lam_up_p(2,4) lam_up_p(4,3) = lam_up_p(3,4) #ifdef TEST write(*,*) 'check that coordinate transformations matrices' write(*,*) 'are consistent....' do i=1,4 do j=1,4 test = 0.0d0 do k = 1,4 test = test + lam_up_p(i,k) * lam_p_up(k,j) enddo write(*,*) ' ',i,j,test enddo enddo #endif c get unprimed (rest frame) coordinates coord_prime(1) = t_prime coord_prime(2) = x_prime coord_prime(3) = y_prime coord_prime(4) = z_prime do i=1,4 coord(i) = 0.0d0 do j=1,4 coord(i) = coord(i) + lam_up_p(i,j) * coord_prime(j) enddo enddo #ifdef TEST write(*,*) 'primed, unprimed coordinates: ' do i=1,4 write(*,*) coord_prime(i),coord(i) enddo #endif c get metric in unprimed (rest frame) coordinates call gettov_cart_iso(coord(2),coord(3),coord(4),press,rho,eps,gxx_up, . gxy_up,gxz_up,gyy_up,gyz_up,gzz_up, . dxgxx_up,dxgxy_up,dxgxz_up,dxgyy_up,dxgyz_up,dxgzz_up, . dygxx_up,dygxy_up,dygxz_up,dygyy_up,dygyz_up,dygzz_up, . dzgxx_up,dzgxy_up,dzgxz_up,dzgyy_up,dzgyz_up,dzgzz_up, . alp_up,ax_up,ay_up,az_up,eos_gamma,eos_k,zero, . presstov,rhotov,epstov,mtov,phitov,r_schwartov,r_isotov,ntov, . imaxtov,rmaxtov,mahc_init_eos_handle) fourmetric_up(1,1) = - alp_up**2 fourmetric_up(1,2) = 0.0d0 fourmetric_up(1,3) = 0.0d0 fourmetric_up(1,4) = 0.0d0 fourmetric_up(2,1) = 0.0d0 fourmetric_up(3,1) = 0.0d0 fourmetric_up(4,1) = 0.0d0 fourmetric_up(2,2) = gxx_up fourmetric_up(2,3) = gxy_up fourmetric_up(2,4) = gxz_up fourmetric_up(3,3) = gyy_up fourmetric_up(3,4) = gyz_up fourmetric_up(4,4) = gzz_up fourmetric_up(3,2) = fourmetric_up(2,3) fourmetric_up(4,2) = fourmetric_up(2,4) fourmetric_up(4,3) = fourmetric_up(3,4) c get primed metric do i=1,4 do j=1,4 fourmetric(i,j) = 0.0d0 do k=1,4 do l=1,4 fourmetric(i,j) = fourmetric(i,j) + . lam_up_p(k,i) * lam_up_p(l,j) * fourmetric_up(k,l) enddo enddo enddo enddo c return the metric gxx = fourmetric(2,2) gxy = fourmetric(2,3) gxz = fourmetric(2,4) gyy = fourmetric(3,3) gyz = fourmetric(3,4) gzz = fourmetric(4,4) c now, need to calculate lapse and shift betalowx = fourmetric(1,2) betalowy = fourmetric(1,3) betalowz = fourmetric(1,4) det= -(gxz**2*gyy) + two*gxy*gxz*gyz - gxx*gyz**2 & - gxy**2*gzz + gxx*gyy*gzz c invert metric uxx=(-gyz**2 + gyy*gzz)/det uxy=(gxz*gyz - gxy*gzz)/det uyy=(-gxz**2 + gxx*gzz)/det uxz=(-gxz*gyy + gxy*gyz)/det uyz=(gxy*gxz - gxx*gyz)/det uzz=(-gxy**2 + gxx*gyy)/det betax = uxx*betalowx + uxy*betalowy + uxz*betalowz betay = uxy*betalowx + uyy*betalowy + uyz*betalowz betaz = uxz*betalowx + uyz*betalowy + uzz*betalowz beta2 = betax*betalowx + betay*betalowy + betaz*betalowz c put square of lapse in alp alp = beta2 - fourmetric(1,1) if(alp < 0.0d0) then write(*,*) 'boost_ns_pt: error: alp not properly defined!' STOP endif alp = sqrt(alp) c now, calculate 4-velocity in each frame fourvelocity_up(1) = 1/alp_up fourvelocity_up(2) = 0.0d0 fourvelocity_up(3) = 0.0d0 fourvelocity_up(4) = 0.0d0 do i=1,4 fourvelocity(i) = 0.0d0 do j=1,4 fourvelocity(i) = fourvelocity(i) + . lam_p_up(i,j) * fourvelocity_up(j) enddo enddo #ifdef TEST c check normalization of u test = 0.0d0 do i=1,4 do j=1,4 test = test + fourvelocity(i) * fourvelocity(j) * fourmetric(i,j) enddo enddo write(*,*) 'norm of 4-velocity: ',test #endif c get lorentz factor (w) and primative velocities w = alp * fourvelocity(1) velx = (fourvelocity(2) + fourvelocity(1)*betax)/w vely = (fourvelocity(3) + fourvelocity(1)*betay)/w velz = (fourvelocity(4) + fourvelocity(1)*betaz)/w #ifdef TEST c check the consistency of veli with w v2 = gxx*velx**2 + gyy*vely**2 + gzz*velz**2 + . two*gxy*velx*vely + two*gxz*velx*velz + two*gyz*vely*velz write(*,*) 'comparing two lorentz factors: ',w, . one/sqrt(one - v2) #endif c now, get time derivative of the 3-metric, all in the primed c (moving frame) coordinates. All this is done in the c unprimed frame d_alp(1) = alp_up * ax_up d_alp(2) = alp_up * ay_up d_alp(3) = alp_up * az_up d_g(1,1,1) = two * dxgxx_up d_g(1,1,2) = two * dxgxy_up d_g(1,1,3) = two * dxgxz_up d_g(1,2,1) = two * dxgxy_up d_g(1,2,2) = two * dxgyy_up d_g(1,2,3) = two * dxgyz_up d_g(1,3,1) = two * dxgxz_up d_g(1,3,2) = two * dxgyz_up d_g(1,3,3) = two * dxgzz_up d_g(2,1,1) = two * dygxx_up d_g(2,1,2) = two * dygxy_up d_g(2,1,3) = two * dygxz_up d_g(2,2,1) = two * dygxy_up d_g(2,2,2) = two * dygyy_up d_g(2,2,3) = two * dygyz_up d_g(2,3,1) = two * dygxz_up d_g(2,3,2) = two * dygyz_up d_g(2,3,3) = two * dygzz_up d_g(3,1,1) = two * dzgxx_up d_g(3,1,2) = two * dzgxy_up d_g(3,1,3) = two * dzgxz_up d_g(3,2,1) = two * dzgxy_up d_g(3,2,2) = two * dzgyy_up d_g(3,2,3) = two * dzgyz_up d_g(3,3,1) = two * dzgxz_up d_g(3,3,2) = two * dzgyz_up d_g(3,3,3) = two * dzgzz_up boost_vec(1) = boost_x boost_vec(2) = boost_y boost_vec(3) = boost_z do i=1,3 do j=i,3 qdot(i,j) = 0.0d0 do k=1,3 temp = 0.0d0 do l=1,3 do m=1,3 temp = temp + . lam_up_p(l+1,i+1) * lam_up_p(m+1,j+1) * d_g(k,l,m) enddo enddo qdot(i,j) = qdot(i,j) - (boost_vec(k) * boost_gamma) * . (temp - two * alp_up * boost_gamma**2 * boost_vec(i) * . boost_vec(j) * d_alp(k)) enddo enddo enddo kxx = qdot(1,1)/(- two * alp) kxy = qdot(1,2)/(- two * alp) kxz = qdot(1,3)/(- two * alp) kyy = qdot(2,2)/(- two * alp) kyz = qdot(2,3)/(- two * alp) kzz = qdot(3,3)/(- two * alp) c if(abs(rcoord) .gt. rcenttov(ntov-2)) then c write(*,*) 'gettov:error: rcoord out of bounds. ' c write(*,*) ' ntov is ',ntov c write(*,*) ' rcenttov(ntov-2) ',rcenttov(ntov-2) c write(*,*) ' rcoord ',rcoord c write(*,*) ' abs(rcoord) ',abs(rcoord) c stop c endif c i = ( abs(rcoord) * dble(ntov - 1)/rmaxtov + 1.5d0) c ilow = i c if(rcenttov(ilow) .le. rcoord) then c ihigh = ilow+1 c if(rcenttov(ihigh) .lt. rcoord) then c write(*,*) 'gettov: alg error 1 ' c stop c endif c else c ihigh = ilow c ilow = ilow-1 c if(rcoord .lt. rcenttov(ilow)) then c write(*,*) 'gettov: alg error 2 ' c stop c endif c endif c if(i .eq. 1) then c rho = rhotov(i) c m = mtov(i) c phi = phitov(i) c else c rho = ( rhotov(ihigh)*(rcoord - rcenttov(ilow)) + c . rhotov(ilow) *(rcenttov(ihigh) - rcoord))/ c . (rcenttov(ihigh) - rcenttov(ilow)) c m = ( mtov(ihigh)*(rcoord - rcenttov(ilow)) + c . mtov(ilow) *(rcenttov(ihigh) - rcoord))/ c . (rcenttov(ihigh) - rcenttov(ilow)) c phi = ( phitov(ihigh)*(rcoord - rcenttov(ilow)) + c . phitov(ilow) *(rcenttov(ihigh) - rcoord))/ c . (rcenttov(ihigh) - rcenttov(ilow)) c endif return end