Added parameters and code to rotate the initial ellipsoidal data for f and

and to translate the origin for spherical and ellipsoidal initial data.

The rotation is done in sequence first around the z-axis, then
the rotated y-axis and finally around the (twice) rotated x-axis.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@66 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 39 insertions, 6 deletions

diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 74103f1..172880d 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -17,6 +17,9 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
 
   CCTK_INT :: i, j, k
   CCTK_REAL :: last_time
+  CCTK_REAL, dimension(3) :: xp, xpt
+  CCTK_REAL, dimension(3,3) :: txyz
+  CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc
 
   ex_value = - ( one + shell_width ) * maxval(cctk_delta_space)
 
@@ -26,11 +29,45 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
   
   cctk_time = last_time
 
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+
   if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then
-    f = sqrt(x**2+y**2+z**2) - initial_rad ! + &
+    f = sqrt( ( x - translate_x )**2 + &
+              ( y - translate_y )**2 + &
+              ( z - translate_z )**2 ) - initial_rad
   end if
   if ( CCTK_EQUALS( initial_f, 'ellipsoid' ) ) then
-    f = sqrt(x**2/initial_a**2+y**2/initial_b**2+z**2/initial_c**2) - 1.0
+    cosa = cos(rotation_alpha)
+    sina = sin(rotation_alpha)
+    cosb = cos(rotation_beta)
+    sinb = sin(rotation_beta)
+    cosc = cos(rotation_gamma)
+    sinc = sin(rotation_gamma)
+    txyz(1,1) = cosa * cosb
+    txyz(1,2) = sina * cosb
+    txyz(1,3) = -sinb
+    txyz(2,1) = cosa * sinb * sinc - sina * cosc
+    txyz(2,2) = sina * sinb * sinc + cosa * cosc
+    txyz(2,3) = cosb * sinc
+    txyz(3,1) = cosa * sinb * cosc + sina * sinc
+    txyz(3,2) = sina * sinb * cosc - cosa * sinc
+    txyz(3,3) = cosb * cosc
+    print*,txyz
+    do k = 1, nz
+      do j = 1, ny
+        do i = 1, nx
+          xp(1) = x(i,j,k) - translate_x
+          xp(2) = y(i,j,k) - translate_y
+          xp(3) = z(i,j,k) - translate_z
+          xpt = matmul ( txyz, xp )
+          f(i,j,k) = sqrt( xpt(1)**2 / initial_a**2 + &
+                           xpt(2)**2 / initial_b**2 + &
+                           xpt(3)**2 / initial_c**2) - 1.0
+        end do
+      end do
+    end do
   end if
   if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
     f = (x**2+y**2+z**2)**2 + cas_a**4 - &
@@ -40,10 +77,6 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
   eh_mask = 0
   delta = maxval ( cctk_delta_space )
 
-  nx = cctk_lsh(1)
-  ny = cctk_lsh(2)
-  nz = cctk_lsh(3)
-
   ngx = cctk_nghostzones(1)
   ngy = cctk_nghostzones(2)
   ngz = cctk_nghostzones(3)