1 files changed, 116 insertions, 0 deletions
diff --git a/src/cartesian_to_spherical.F b/src/cartesian_to_spherical.F
new file mode 100644
index 0000000..d55f84f
--- /dev/null
+++ b/src/cartesian_to_spherical.F
@@ -0,0 +1,116 @@
+
+#include "cctk.h"
+
+c     ==================================================================
+
+      SUBROUTINE cartesian_to_spherical(theta,phi,r,gxxs,gxys,gxzs,
+     &     gyys,gyzs, gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys,
+     &     dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps)
+
+c     ------------------------------------------------------------------
+c
+c     Convert components of a symmetric 2nd rank 3D tensor on a 2-sphere 
+c     from Cartesian to spherical polar coordinates
+c
+c                (x,y,z)        ---->         (r,t,p)
+c
+c           ( gxx  gxy gxz )             ( grr  grt grp )
+c           (  .   gyy gyz )    ---->    (  .   gpp gtp )   
+c           (  .    .  gzz )             (  .    .  gpp ) 
+c
+c     Also convert radial derivatives of the tensor components
+c
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      CCTK_REAL,INTENT(IN) ::  
+     &    r
+      CCTK_REAL,INTENT(IN),DIMENSION (:) ::  
+     &    theta,phi
+      CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: 
+     &    gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,dgyzs,
+     &    dgzzs
+
+c     Output variables
+      CCTK_REAL,INTENT(OUT),DIMENSION (:,:) :: 
+     &    grrs,grts,grps,gtts,gtps,
+     &    gpps,dgtts,dgtps,dgpps
+
+c     Local variables
+      INTEGER ::
+     &    it,ip,Nt,Np
+      CCTK_REAL :: 
+     &    ct,st,cp,sp,ct2,st2,cp2,sp2,r2,gxx,gxy,gxz,gyy,gyz,gzz,dgxx,
+     &    dgxy,dgxz,dgyy,dgyz,dgzz
+      CCTK_REAL,PARAMETER ::
+     &    two = 2.0D0
+
+c     ------------------------------------------------------------------
+
+      Nt = SIZE(theta)
+      Np = SIZE(phi)
+
+      r2 = r*r
+
+      DO it = 1, Nt
+ 
+	ct = COS(theta(it)) ; ct2 = ct*ct
+        st = SIN(theta(it)) ; st2 = st*st
+
+          DO ip = 1, Np
+
+          cp = COS(phi(ip)) ; cp2 = cp*cp
+          sp = SIN(phi(ip)) ; sp2 = sp*sp
+
+	  gxx = gxxs(it,ip) ; dgxx = dgxxs(it,ip)
+	  gxy = gxys(it,ip) ; dgxy = dgxys(it,ip)
+	  gxz = gxzs(it,ip) ; dgxz = dgxzs(it,ip)
+          gyy = gyys(it,ip) ; dgyy = dgyys(it,ip)
+          gyz = gyzs(it,ip) ; dgyz = dgyzs(it,ip)
+          gzz = gzzs(it,ip) ; dgzz = dgzzs(it,ip)
+
+          grrs(it,ip) = st2*cp2*gxx+st2*sp2*gyy+ct2*gzz
+     &       +two*st2*cp*sp*gxy+two*st*cp*ct*gxz+two*st*ct*sp*gyz
+
+          grts(it,ip) = r*(st*cp2*ct*gxx+two*st*ct*sp*cp*gxy
+     &       +cp*(ct2-st2)*gxz+st*sp2*ct*gyy+sp*(ct2-st2)*gyz-ct*st*gzz)
+
+          grps(it,ip) = r*st*(-st*sp*cp*gxx-st*(sp2-cp2)*gxy-sp*ct*gxz
+     &       +st*sp*cp*gyy+ct*cp*gyz)
+     
+          gtts(it,ip) = r2*(ct2*cp2*gxx+two*ct2*sp*cp*gxy
+     &       -two*st*ct*cp*gxz+ct2*sp2*gyy-two*st*sp*ct*gyz+st2*gzz)
+     
+          gtps(it,ip) = r2*st*(-cp*sp*ct*gxx-ct*(sp2-cp2)*gxy
+     &       +st*sp*gxz+cp*sp*ct*gyy-st*cp*gyz)
+     
+          gpps(it,ip) = r2*st2*(sp2*gxx-two*cp*sp*gxy+cp2*gyy)
+     
+          dgtts(it,ip) = two/r*gtts(it,ip)+r2*(ct2*cp2*dgxx
+     &       +two*ct2*sp*cp*dgxy-two*st*ct*cp*dgxz+ct2*sp2*dgyy
+     &       -two*st*sp*ct*dgyz+st2*dgzz)
+
+          dgtps(it,ip) = two/r*gtps(it,ip)+r2*st*(-cp*sp*ct*dgxx
+     &       -ct*(sp2-cp2)*dgxy+st*sp*dgxz+cp*sp*ct*dgyy-st*cp*dgyz)
+     
+          dgpps(it,ip) = two/r*gpps(it,ip) +r2*st2*(sp2*dgxx
+     &       -two*cp*sp*dgxy+cp2*dgyy)
+
+        END DO
+	  				 
+      END DO
+
+
+      END SUBROUTINE cartesian_to_spherical
+
+
+
+
+
+
+
+
+
+