Added more explanatory comments in the routine EHFinder_Init_Weights.

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

1 files changed, 29 insertions, 0 deletions

diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index 24ad3e9..f1211d4 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -203,6 +203,9 @@ subroutine EHFinder_Integrate2(CCTK_ARGUMENTS)
   
 end subroutine EHFinder_Integrate2
 
+
+! This routine sets up the weights for the Simpsons rule integration
+! over the surface.
 subroutine EHFinder_Init_Weights(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -218,17 +221,30 @@ subroutine EHFinder_Init_Weights(CCTK_ARGUMENTS)
   CCTK_INT, dimension(4) :: bbox
   CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
 
+  ! Find out the lower bounds of the distributed integration grid arrays.
   call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, "ehfinder::surface_arrays" )
   if ( ierr .lt. 0 ) then
     call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" )
   end if
+
+  ! Find out the local size of the distributed integration grid arrays
   call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, "ehfinder::surface_arrays" )
   if ( ierr .lt. 0 ) then
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
 
   weights = one
+
+  ! Initialise the weight grid array for the 2D Simpsons rule integration.
+  ! To do this I need to figure out the global location of the given point.
+  ! There are 3 cases in the one dimensional case. If the point is on the
+  !  boundary the weight is 1/3. If it is at an even position the weight
+  ! is 4/3 and if it is at an odd position the weight is 2/3. 
+
   do j = 1, lsh(2)
+
+    ! This is first done in the phi direction. Meaning that all points with
+    ! the same theta coordinate are set to the same weight.
     if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then
       weights(:,j) = onethird
     else if ( mod(lbnd(2)+j,2) .eq. 0 ) then
@@ -236,6 +252,9 @@ subroutine EHFinder_Init_Weights(CCTK_ARGUMENTS)
     else
       weights(:,j) = twothirds
     end if
+
+    ! Then it is done in the theta direction with the one-dimensional
+    ! weights beeing multiplied.
     do i = 1, lsh(1)
       if ( ( lbnd(1)+i .eq. 1 ) .or. ( lbnd(1)+i .eq. ntheta ) ) then
         weights(i,j) = onethird * weights(i,j)
@@ -246,4 +265,14 @@ subroutine EHFinder_Init_Weights(CCTK_ARGUMENTS)
       end if
     end do
   end do
+
+  ! The end result is a 2D array with the weights in the following pattern.
+
+  !  1/9   4/9  2/9   4/9  2/9   4/9  1/9
+  !  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+  !  2/9   8/9  4/9   8/9  4/9   8/9  2/9
+  !  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+  !  2/9   8/9  4/9   8/9  4/9   8/9  2/9
+  !  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+  !  1/9   4/9  2/9   4/9  2/9   4/9  1/9
 end subroutine EHFinder_Init_Weights