code to Lorentz-boost g_ab and g^ab by an arbitrary 3-vector velocity

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@166 e296648e-0e4f-0410-bd07-d597d9acff87

1 files changed, 327 insertions, 0 deletions

diff --git a/src/boost.F77 b/src/boost.F77
new file mode 100644
index 0000000..1016305
--- /dev/null
+++ b/src/boost.F77
@@ -0,0 +1,327 @@
+c This subroutine calculates the 4-metric and its inverse at an event,
+c taking into account an optional Lorentz boost.
+c $Header$
+c
+c The coordinates are
+c   Cx(a) = Cactus $x^a$
+c   Mx(a) = Model  $X^a$
+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 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. 
+c
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include "param_defs.inc"
+
+	subroutine Exact__metric(
+     $     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, rama)
+
+	implicit none
+	DECLARE_CCTK_FUNCTIONS
+	DECLARE_CCTK_PARAMETERS
+
+c input arguments
+	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, rama
+
+c intrinsic functions called
+	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
+
+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)
+
+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
+
+c flags, array indices, etc
+	logical Tmunu_flag
+	integer i, j
+	integer Ca, Cb, MA, MB
+
+c constants
+	integer n
+	parameter (n = 3)
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c optimized fast-path if no Lorentz boost
+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,
+     $			rama)
+		return
+	endif
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c the rest of this function is the Lorentz-boost case:
+c - Lorentz-transform Cactus coordinates --> Model coordinates
+c - compute Model 4-metric and inverse at Model coordinates
+c - tensor-transform 4-metric and inverse from Model coordinates
+c   --> Cactus coordinates
+c
+c All the equations used are given in ../doc/documentation.tex
+c
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c compute Lorentz transformation information on first call
+c
+	if (firstcall) then
+		firstcall = .false.
+
+c boost velocity
+		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
+
+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
+
+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
+
+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
+
+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
+
+c
+c Lorentz-transform the Cactus coordinate 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
+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,
+     $			rama)
+
+	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
+
+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)
+
+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
+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