Changing the code so that arryas start at 1 and not at 0. This is better

if we want to transform them later to Cactus arrays.


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

1 files changed, 124 insertions, 68 deletions

diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index d953a90..37be435 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -198,45 +198,45 @@
 
       if (myproc.eq.0) then
 
-         allocate(rr(0:ntheta,0:nphi))
+         allocate(rr(1:ntheta+1,1:nphi+1))
 
-         allocate(xa(0:ntheta,0:nphi))
-         allocate(ya(0:ntheta,0:nphi))
-         allocate(za(0:ntheta,0:nphi))
+         allocate(xa(1:ntheta+1,1:nphi+1))
+         allocate(ya(1:ntheta+1,1:nphi+1))
+         allocate(za(1:ntheta+1,1:nphi+1))
 
-         allocate(txx(0:ntheta,0:nphi))
-         allocate(tyy(0:ntheta,0:nphi))
-         allocate(tzz(0:ntheta,0:nphi))
-         allocate(txy(0:ntheta,0:nphi))
-         allocate(txz(0:ntheta,0:nphi))
-         allocate(tyz(0:ntheta,0:nphi))
+         allocate(txx(1:ntheta+1,1:nphi+1))
+         allocate(tyy(1:ntheta+1,1:nphi+1))
+         allocate(tzz(1:ntheta+1,1:nphi+1))
+         allocate(txy(1:ntheta+1,1:nphi+1))
+         allocate(txz(1:ntheta+1,1:nphi+1))
+         allocate(tyz(1:ntheta+1,1:nphi+1))
 
-         allocate(g11(0:ntheta,0:nphi))
-         allocate(g22(0:ntheta,0:nphi))
-         allocate(g12(0:ntheta,0:nphi))
+         allocate(g11(1:ntheta+1,1:nphi+1))
+         allocate(g22(1:ntheta+1,1:nphi+1))
+         allocate(g12(1:ntheta+1,1:nphi+1))
 
-         allocate(gaussian(0:ntheta,0:nphi))
+         allocate(gaussian(1:ntheta+1,1:nphi+1))
 
       else
 
-         allocate(rr(0:0,0:0))
+         allocate(rr(1,1))
 
-         allocate(xa(0:0,0:0))
-         allocate(ya(0:0,0:0))
-         allocate(za(0:0,0:0))
+         allocate(xa(1,1))
+         allocate(ya(1,1))
+         allocate(za(1,1))
 
-         allocate(txx(0:0,0:0))
-         allocate(tyy(0:0,0:0))
-         allocate(tzz(0:0,0:0))
-         allocate(txy(0:0,0:0))
-         allocate(txz(0:0,0:0))
-         allocate(tyz(0:0,0:0))
+         allocate(txx(1,1))
+         allocate(tyy(1,1))
+         allocate(tzz(1,1))
+         allocate(txy(1,1))
+         allocate(txz(1,1))
+         allocate(tyz(1,1))
 
-         allocate(g11(0:0,0:0))
-         allocate(g22(0:0,0:0))
-         allocate(g12(0:0,0:0))
+         allocate(g11(1,1))
+         allocate(g22(1,1))
+         allocate(g12(1,1))
 
-         allocate(gaussian(0:0,0:0))
+         allocate(gaussian(1,1))
 
       end if
 
@@ -491,13 +491,16 @@
 
       if (myproc.eq.0) then
 
-         do j=0,nphi
-            do i=0,ntheta
+         avgx = 0.0D0
+         avgy = 0.0D0
+
+         do j=1,nphi+1
+            do i=1,ntheta+1
 
 !              Find {theta,phi}.
 
-               theta = dtheta*dble(i)
-               phi   = dphi*dble(j) + phistart
+               theta = dtheta*dble(i-1)
+               phi   = dphi*dble(j-1) + phistart
 
 !              Find sines and cosines.
 
@@ -558,6 +561,9 @@
                xp = xa(i,j)
                yp = ya(i,j)
                zp = za(i,j)
