added radial shift drift-correction, changed 1.1 to cactus parameter

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@266 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 58 insertions, 12 deletions

diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 803c7da..063dd22 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -46,7 +46,7 @@
       CCTK_REAL trxi,xi2
       CCTK_REAL zero,half,one,two,three,four,pi
       CCTK_REAL aux,sina,cosa
-      CCTK_REAL omega_corr_all
+      CCTK_REAL omega_corr_all, radvel_corr_all
 
       CCTK_REAL, dimension(3,3)     :: ug,xi
       CCTK_REAL, dimension(2,2)     :: ga,ua
@@ -583,18 +583,17 @@
 
          ahf_centroid_x = avgx
          ahf_centroid_y = avgy
+         write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x
+         write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y
 
 !        Calculate correction in omega to stop "drift" of apparent horizon.
-
-         if (drift_correct_on.eq.1) then
-
-            if ((.not.find3).or.
+         if ((.not.find3).or.
      .          (find3.and.(mfind-1.eq.drift_correct_horizon))) then
 
+            if (drift_correct_on.eq.1) then
+
                xi_new = atan2(avgx,avgy)
 
-               write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x
-               write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y
                write(0,*) 'AHFinder_gau: xi_new = ',xi_new
 
                if (ahf_ncall.gt.1) then
@@ -603,7 +602,7 @@
 
                   write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
 
-                  omega_corr = - 1.1D0 * xi_drift
+                  omega_corr = - drift_correct_factor * xi_drift
      .               /(dble(ahf_findevery)*cctk_delta_time)
                   omega_cum  = omega_cum + omega_corr
 
@@ -622,12 +621,31 @@
 
                end if
 
-               write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum
+               write(0,*) 'AHFinder_gau: omega_cum = ',omega_cum
 
                xi_old = xi_new
 
             end if
-         end if
+
+!        Do drift correction in the radial direction
+           if (rad_drift_correct.eq.1) then
+
+                r_new = sqrt(avgx**2 + avgy**2)
+
+                if (ahf_ncall.gt.1) then
+
+                   radvel_corr = - rad_drift_correct_factor * 
+     .                (r_new - r_old) /(dble(ahf_findevery)*cctk_delta_time)
+
+                   write(0,*) 'AHFinder_gau: radvel_corr = ',radvel_corr
+                else
+                   radvel_corr = 0
+                endif
+
+                r_old = r_new
+
+            endif
+        endif
 
 !     Other processors.
 
@@ -643,14 +661,15 @@
 !        (Recommended by Thomas Radke)
 
          omega_corr = zero
+         radvel_corr = zero
 
       end if
 
 !     Actually apply drift correction to shift.
 
-      if (drift_correct_on.eq.1) then
-        if ((.not.find3).or.
+      if ((.not.find3).or.
      .       (find3.and.(mfind-1.eq.drift_correct_horizon))) then
+        if (drift_correct_on.eq.1) then
 
             call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
      .         omega_corr,omega_corr_all,CCTK_VARIABLE_REAL)
@@ -671,8 +690,35 @@
             end if
 
          end if
+
+        if (rad_drift_correct.eq.1) then
+
+            call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
+     .         radvel_corr,radvel_corr_all,CCTK_VARIABLE_REAL)
+            if (ierror.ne.0) then
+               call CCTK_WARN(1,"Reduction failed!")
+            end if
+
+            radvel_corr = radvel_corr_all
+
+            ahf_radvel_corr = radvel_corr
+
+            if (ahf_ncall.gt.1) then
+               
+               where (x**2+y**2 .gt. 0.001)
+                  betax = betax + radvel_corr * x/sqrt(x**2+y**2) /
+     .                     psi**rotation_psipower
+                  betay = betay + radvel_corr * y/sqrt(x**2+y**2) / 
+     .                     psi**rotation_psipower
+               end where
+            end if
+
+         end if
       end if
 
+
+
+
 !     Reduce the errors across processors (all processors must
 !     know about this since all will participate on the interpolation
 !     below).