Added support for using all the different kinds of initial guesses for the

level set function for all the different level set functions.


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

1 files changed, 28 insertions, 28 deletions

diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 8fc9587..c291376 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -42,10 +42,10 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
     deallocate ( reparam_this_level_set )
   end if
   allocate ( reparam_this_level_set(eh_number_level_sets) )
-  if ( allocated(reparam_undone) ) then
-    deallocate ( reparam_undone )
+  if ( allocated(re_initialize_undone) ) then
+    deallocate ( re_initialize_undone )
   end if
-  allocate ( reparam_undone(eh_number_level_sets) )
+  allocate ( re_initialize_undone(eh_number_level_sets) )
 
   if ( evolve_generators .gt. 0 ) then
 
@@ -94,16 +94,16 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
       ! Set up a sphere of radius initial_rad and translated 
       ! (translate_x,translate_y,translate_z) away from the origin.
-      f(:,:,:,l) = sqrt( ( x - translate_x )**2 + &
-                         ( y - translate_y )**2 + &
-                         ( z - translate_z )**2 ) - initial_rad(l)
+      f(:,:,:,l) = sqrt( ( x - translate_x(l) )**2 + &
+                         ( y - translate_y(l) )**2 + &
+                         ( z - translate_z(l) )**2 ) - initial_rad(l)
       if ( evolve_generators .gt. 0 ) then
         if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
           do i = 1, lsh(1)
              theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-             xg(i,l) = initial_rad(l) * sin(theta) + translate_x
-             yg(i,l) = translate_y
-             zg(i,l) = initial_rad(l) * cos(theta) + translate_z
+             xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l)
+             yg(i,l) = translate_y(l)
+             zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l)
           end do
         end if
       end if
@@ -113,12 +113,12 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
     if ( CCTK_EQUALS( initial_f(l), 'ellipsoid' ) ) then
 
       ! Calculate sines and cosines of the rotation parameters.
-      cosa = cos(rotation_alpha)
-      sina = sin(rotation_alpha)
-      cosb = cos(rotation_beta)
-      sinb = sin(rotation_beta)
-      cosc = cos(rotation_gamma)
-      sinc = sin(rotation_gamma)
+      cosa = cos(rotation_alpha(l))
+      sina = sin(rotation_alpha(l))
+      cosb = cos(rotation_beta(l))
+      sinb = sin(rotation_beta(l))
+      cosc = cos(rotation_gamma(l))
+      sinc = sin(rotation_gamma(l))
 
       ! Set up the rotation matrix. The order is alpha around the z-axis,
       ! beta around the y-axis and finally gamma around the x-axis.
@@ -139,13 +139,13 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       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
+            xp(1) = x(i,j,k) - translate_x(l)
+            xp(2) = y(i,j,k) - translate_y(l)
+            xp(3) = z(i,j,k) - translate_z(l)
             xpt = matmul ( txyz, xp )
-            f(i,j,k,l) = sqrt( xpt(1)**2 / initial_a**2 + &
-                             xpt(2)**2 / initial_b**2 + &
-                             xpt(3)**2 / initial_c**2) - 1.0
+            f(i,j,k,l) = sqrt( xpt(1)**2 / initial_a(l)**2 + &
+                             xpt(2)**2 / initial_b(l)**2 + &
+                             xpt(3)**2 / initial_c(l)**2) - 1.0
           end do
         end do
       end do
@@ -154,11 +154,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
         if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
           do i = 1, lsh(1)
             theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-            r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + &
-                           cos(theta)**2 / initial_c**2 ) )
-            xp(1) = r_el * sin(theta) + translate_x
-            xp(2) = translate_y
-            xp(3) = r_el * cos(theta) + translate_z
+            r_el = sqrt ( one / ( sin(theta)**2 / initial_a(l)**2 + &
+                           cos(theta)**2 / initial_c(l)**2 ) )
+            xp(1) = r_el * sin(theta) + translate_x(l)
+            xp(2) = translate_y(l)
+            xp(3) = r_el * cos(theta) + translate_z(l)
             xpt = matmul ( txyz, xp )
             xg(i,l) = xpt(1)
             yg(i,l) = xpt(2)
@@ -170,8 +170,8 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
     ! if an ovaloid of Cassini is requested...
     if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
-      f(:,:,:,l) = (x**2+y**2+z**2)**2 + cas_a**4 - &
-            2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4
+      f(:,:,:,l) = (x**2+y**2+z**2)**2 + cas_a(l)**4 - &
+            2*cas_a(l)**2*(x**2 - (y**2+z**2)) - cas_b(l)**4
     end if
   end do