diff options
Diffstat (limited to 'src/conversion.F90')
| -rw-r--r-- | src/conversion.F90 | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/src/conversion.F90 b/src/conversion.F90 new file mode 100644 index 0000000..d6e30e9 --- /dev/null +++ b/src/conversion.F90 @@ -0,0 +1,185 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module conversion + implicit none + DECLARE_CCTK_PARAMETERS + private + + public make_cylindrical2cartesian + public make_spherical2cartesian + public make_cartesian2spherical + +contains + + subroutine make_cylindrical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho + + ! cylindrical coordinates: rho phi z + ! Cartesian coordinates: x y z + ! rho^2 = x^2 + y^2 + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[cyl](p)_j + ! where p = p[cart] + + ! Transformation: + ! x = rho cos(phi) + ! y = rho sin(phi) + ! z = z + ! rho = sqrt(x^2+y^2) + ! phi = atan(y/x) + ! z = z + + ! Jacobian: + ! drho/dx = x/rho + ! drho/dy = y/rho + ! drho/dz = 0 + ! dphi/dx = -y/rho^2 + ! dphi/dy = x/rho^2 + ! dphi/dz = 0 + ! dz /dx = 0 + ! dz /dy = 0 + ! dz /dz = 1 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + + tt(1,1) = x/(rho+eps) + tt(2,1) = y/(rho+eps) + tt(3,1) = 0 + tt(1,2) = -y/(rho**2+eps) + tt(2,2) = x/(rho**2+eps) + tt(3,2) = 0 + tt(1,3) = 0 + tt(2,3) = 0 + tt(3,3) = 1 + + end subroutine make_cylindrical2cartesian + + subroutine make_spherical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[spher](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dr /dx = x/r + ! dr /dy = y/r + ! dr /dz = z/r + ! dtheta/dx = xz/r^2rho ! 0 + ! dtheta/dy = yz/r^2rho ! 0 + ! dtheta/dz = -rho/r^2 ! 1/rho + ! dphi /dx = -y/rho^2 + ! dphi /dy = x/rho^2 + ! dphi /dz = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = y/(r+eps) + tt(3,1) = z/(r+eps) + tt(1,2) = x*z/(r**2*rho+eps) + tt(2,2) = y*z/(r**2*rho+eps) + tt(3,2) = -rho**2/(r**2*rho+eps) + tt(1,3) = -y/(rho**2+eps) + tt(2,3) = x/(rho**2+eps) + tt(3,3) = 0 + + end subroutine make_spherical2cartesian + + subroutine make_cartesian2spherical (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[spher](p)_i = TT(p)_i^j A[cart](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dx/dr = sin(theta) cos(phi) = x/r + ! dx/dtheta = r cos(theta) cos(phi) = x z/rho + ! dx/dphi = - r sin(theta) sin(phi) = -y + ! dy/dr = sin(theta) sin(phi) = y/r + ! dy/dtheta = r cos(theta) sin(phi) = y z/rho + ! dy/dphi = r sin(theta) cos(phi) = x + ! dz/dr = cos(theta) = z/r + ! dz/dtheta = - r sin(theta) = -rho + ! dz/dphi = 0 = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = x * z/(rho+eps) + tt(3,1) = -y + tt(1,2) = y/(r+eps) + tt(2,2) = y * z/(rho+eps) + tt(3,2) = x + tt(1,3) = z/(r+eps) + tt(2,3) = -rho + tt(3,3) = 0 + + end subroutine make_cartesian2spherical + +end module conversion |