+ 
+               avgx = avgx + xp
+               avgy = avgy + yp
 
                if ((xp.gt.xmx).or.(xp.lt.xmn).or.
      .             (yp.gt.ymx).or.(yp.lt.ymn).or.
@@ -568,6 +574,56 @@
             end do
          end do
 
+         avgx = avgx/(dble(nphi+1)*(ntheta+1))
+         avgy = avgy/(dble(nphi+1)*(ntheta+1))
+
+!        Alter shift to correct for rotational "drift" of apparent horizon.
+         
+         if (drift_correct_on.eq.1) then
+
+           write(0,*) '**** AHFinder_gau: avgx = ',avgx
+           write(0,*) '**** AHFinder_gau: avgy = ',avgy
+
+           xi_new = atan2(avgx,avgy)
+
+           if (ahf_ncall .gt. 1) then
+              xi_drift = xi_new - xi_old
+           else
+              xi_drift = xi_new
+           endif
+
+           omega_corr = - xi_drift / (ahf_findevery * cctk_delta_time)
+
+           if (ahf_ncall .gt. 1) then
+              omega_cum = omega_cum + omega_corr
+           else
+              omega_cum = omega_corr
+           endif
+
+           write(0,*) '**** AHFinder_gau: ncall = ',ahf_ncall
+           write(0,*) '**** AHFinder_gau: xi_drift = ',xi_drift
+           write(0,*) '**** AHFinder_gau: ahf_findevery = ',
+     .         ahf_findevery
+           write(0,*) '**** AHFinder_gau: cctk_delta_time = ',
+     .            cctk_delta_time
+           write(0,*) '**** AHFinder_gau: omega_corr = ',omega_corr
+           write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum
+
+           do k=1,nz
+             do j=1,ny
+               do i=1,nx
+                 betax(i,j,k) = betax(i,j,k) 
+     .                    - y(i,j,k)*omega_corr/psi(i,j,k)**2
+                 betay(i,j,k) = betay(i,j,k) 
+     .                    + x(i,j,k)*omega_corr/psi(i,j,k)**2
+               end do
+             end do
+           end do
+
+           xi_old = xi_new
+
+         end if
+
       else
 
          xa = half*(xmx+xmn)
@@ -679,13 +735,13 @@
 
       if (myproc.eq.0) then
 
-         do j=0,nphi
-            do i=0,ntheta
+         do j=1,nphi+1
+            do i=1,ntheta+1
 
 !              Find {theta,phi}.
 
-               theta = dtheta*dble(i)
-               phi   = dphi*dble(j) + phistart
+               theta = dtheta*dble(i-1)
+               phi   = dphi*dble(j-1) + phistart
 
 !              Find {rp}.
 
@@ -730,25 +786,25 @@
 !              use centered differences, and for the boundary points
 !              I use second order one-sided differences.
 
-               if ((i.ne.0).and.(i.ne.ntheta)) then
+               if ((i.ne.1).and.(i.ne.ntheta+1)) then
                   ft = idtheta*(rr(i+1,j) - rr(i-1,j))
-               else if (i.eq.0) then
-                  ft = - idtheta*(three*rr(0,j)
-     .               - four*rr(1,j) + rr(2,j))
+               else if (i.eq.1) then
+                  ft = - idtheta*(three*rr(1,j)
+     .               - four*rr(2,j) + rr(3,j))
                else
                   ft = + idtheta*(three*rr(ntheta,j)
-     .               - four*rr(ntheta-1,j) + rr(ntheta-2,j))
+     .               - four*rr(ntheta,j) + rr(ntheta-1,j))
                end if
 
                if (nonaxi) then
-                  if ((j.ne.0).and.(j.ne.nphi)) then
+                  if ((j.ne.1).and.(j.ne.nphi+1)) then
                      fp = idphi*(rr(i,j+1) - rr(i,j-1))
-                  else if (j.eq.0) then
-                     fp = - idphi*(three*rr(i,0)
-     .                  - four*rr(i,1) + rr(i,2))
+                  else if (j.eq.1) then
+                     fp = - idphi*(three*rr(i,1)
+     .                  - four*rr(i,2) + rr(i,3))
                   else
                      fp = + idphi*(three*rr(i,nphi)
-     .                  - four*rr(i,nphi-1) + rr(i,nphi-2))
+     .                  - four*rr(i,nphi) + rr(i,nphi-1))
                   end if
                else
                   fp = zero
@@ -776,9 +832,9 @@
 
          if (refz) then
 
-            do j=0,nphi-1
+            do j=1,nphi
                circ_eq = circ_eq + dphi
-     .            *dsqrt(dabs(half*(g22(ntheta,j) + g22(ntheta,j+1))))
+     .            *dsqrt(dabs(half*(g22(ntheta+1,j) + g22(ntheta+1,j+1))))
             end do
 
          else
@@ -790,7 +846,7 @@
             if (cartoon) then
                circ_eq = 6.283185307D0*dsqrt(g22(i,0))
             else
