File:/home/sbrandt/cactus/Cactus/configs/sim2/build/Extract/cartesian_to_spherical.f
1:# 4 "/home/sbrandt/cactus/Cactus/arrangements/EinsteinAnalysis/Extract/src/cartesian_to_spherical.F"
2:c     ==================================================================
3:
4:      SUBROUTINE cartesian_to_spherical(theta,phi,r,gxxs,gxys,gxzs,
5:     &     gyys,gyzs, gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys,
6:     &     dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps)
7:
8:c     ------------------------------------------------------------------
9:c
10:c     Convert components of a symmetric 2nd rank 3D tensor on a 2-sphere 
11:c     from Cartesian to spherical polar coordinates
12:c
13:c                (x,y,z)        ---->         (r,t,p)
14:c
15:c           ( gxx  gxy gxz )             ( grr  grt grp )
16:c           (  .   gyy gyz )    ---->    (  .   gpp gtp )   
17:c           (  .    .  gzz )             (  .    .  gpp ) 
18:c
19:c     Also convert radial derivatives of the tensor components
20:c
21:c     ------------------------------------------------------------------
22:
23:      IMPLICIT NONE
24:
25:c     Input variables
26:      REAL*8,INTENT(IN) ::
27:     &    r
28:      REAL*8,INTENT(IN),DIMENSION (:) ::
29:     &    theta,phi
30:      REAL*8,INTENT(IN),DIMENSION (:,:) ::
31:     &    gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,dgyzs,
32:     &    dgzzs
33:
34:c     Output variables
35:      REAL*8,INTENT(OUT),DIMENSION (:,:) ::
36:     &    grrs,grts,grps,gtts,gtps,
37:     &    gpps,dgtts,dgtps,dgpps
38:
39:c     Local variables
40:      INTEGER ::
41:     &    it,ip,Nt,Np
42:      REAL*8 ::
43:     &    ct,st,cp,sp,ct2,st2,cp2,sp2,r2,gxx,gxy,gxz,gyy,gyz,gzz,dgxx,
44:     &    dgxy,dgxz,dgyy,dgyz,dgzz
45:      REAL*8,PARAMETER ::
46:     &    two = 2.0D0
47:
48:c     ------------------------------------------------------------------
49:
50:      Nt = SIZE(theta)
51:      Np = SIZE(phi)
52:
53:      r2 = r*r
54:
55:      DO it = 1, Nt
56:
57:        ct = COS(theta(it)) ; ct2 = ct*ct
58:        st = SIN(theta(it)) ; st2 = st*st
59:
60:          DO ip = 1, Np
61:
62:          cp = COS(phi(ip)) ; cp2 = cp*cp
63:          sp = SIN(phi(ip)) ; sp2 = sp*sp
64:
65:          gxx = gxxs(it,ip) ; dgxx = dgxxs(it,ip)
66:          gxy = gxys(it,ip) ; dgxy = dgxys(it,ip)
67:          gxz = gxzs(it,ip) ; dgxz = dgxzs(it,ip)
68:          gyy = gyys(it,ip) ; dgyy = dgyys(it,ip)
69:          gyz = gyzs(it,ip) ; dgyz = dgyzs(it,ip)
70:          gzz = gzzs(it,ip) ; dgzz = dgzzs(it,ip)
71:
72:          grrs(it,ip) = st2*cp2*gxx+st2*sp2*gyy+ct2*gzz
73:     &       +two*st2*cp*sp*gxy+two*st*cp*ct*gxz+two*st*ct*sp*gyz
74:
75:          grts(it,ip) = r*(st*cp2*ct*gxx+two*st*ct*sp*cp*gxy
76:     &       +cp*(ct2-st2)*gxz+st*sp2*ct*gyy+sp*(ct2-st2)*gyz-ct*st*gzz)
77:# 78
78:     &
79:
80:          grps(it,ip) = r*st*(-st*sp*cp*gxx-st*(sp2-cp2)*gxy-sp*ct*gxz
81:     &       +st*sp*cp*gyy+ct*cp*gyz)
82:
83:          gtts(it,ip) = r2*(ct2*cp2*gxx+two*ct2*sp*cp*gxy
84:     &       -two*st*ct*cp*gxz+ct2*sp2*gyy-two*st*sp*ct*gyz+st2*gzz)
85:
86:          gtps(it,ip) = r2*st*(-cp*sp*ct*gxx-ct*(sp2-cp2)*gxy
87:     &       +st*sp*gxz+cp*sp*ct*gyy-st*cp*gyz)
88:
89:          gpps(it,ip) = r2*st2*(sp2*gxx-two*cp*sp*gxy+cp2*gyy)
90:
91:          dgtts(it,ip) = two/r*gtts(it,ip)+r2*(ct2*cp2*dgxx
92:     &       +two*ct2*sp*cp*dgxy-two*st*ct*cp*dgxz+ct2*sp2*dgyy
93:     &       -two*st*sp*ct*dgyz+st2*dgzz)
94:
95:          dgtps(it,ip) = two/r*gtps(it,ip)+r2*st*(-cp*sp*ct*dgxx
96:     &       -ct*(sp2-cp2)*dgxy+st*sp*dgxz+cp*sp*ct*dgyy-st*cp*dgyz)
97:
98:          dgpps(it,ip) = two/r*gpps(it,ip) +r2*st2*(sp2*dgxx
99:     &       -two*cp*sp*dgxy+cp2*dgyy)
100:
101:        END DO
102:
103:      END DO
104:
105:
106:      END SUBROUTINE cartesian_to_spherical