diff options
Diffstat (limited to 'src/boost.F77')
| -rw-r--r-- | src/boost.F77 | 577 |
1 files changed, 357 insertions, 220 deletions
diff --git a/src/boost.F77 b/src/boost.F77 index d81b100..dc26340 100644 --- a/src/boost.F77 +++ b/src/boost.F77 @@ -1,5 +1,7 @@ c This subroutine calculates the 4-metric and its inverse at an event, -c taking into account an optional Lorentz boost. +c taking into account an optional Lorentz boost and an optional rotation. +c The model is first rotated and then boosted, such that the boost is +c applied to the rotated model. c $Header$ c c The coordinates are @@ -9,6 +11,12 @@ c The 4-metrics are c Cgdd(a,b) = Cactus $g_{ab}$ Cguu(a,b) = Cactus $g^{ab}$ c Mgdd(a,b) = Model $g_{ab}$ Mguu(a,b) = Model $g^{ab}$ c +c For a definition of the Euler angles in the conventions used below, see +c http://mathworld.wolfram.com/EulerAngles.html +c Another useful resource may be +c http://en.wikipedia.org/wiki/Euler_angles +c although this uses (on 2006-11-29) different conventions. +c c This file is copyright (c) 2003 by Jonathan Thornburg <jthorn@aei.mpg.de>. c This file is covered by the GNU GPL license; see the files ../README c and ../COPYING for details. @@ -21,7 +29,7 @@ c #include "param_defs.inc" - subroutine Exact__metric( + subroutine Exact__metric( $ decoded_exact_model, $ x, y, z, t, $ gdtt, gdtx, gdty, gdtz, @@ -30,83 +38,94 @@ c $ guxx, guyy, guzz, guxy, guyz, guxz, $ psi) - implicit none - DECLARE_CCTK_FUNCTIONS - DECLARE_CCTK_PARAMETERS + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS c input arguments - CCTK_INT decoded_exact_model - CCTK_REAL x, y, z, t + CCTK_INT decoded_exact_model + CCTK_REAL x, y, z, t c output arguments - CCTK_REAL gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, - $ psi + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi c intrinsic functions called - CCTK_REAL sqrt + CCTK_REAL sqrt c static local variables describing Lorentz transformation - logical firstcall - data firstcall /.true./ - CCTK_REAL gamma - CCTK_REAL vv(3), nn(3) - CCTK_REAL parallel(3,3), perp(3,3) - CCTK_REAL Cx_par(3), Cx_perp(3) - CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3) - CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3) - save firstcall - save gamma - save vv, nn - save parallel, perp - save Cx_par, Cx_perp - save partial_Mx_wrt_Cx - save partial_Cx_wrt_Mx + logical firstcall + data firstcall /.true./ + CCTK_REAL gamma + CCTK_REAL vv(3), nn(3) + CCTK_REAL parallel(3,3), perp(3,3) + CCTK_REAL Cx_par(3), Cx_perp(3) + CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3) + CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3) + CCTK_REAL R(0:3,0:3) + save firstcall + save gamma + save vv, nn + save parallel, perp + save Cx_par, Cx_perp + save partial_Mx_wrt_Cx + save partial_Cx_wrt_Mx + save R c coordinates and 4-metric - CCTK_REAL Ct, Cx(3) - CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3) - CCTK_REAL Mt, Mx(3) - CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3) + CCTK_REAL Cx(0:3) + CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3) + CCTK_REAL Mx(0:3) + CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3) + CCTK_REAL Nx(0:3) + CCTK_REAL Ngdd(0:3,0:3), Nguu(0:3,0:3) c misc temps - CCTK_REAL vnorm, vnormsq - CCTK_REAL delta_ij - CCTK_REAL Cx_par_i, Cx_perp_i - CCTK_REAL vdotCx - CCTK_REAL Cgdd_ab, Cguu_ab - character*100 warn_buffer + CCTK_REAL vnorm, vnormsq + CCTK_REAL delta_ij + CCTK_REAL Cx_par_i, Cx_perp_i + CCTK_REAL vdotCx + CCTK_REAL Cgdd_ab, Cguu_ab + CCTK_REAL cos_phi, sin_phi + CCTK_REAL cos_theta, sin_theta + CCTK_REAL cos_psi, sin_psi + CCTK_REAL R_phi(0:3,0:3), R_theta(0:3,0:3), R_psi(0:3,0:3) + character*1000 warn_buffer c flags, array indices, etc - logical Tmunu_flag - integer i, j - integer Ca, Cb, MA, MB + logical Tmunu_flag + integer i, j, k, l + integer Ca, Cb, MA, MB c constants - integer n - parameter (n = 3) + integer n + parameter (n = 3) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c -c optimized fast-path if no Lorentz boost +c optimized fast-path if no Lorentz boost and no rotation c - if ( (boost_vx .eq. 0.0) - $ .and. (boost_vy .eq. 0.0) - $ .and. (boost_vz .eq. 0.0) ) then - call Exact__metric_for_model( - $ decoded_exact_model, - $ x, y, z, t, - $ gdtt, gdtx, gdty, gdtz, - $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, - $ gutt, gutx, guty, gutz, - $ guxx, guyy, guzz, guxy, guyz, guxz, - $ psi, - $ Tmunu_flag) - return - endif + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0) + $ .and. (rotation_euler_phi .eq. 0.0) + $ .and. (rotation_euler_theta .eq. 0.0) + $ .and. (rotation_euler_psi .eq. 0.0)) then + call Exact__metric_for_model( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi, + $ Tmunu_flag) + return + end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -125,203 +144,321 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c compute Lorentz transformation information on first call c - if (firstcall) then - firstcall = .false. + if (firstcall) then + firstcall = .false. c boost velocity - vv(1) = boost_vx - vv(2) = boost_vy - vv(3) = boost_vz + vv(1) = boost_vx + vv(2) = boost_vy + vv(3) = boost_vz c Lorentz gamma factor, unit vector in direction of boost velocity - vnormsq = 0.0d0 - do 100 i = 1,n - vnormsq = vnormsq + vv(i)*vv(i) -100 continue - gamma = 1.0 / sqrt(1.0 - vnormsq) - vnorm = sqrt(vnormsq) - do 110 i = 1,n - nn(i) = vv(i) / vnorm -110 continue + vnormsq = 0.0 + do i = 1,n + vnormsq = vnormsq + vv(i)*vv(i) + end do + gamma = 1.0 / sqrt(1.0 - vnormsq) + vnorm = sqrt(vnormsq) + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0)) then + nn(1) = 1 + nn(2) = 0 + nn(3) = 0 + else + do i = 1,n + nn(i) = vv(i) / vnorm + end do + end if c projection operators parallel(*,*) and perp(*,*) - do 210 j = 1,n - do 200 i = 1,n - parallel(i,j) = nn(i) * nn(j) - if (i .eq. j) then - delta_ij = 1.0d0 - else - delta_ij = 0.0d0 - endif - perp(i,j) = delta_ij - parallel(i,j) -200 continue -210 continue + do j = 1,n + do i = 1,n + parallel(i,j) = nn(i) * nn(j) + if (i .eq. j) then + delta_ij = 1.0 + else + delta_ij = 0.0 + end if + perp(i,j) = delta_ij - parallel(i,j) + end do + end do c partial derivatives of Model coordinates with respect to Cactus coordinates - partial_Mx_wrt_Cx(0,0) = gamma - do 300 i = 1,n - partial_Mx_wrt_Cx(0,i) = -gamma*vv(i) -300 continue - do 320 i = 1,n - partial_Mx_wrt_Cx(i,0) = -gamma*vv(i) - do 310 j=1,n - partial_Mx_wrt_Cx(i,j) - $ = gamma*parallel(i,j) - $ + perp(i,j) -310 continue -320 continue + partial_Mx_wrt_Cx(0,0) = gamma + do i = 1,n + partial_Mx_wrt_Cx(0,i) = -gamma*vv(i) + end do + do i = 1,n + partial_Mx_wrt_Cx(i,0) = -gamma*vv(i) + do j=1,n + partial_Mx_wrt_Cx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do c partial derivatives of Cactus coordinates with respect to Model coordinates - partial_Cx_wrt_Mx(0,0) = gamma - do 400 i = 1,n - partial_Cx_wrt_Mx(0,i) = + gamma*vv(i) -400 continue - do 420 i = 1,n - partial_Cx_wrt_Mx(i,0) = + gamma*vv(i) - do 410 j=1,n - partial_Cx_wrt_Mx(i,j) - $ = gamma*parallel(i,j) - $ + perp(i,j) -410 continue -420 continue - - endif + partial_Cx_wrt_Mx(0,0) = gamma + do i = 1,n + partial_Cx_wrt_Mx(0,i) = + gamma*vv(i) + end do + do i = 1,n + partial_Cx_wrt_Mx(i,0) = + gamma*vv(i) + do j=1,n + partial_Cx_wrt_Mx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do + +c Sines and cosines of rotation angles + cos_phi = cos (rotation_euler_phi) + sin_phi = sin (rotation_euler_phi) + cos_theta = cos (rotation_euler_theta) + sin_theta = sin (rotation_euler_theta) + cos_psi = cos (rotation_euler_psi) + sin_psi = sin (rotation_euler_psi) + +c Set up individual rotation matrices + R_phi(0,0) = 1 + R_phi(0,1) = 0 + R_phi(0,2) = 0 + R_phi(0,3) = 0 + R_phi(1,0) = 0 + R_phi(1,1) = + cos_phi + R_phi(1,2) = + sin_phi + R_phi(1,3) = 0 + R_phi(2,0) = 0 + R_phi(2,1) = - sin_phi + R_phi(2,2) = + cos_phi + R_phi(2,3) = 0 + R_phi(3,0) = 0 + R_phi(3,1) = 0 + R_phi(3,2) = 0 + R_phi(3,3) = 1 + + R_theta(0,0) = 1 + R_theta(0,1) = 0 + R_theta(0,2) = 0 + R_theta(0,3) = 0 + R_theta(1,0) = 0 + R_theta(1,1) = 1 + R_theta(1,2) = 0 + R_theta(1,3) = 0 + R_theta(2,0) = 0 + R_theta(2,1) = 0 + R_theta(2,2) = + cos_theta + R_theta(2,3) = + sin_theta + R_theta(3,0) = 0 + R_theta(3,1) = 0 + R_theta(3,2) = - sin_theta + R_theta(3,3) = + cos_theta + + R_psi(0,0) = 1 + R_psi(0,1) = 0 + R_psi(0,2) = 0 + R_psi(0,3) = 0 + R_psi(1,0) = 0 + R_psi(1,1) = + cos_psi + R_psi(1,2) = + sin_psi + R_psi(1,3) = 0 + R_psi(2,0) = 0 + R_psi(2,1) = - sin_psi + R_psi(2,2) = + cos_psi + R_psi(2,3) = 0 + R_psi(3,0) = 0 + R_psi(3,1) = 0 + R_psi(3,2) = 0 + R_psi(3,3) = 1 + +c Combine individual rotation matrices + do i = 0,n + do j = 0,n + R(i,j) = 0 + do k = 0,n + do l = 0,n + R(i,j) = R(i,j) + R_psi(i,k) * R_theta(k,l) * R_phi(l,j) + end do + end do + end do + end do + +c Notes that help me (Erik Schnetter) think: +c This considers a rotation with phi=0, theta=pi/2, psi=0. +c Nx(i) = Nx(i) + R(j,i) * Mx(j) +c Nx(1) = - Mx(3) +c Nx(3) = Mx(1) +c Mgxx(1,1) = Ngxx(3,3) +c Mgxx(1,3) = - Ngxx(1,3) +c Mgxx(3,3) = Ngxx(1,1) +c Mgxx(i,j) = R(i,k) R(j,l) Ngxx(k,l) [correct] +c Mgxx(i,j) = R(k,i) R(l,j) Ngxx(k,l) [correct] +c Mbetax(1) = - Nbetax(3) +c Mbetax(3) = Nbetax(1) +c Mbetax(i) + R(i,j) Nbetax(j) [wrong] +c Mbetax(i) + R(j,i) Nbetax(j) [correct] + + end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c compute flat-space components of Cx(*) parallel and perpendicular to vv(*) c - Ct = t - Cx(1) = x - Cx(2) = y - Cx(3) = z - - do 530 i=1,n - Cx_par_i = 0.0d0 - Cx_perp_i = 0.0d0 - do 520 j=1,n - Cx_par_i = Cx_par_i - $ + parallel(i,j)*Cx(j) - Cx_perp_i = Cx_perp_i - $ + perp (i,j)*Cx(j) -520 continue - Cx_par (i) = Cx_par_i - Cx_perp(i) = Cx_perp_i -530 continue + Cx(0) = t + Cx(1) = x + Cx(2) = y + Cx(3) = z + + do i=1,n + Cx_par_i = 0.0 + Cx_perp_i = 0.0 + do j=1,n + Cx_par_i = Cx_par_i + parallel(i,j)*Cx(j) + Cx_perp_i = Cx_perp_i + perp (i,j)*Cx(j) + end do + Cx_par (i) = Cx_par_i + Cx_perp(i) = Cx_perp_i + end do c -c Lorentz-transform the Cactus coordinate to get the Model coordinates +c Lorentz-transform and rotate the Cactus coordinate +c to get the Model coordinates c - vdotCx = 0.0 - do 600 i = 1,n - vdotCx = vdotCx + vv(i)*Cx(i) -600 continue - Mt = gamma * (Ct - vdotCx) - do 610 i=1,n - Mx(i) = gamma * (Cx_par(i) - vv(i)*Ct) - $ + Cx_perp(i) -610 continue +c Boost + vdotCx = 0.0 + do i = 1,n + vdotCx = vdotCx + vv(i)*Cx(i) + end do + + Mx(0) = gamma * (Cx(0) - vdotCx) + do i=1,n + Mx(i) = gamma * (Cx_par(i) - vv(i)*Cx(0)) + Cx_perp(i) + end do + +c Rotation + do i=0,n + Nx(i) = 0 + do j = 0,n + Nx(i) = Nx(i) + R(i,j) * Mx(j) + end do + end do + c c compute the Model 4-metric and inverse 4-metric at the Model coordinates c - call Exact__metric_for_model( - $ decoded_exact_model, - $ Mx(1), Mx(2), Mx(3), Mt, - $ Mgdd(0,0), Mgdd(0,1), Mgdd(0,2), Mgdd(0,3), - $ Mgdd(1,1), Mgdd(2,2), Mgdd(3,3), - $ Mgdd(1,2), Mgdd(2,3), Mgdd(1,3), - $ Mguu(0,0), Mguu(0,1), Mguu(0,2), Mguu(0,3), - $ Mguu(1,1), Mguu(2,2), Mguu(3,3), - $ Mguu(1,2), Mguu(2,3), Mguu(1,3), - $ psi, Tmunu_flag) - - if (Tmunu_flag) then - write (warn_buffer, '(a,i8,a,a)') - $ 'exact_model = ', decoded_exact_model, - $ 'sets the stress-energy tensor', - $ ' ==> we cannot Lorentz-boost it! :(' - call CCTK_WARN(0, warn_buffer) - endif + call Exact__metric_for_model( + $ decoded_exact_model, + $ Nx(1), Nx(2), Nx(3), Nx(0), + $ Ngdd(0,0), Ngdd(0,1), Ngdd(0,2), Ngdd(0,3), + $ Ngdd(1,1), Ngdd(2,2), Ngdd(3,3), + $ Ngdd(1,2), Ngdd(2,3), Ngdd(1,3), + $ Nguu(0,0), Nguu(0,1), Nguu(0,2), Nguu(0,3), + $ Nguu(1,1), Nguu(2,2), Nguu(3,3), + $ Nguu(1,2), Nguu(2,3), Nguu(1,3), + $ psi, Tmunu_flag) + + if (Tmunu_flag) then + write (warn_buffer, '(a,i8,a,a)') + $ 'exact_model = ', decoded_exact_model, + $ 'sets the stress-energy tensor', + $ ' ==> we cannot Lorentz-boost or rotate it!' + call CCTK_WARN(0, warn_buffer) + end if c c symmetrize the Model 4-metric and inverse 4-metric arrays c (the Exact__metric_for_model() call only set the upper triangles) c - Mgdd(1,0) = Mgdd(0,1) - Mgdd(2,0) = Mgdd(0,2) - Mgdd(2,1) = Mgdd(1,2) - Mgdd(3,0) = Mgdd(0,3) - Mgdd(3,1) = Mgdd(1,3) - Mgdd(3,2) = Mgdd(2,3) - - Mguu(1,0) = Mguu(0,1) - Mguu(2,0) = Mguu(0,2) - Mguu(2,1) = Mguu(1,2) - Mguu(3,0) = Mguu(0,3) - Mguu(3,1) = Mguu(1,3) - Mguu(3,2) = Mguu(2,3) + Ngdd(1,0) = Ngdd(0,1) + Ngdd(2,0) = Ngdd(0,2) + Ngdd(2,1) = Ngdd(1,2) + Ngdd(3,0) = Ngdd(0,3) + Ngdd(3,1) = Ngdd(1,3) + Ngdd(3,2) = Ngdd(2,3) + + Nguu(1,0) = Nguu(0,1) + Nguu(2,0) = Nguu(0,2) + Nguu(2,1) = Nguu(1,2) + Nguu(3,0) = Nguu(0,3) + Nguu(3,1) = Nguu(1,3) + Nguu(3,2) = Nguu(2,3) c c tensor-transorm (the upper triangle of) the 4-metric and inverse 4-metric c - do 730 Ca = 0,n - do 720 Cb = Ca,n - Cgdd_ab = 0.0d0 - do 710 Ma = 0,n - do 700 Mb = 0,n - Cgdd_ab = Cgdd_ab - $ + Mgdd(Ma,Mb) - $ * partial_Mx_wrt_Cx(Ma,Ca) - $ * partial_Mx_wrt_Cx(Mb,Cb) -700 continue -710 continue - Cgdd(Ca,Cb) = Cgdd_ab -720 continue -730 continue - - do 830 Ca = 0,n - do 820 Cb = Ca,n - Cguu_ab = 0.0d0 - do 810 Ma = 0,n - do 800 Mb = 0,n - Cguu_ab = Cguu_ab - $ + Mguu(Ma,Mb) - $ * partial_Cx_wrt_Mx(Ca,Ma) - $ * partial_Cx_wrt_Mx(Cb,Mb) -800 continue -810 continue - Cguu(Ca,Cb) = Cguu_ab -820 continue -830 continue + +c Rotations + do i = 0,n + do j = 0,n + Mgdd(i,j) = 0 + Mguu(i,j) = 0 + do k = 0,n + do l = 0,n + Mgdd(i,j) = Mgdd(i,j) + R(k,i) * R(l,j) * Ngdd(k,l) +c The inverse of R is also its transpose. That means that the +c transpose of the inverse, which you would use for g^ij, is just R +c again. + Mguu(i,j) = Mguu(i,j) + R(k,i) * R(l,j) * Nguu(k,l) + end do + end do + end do + end do + +c Boost + do Ca = 0,n + do Cb = Ca,n + Cgdd_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cgdd_ab = Cgdd_ab + $ + Mgdd(Ma,Mb) + $ * partial_Mx_wrt_Cx(Ma,Ca) + $ * partial_Mx_wrt_Cx(Mb,Cb) + end do + end do + Cgdd(Ca,Cb) = Cgdd_ab + end do + end do + + do Ca = 0,n + do Cb = Ca,n + Cguu_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cguu_ab = Cguu_ab + $ + Mguu(Ma,Mb) + $ * partial_Cx_wrt_Mx(Ca,Ma) + $ * partial_Cx_wrt_Mx(Cb,Mb) + end do + end do + Cguu(Ca,Cb) = Cguu_ab + end do + end do c c unpack the Cactus-coordinates 4-metric and inverse 4-metric c into the corresponding output arguments c - gdtt = Cgdd(0,0) - gdtx = Cgdd(0,1) - gdty = Cgdd(0,2) - gdtz = Cgdd(0,3) - gdxx = Cgdd(1,1) - gdxy = Cgdd(1,2) - gdxz = Cgdd(1,3) - gdyy = Cgdd(2,2) - gdyz = Cgdd(2,3) - gdzz = Cgdd(3,3) - - gutt = Cguu(0,0) - gutx = Cguu(0,1) - guty = Cguu(0,2) - gutz = Cguu(0,3) - guxx = Cguu(1,1) - guxy = Cguu(1,2) - guxz = Cguu(1,3) - guyy = Cguu(2,2) - guyz = Cguu(2,3) - guzz = Cguu(3,3) - - return - end + gdtt = Cgdd(0,0) + gdtx = Cgdd(0,1) + gdty = Cgdd(0,2) + gdtz = Cgdd(0,3) + gdxx = Cgdd(1,1) + gdxy = Cgdd(1,2) + gdxz = Cgdd(1,3) + gdyy = Cgdd(2,2) + gdyz = Cgdd(2,3) + gdzz = Cgdd(3,3) + + gutt = Cguu(0,0) + gutx = Cguu(0,1) + guty = Cguu(0,2) + gutz = Cguu(0,3) + guxx = Cguu(1,1) + guxy = Cguu(1,2) + guxz = Cguu(1,3) + guyy = Cguu(2,2) + guyz = Cguu(2,3) + guzz = Cguu(3,3) + + end |