Add optional argument "order" to set the order of differentiation.

Use an explicit path to include another file.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@39 b716e942-a2de-43ad-8f52-f3dfe468e4e7

2 files changed, 172 insertions, 119 deletions

diff --git a/src/pointwise.F90 b/src/pointwise.F90
index 6f435cb..3b55f0e 100644
--- a/src/pointwise.F90
+++ b/src/pointwise.F90
@@ -27,7 +27,7 @@ module pointwise
   public get_tensorderivs3
 contains
 #define TGR_INCLUDED
-#include "derivs.F90"
+#include "TAT/TGRtensor/src/derivs.F90"
 #undef TGR_INCLUDED
   subroutine calc_position (shape, i,j,k, pos)
     integer, intent(in)  :: shape(3)
@@ -166,38 +166,41 @@ contains
     gammazyz(pos) = gamma(3,2,3)
     gammazzz(pos) = gamma(3,3,3)
   end subroutine set_connections
-  subroutine get_scalarderivs (a, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: a(*)
-    CCTK_REAL, intent(out) :: f(3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
-    call get_derivs (a, f, pos, off, dx)
+  subroutine get_scalarderivs (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
+    call get_derivs (a, f, pos, off, dx, order)
   end subroutine get_scalarderivs
-  subroutine get_vectorderivs (ax,ay,az, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
-    CCTK_REAL, intent(out) :: f(3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_vectorderivs (ax,ay,az, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out)           :: f(3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fx(3),fy(3),fz(3)
-    call get_derivs (ax, fx, pos, off, dx)
-    call get_derivs (ay, fy, pos, off, dx)
-    call get_derivs (az, fz, pos, off, dx)
+    call get_derivs (ax, fx, pos, off, dx, order)
+    call get_derivs (ay, fy, pos, off, dx, order)
+    call get_derivs (az, fz, pos, off, dx, order)
     f(1,:) = fx
     f(2,:) = fy
     f(3,:) = fz
   end subroutine get_vectorderivs
-  subroutine get_tensorderivs (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
-    CCTK_REAL, intent(out) :: f(3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_tensorderivs (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fxx(3),fxy(3),fxz(3),fyy(3),fyz(3),fzz(3)
-    call get_derivs (axx, fxx, pos, off, dx)
-    call get_derivs (axy, fxy, pos, off, dx)
-    call get_derivs (axz, fxz, pos, off, dx)
-    call get_derivs (ayy, fyy, pos, off, dx)
-    call get_derivs (ayz, fyz, pos, off, dx)
-    call get_derivs (azz, fzz, pos, off, dx)
+    call get_derivs (axx, fxx, pos, off, dx, order)
+    call get_derivs (axy, fxy, pos, off, dx, order)
+    call get_derivs (axz, fxz, pos, off, dx, order)
+    call get_derivs (ayy, fyy, pos, off, dx, order)
+    call get_derivs (ayz, fyz, pos, off, dx, order)
+    call get_derivs (azz, fzz, pos, off, dx, order)
     f(1,1,:) = fxx
     f(1,2,:) = fxy
     f(1,3,:) = fxz
@@ -208,38 +211,41 @@ contains
     f(3,2,:) = fyz
     f(3,3,:) = fzz
   end subroutine get_tensorderivs
-  subroutine get_scalarderivs2 (a, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: a(*)
-    CCTK_REAL, intent(out) :: f(3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
-    call get_derivs2 (a, f, pos, off, dx)
+  subroutine get_scalarderivs2 (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
+    call get_derivs2 (a, f, pos, off, dx, order)
   end subroutine get_scalarderivs2
-  subroutine get_vectorderivs2 (ax,ay,az, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
-    CCTK_REAL, intent(out) :: f(3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_vectorderivs2 (ax,ay,az, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fx(3,3),fy(3,3),fz(3,3)
-    call get_derivs2 (ax, fx, pos, off, dx)
-    call get_derivs2 (ay, fy, pos, off, dx)
-    call get_derivs2 (az, fz, pos, off, dx)
+    call get_derivs2 (ax, fx, pos, off, dx, order)
+    call get_derivs2 (ay, fy, pos, off, dx, order)
+    call get_derivs2 (az, fz, pos, off, dx, order)
     f(1,:,:) = fx
     f(2,:,:) = fy
     f(3,:,:) = fz
   end subroutine get_vectorderivs2
-  subroutine get_tensorderivs2 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
-    CCTK_REAL, intent(out) :: f(3,3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_tensorderivs2 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fxx(3,3),fxy(3,3),fxz(3,3),fyy(3,3),fyz(3,3),fzz(3,3)
-    call get_derivs2 (axx, fxx, pos, off, dx)
-    call get_derivs2 (axy, fxy, pos, off, dx)
-    call get_derivs2 (axz, fxz, pos, off, dx)
-    call get_derivs2 (ayy, fyy, pos, off, dx)
-    call get_derivs2 (ayz, fyz, pos, off, dx)
-    call get_derivs2 (azz, fzz, pos, off, dx)
+    call get_derivs2 (axx, fxx, pos, off, dx, order)
+    call get_derivs2 (axy, fxy, pos, off, dx, order)
+    call get_derivs2 (axz, fxz, pos, off, dx, order)
+    call get_derivs2 (ayy, fyy, pos, off, dx, order)
+    call get_derivs2 (ayz, fyz, pos, off, dx, order)
+    call get_derivs2 (azz, fzz, pos, off, dx, order)
     f(1,1,:,:) = fxx
     f(1,2,:,:) = fxy
     f(1,3,:,:) = fxz
@@ -250,39 +256,42 @@ contains
     f(3,2,:,:) = fyz
     f(3,3,:,:) = fzz
   end subroutine get_tensorderivs2
-  subroutine get_scalarderivs3 (a, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: a(*)
-    CCTK_REAL, intent(out) :: f(3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
-    call get_derivs3 (a, f, pos, off, dx)
+  subroutine get_scalarderivs3 (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
+    call get_derivs3 (a, f, pos, off, dx, order)
   end subroutine get_scalarderivs3
-  subroutine get_vectorderivs3 (ax,ay,az, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
-    CCTK_REAL, intent(out) :: f(3,3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_vectorderivs3 (ax,ay,az, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fx(3,3,3),fy(3,3,3),fz(3,3,3)
-    call get_derivs3 (ax, fx, pos, off, dx)
-    call get_derivs3 (ay, fy, pos, off, dx)
-    call get_derivs3 (az, fz, pos, off, dx)
+    call get_derivs3 (ax, fx, pos, off, dx, order)
+    call get_derivs3 (ay, fy, pos, off, dx, order)
+    call get_derivs3 (az, fz, pos, off, dx, order)
     f(1,:,:,:) = fx
     f(2,:,:,:) = fy
     f(3,:,:,:) = fz
   end subroutine get_vectorderivs3
-  subroutine get_tensorderivs3 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
-    CCTK_REAL, intent(out) :: f(3,3,3,3,3)
-    integer,   intent(in)  :: pos, off(3)
-    CCTK_REAL, intent(in)  :: dx(3)
+  subroutine get_tensorderivs3 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out)           :: f(3,3,3,3,3)
+    integer,   intent(in)            :: pos, off(3)
+    CCTK_REAL, intent(in)            :: dx(3)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fxx(3,3,3),fxy(3,3,3),fxz(3,3,3),&
          &       fyy(3,3,3),fyz(3,3,3),fzz(3,3,3)
-    call get_derivs3 (axx, fxx, pos, off, dx)
-    call get_derivs3 (axy, fxy, pos, off, dx)
-    call get_derivs3 (axz, fxz, pos, off, dx)
-    call get_derivs3 (ayy, fyy, pos, off, dx)
-    call get_derivs3 (ayz, fyz, pos, off, dx)
-    call get_derivs3 (azz, fzz, pos, off, dx)
+    call get_derivs3 (axx, fxx, pos, off, dx, order)
+    call get_derivs3 (axy, fxy, pos, off, dx, order)
+    call get_derivs3 (axz, fxz, pos, off, dx, order)
+    call get_derivs3 (ayy, fyy, pos, off, dx, order)
+    call get_derivs3 (ayz, fyz, pos, off, dx, order)
+    call get_derivs3 (azz, fzz, pos, off, dx, order)
     f(1,1,:,:,:) = fxx
     f(1,2,:,:,:) = fxy
     f(1,3,:,:,:) = fxz
diff --git a/src/pointwise2.F90 b/src/pointwise2.F90
index 85e10c8..d337721 100644
--- a/src/pointwise2.F90
+++ b/src/pointwise2.F90
@@ -20,9 +20,12 @@ module pointwise2
   public get_2scalarderivs2
   public get_2vectorderivs2
   public get_2tensorderivs2
+  public get_2scalarderivs3
+  public get_2vectorderivs3
+  public get_2tensorderivs3
 contains
 #define TGR_INCLUDED
-#include "derivs2.F90"
+#include "TAT/TGRtensor/src/derivs2.F90"
 #undef TGR_INCLUDED
   subroutine calc_2position (shape, i,j, pos)
     integer, intent(in)  :: shape(2)
@@ -79,68 +82,109 @@ contains
     axy(pos) = f(1,2)
     ayy(pos) = f(2,2)
   end subroutine set_2tensor
-  subroutine get_2scalarderivs (a, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: a(*)
-    CCTK_REAL, intent(out) :: f(2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
-    call get_2derivs (a, f, pos, off, dx)
+  subroutine get_2scalarderivs (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
+    call get_2derivs (a, f, pos, off, dx, order)
   end subroutine get_2scalarderivs
-  subroutine get_2vectorderivs (ax,ay, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: ax(*),ay(*)
-    CCTK_REAL, intent(out) :: f(2,2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
+  subroutine get_2vectorderivs (ax,ay, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*)
+    CCTK_REAL, intent(out)           :: f(2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fx(2),fy(2)
-    call get_2derivs (ax, fx, pos, off, dx)
-    call get_2derivs (ay, fy, pos, off, dx)
+    call get_2derivs (ax, fx, pos, off, dx, order)
+    call get_2derivs (ay, fy, pos, off, dx, order)
     f(1,:) = fx
     f(2,:) = fy
   end subroutine get_2vectorderivs
-  subroutine get_2tensorderivs (axx,axy,ayy, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: axx(*),axy(*),ayy(*)
-    CCTK_REAL, intent(out) :: f(2,2,2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
+  subroutine get_2tensorderivs (axx,axy,ayy, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fxx(2),fxy(2),fyy(2)
-    call get_2derivs (axx, fxx, pos, off, dx)
-    call get_2derivs (axy, fxy, pos, off, dx)
-    call get_2derivs (ayy, fyy, pos, off, dx)
+    call get_2derivs (axx, fxx, pos, off, dx, order)
+    call get_2derivs (axy, fxy, pos, off, dx, order)
+    call get_2derivs (ayy, fyy, pos, off, dx, order)
     f(1,1,:) = fxx
     f(1,2,:) = fxy
     f(2,1,:) = fxy
     f(2,2,:) = fyy
   end subroutine get_2tensorderivs
-  subroutine get_2scalarderivs2 (a, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: a(*)
-    CCTK_REAL, intent(out) :: f(2,2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
-    call get_2derivs2 (a, f, pos, off, dx)
+  subroutine get_2scalarderivs2 (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
+    call get_2derivs2 (a, f, pos, off, dx, order)
   end subroutine get_2scalarderivs2
-  subroutine get_2vectorderivs2 (ax,ay, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: ax(*),ay(*)
-    CCTK_REAL, intent(out) :: f(2,2,2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
+  subroutine get_2vectorderivs2 (ax,ay, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fx(2,2),fy(2,2)
-    call get_2derivs2 (ax, fx, pos, off, dx)
-    call get_2derivs2 (ay, fy, pos, off, dx)
+    call get_2derivs2 (ax, fx, pos, off, dx, order)
+    call get_2derivs2 (ay, fy, pos, off, dx, order)
     f(1,:,:) = fx
     f(2,:,:) = fy
   end subroutine get_2vectorderivs2
-  subroutine get_2tensorderivs2 (axx,axy,ayy, f, pos, off, dx)
-    CCTK_REAL, intent(in)  :: axx(*),axy(*),ayy(*)
-    CCTK_REAL, intent(out) :: f(2,2,2,2)
-    integer,   intent(in)  :: pos, off(2)
-    CCTK_REAL, intent(in)  :: dx(2)
+  subroutine get_2tensorderivs2 (axx,axy,ayy, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
     CCTK_REAL :: fxx(2,2),fxy(2,2),fyy(2,2)
-    call get_2derivs2 (axx, fxx, pos, off, dx)
-    call get_2derivs2 (axy, fxy, pos, off, dx)
-    call get_2derivs2 (ayy, fyy, pos, off, dx)
+    call get_2derivs2 (axx, fxx, pos, off, dx, order)
+    call get_2derivs2 (axy, fxy, pos, off, dx, order)
+    call get_2derivs2 (ayy, fyy, pos, off, dx, order)
     f(1,1,:,:) = fxx
     f(1,2,:,:) = fxy
     f(2,1,:,:) = fxy
     f(2,2,:,:) = fyy
   end subroutine get_2tensorderivs2
+  subroutine get_2scalarderivs3 (a, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: a(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
+    call get_2derivs3 (a, f, pos, off, dx, order)
+  end subroutine get_2scalarderivs3
+  subroutine get_2vectorderivs3 (ax,ay, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: ax(*),ay(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
+    CCTK_REAL :: fx(2,2,2),fy(2,2,2)
+    call get_2derivs3 (ax, fx, pos, off, dx, order)
+    call get_2derivs3 (ay, fy, pos, off, dx, order)
+    f(1,:,:,:) = fx
+    f(2,:,:,:) = fy
+  end subroutine get_2vectorderivs3
+  subroutine get_2tensorderivs3 (axx,axy,ayy, f, pos, off, dx, order)
+    CCTK_REAL, intent(in)            :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out)           :: f(2,2,2,2,2)
+    integer,   intent(in)            :: pos, off(2)
+    CCTK_REAL, intent(in)            :: dx(2)
+    integer,   intent(in),  optional :: order
+    CCTK_REAL :: fxx(2,2,2),fxy(2,2,2),fyy(2,2,2)
+    call get_2derivs3 (axx, fxx, pos, off, dx, order)
+    call get_2derivs3 (axy, fxy, pos, off, dx, order)
+    call get_2derivs3 (ayy, fyy, pos, off, dx, order)
+    f(1,1,:,:,:) = fxx
+    f(1,2,:,:,:) = fxy
+    f(2,1,:,:,:) = fxy
+    f(2,2,:,:,:) = fyy
+  end subroutine get_2tensorderivs3
 end module pointwise2