-               do j=0,nphi-1
+               do j=1,nphi
                   circ_eq = circ_eq + dphi
      .               *((one-aux)*dsqrt(dabs(half*(g22(i,j) + g22(i,j+1))))
      .               + aux*dsqrt(dabs(half*(g22(i+1,j) + g22(i+1,j+1)))))
@@ -800,23 +856,23 @@
 
 !        Meridians.
 
-         do i=0,ntheta-1
+         do i=1,ntheta
             meri_p1 = meri_p1 + dtheta
-     .         *dsqrt(dabs(half*(g11(i,0) + g11(i+1,0))))
+     .         *dsqrt(dabs(half*(g11(i,1) + g11(i+1,1))))
          end do
 
          if (cartoon) then
             meri_p2 = meri_p1
          else if (refx.and.refy) then
-            do i=0,ntheta-1
+            do i=1,ntheta
                meri_p2 = meri_p2 + dtheta
-     .            *dsqrt(dabs(half*(g11(i,nphi) + g11(i+1,nphi))))
+     .            *dsqrt(dabs(half*(g11(i,nphi+1) + g11(i+1,nphi+1))))
             end do
          else
             aux = half*pi/dphi
             j = int(aux)
             aux = aux - int(aux)
-            do i=0,ntheta-1
+            do i=1,ntheta
                meri_p2 = meri_p2 + dtheta
      .            *((one-aux)*dsqrt(dabs(half*(g11(i,j) + g11(i,j))))
      .            + aux*dsqrt(dabs(half*(g11(i,j+1) + g11(i,j+1)))))
@@ -847,7 +903,7 @@
 !     second derivatives of the interpolated metric, and this
 !     seems not to behave well numerically (we would need higher
 !     order interpolation).  I did not want to delete this code
-!     since it might be useful in the future, so I here just jump
+!     since it might be useful in the future, so here I just jump
 !     over it.
 
       if (.false.) then
@@ -882,8 +938,8 @@
 
 !        Interior points.
 
-         do j=1,nphi-1
-            do i=1,ntheta-1
+         do j=2,nphi
+            do i=2,ntheta
 
 !              Find angular metric.
 
@@ -977,7 +1033,7 @@
 !              first two terms of the Taylor series (assuming that
 !              for small theta  g_{phi,phi} = f(r,phi) sin(theta)^2).
 
-               theta = dtheta*dble(i)
+               theta = dtheta*dble(i-1)
 
                d1a(1,2,2) = d1a(1,2,2) + (two/three)*sin(theta)
      .                    *(g22(i+1,j) - two*g22(i,j) + g22(i-1,j))
@@ -1053,13 +1109,13 @@
 
 !        Boundaries.
 
-         gaussian(0,:) = two*gaussian(1,:) - gaussian(2,:)
-         gaussian(ntheta,:) = two*gaussian(ntheta-1,:)
-     .                      - gaussian(ntheta-2,:)
+         gaussian(1,:) = two*gaussian(2,:) - gaussian(3,:)
+         gaussian(ntheta+1,:) = two*gaussian(ntheta,:)
+     .                      - gaussian(ntheta-1,:)
 
-         gaussian(:,0) = two*gaussian(:,1) - gaussian(:,2)
-         gaussian(:,nphi) = two*gaussian(:,nphi-1)
-     .                    - gaussian(:,nphi-2)
+         gaussian(:,1) = two*gaussian(:,2) - gaussian(:,3)
+         gaussian(:,nphi+1) = two*gaussian(:,nphi)
+     .                    - gaussian(:,nphi-1)
 
       end if
 
@@ -1116,8 +1172,8 @@
             write(1,"(A1)")  '#'
             write(1,"(A35)") '# The data is written in a loop as:'
             write(1,"(A1)")  '#'
-            write(1,"(A17)") '# do i=0,ntheta-1'
-            write(1,"(A18)") '#    do j=0,nphi-1'
+            write(1,"(A17)") '# do i=1,ntheta'
+            write(1,"(A18)") '#    do j=1,nphi'
             write(1,"(A27)") '#       write gaussian(i,j)'
             write(1,"(A11)") '#    end do'
             write(1,"(A8)")  '# end do'
@@ -1161,8 +1217,8 @@
 
 !        Save data points.
 
-         do i=0,ntheta
-            do j=0,nphi
+         do i=1,ntheta+1
+            do j=1,nphi+1
                write(1,"(ES11.3)") gaussian(i,j)
             end do
          end do