[1765] | 1 | module spherical |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | contains |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | subroutine rectsph(rect,col,azm,r,r2) |
---|
| 9 | |
---|
| 10 | ! Converts rectangular coordinates to spherical coordinates. |
---|
| 11 | ! (Ox) is taken as the polar axis. |
---|
| 12 | ! Azimuth is the angle in the (Oyz) plane, measured from vector y. |
---|
| 13 | ! If neither "r" nor "r2" is present then they are assumed to be |
---|
| 14 | ! equal to 1. |
---|
| 15 | |
---|
| 16 | real, intent(in):: rect(3) |
---|
| 17 | ! (Rectangular coordinates x,y,z. Should be different from (0,0 |
---|
| 18 | ! ,0)) |
---|
| 19 | |
---|
| 20 | real, intent(out):: col ! Colatitude (in radians) |
---|
| 21 | real, intent(out):: azm ! Azimuth (in radians). -pi < azm <= pi |
---|
| 22 | real, intent(out), optional:: r ! Radius |
---|
| 23 | real, intent(out), optional:: r2 ! Square of radius |
---|
| 24 | |
---|
| 25 | ! Local variables: |
---|
| 26 | logical p, p2 |
---|
| 27 | real d ! Radius |
---|
| 28 | real d2 ! Square of radius |
---|
| 29 | |
---|
| 30 | !---------------------------------------------------------------- |
---|
| 31 | |
---|
| 32 | p = present(r) |
---|
| 33 | p2 = present(r2) |
---|
| 34 | |
---|
| 35 | if (.not. p .and. .not. p2) then |
---|
| 36 | ! Distance to the origin is assumed to be 1 |
---|
| 37 | col = acos(rect(1)) |
---|
| 38 | else |
---|
| 39 | ! Compute the distance to the origin: |
---|
| 40 | d2 = dot_product(rect,rect) |
---|
| 41 | d = sqrt(d2) |
---|
| 42 | col = acos(rect(1)/d) |
---|
| 43 | if (p) r = d |
---|
| 44 | if (p2) r2 = d2 |
---|
| 45 | end if |
---|
| 46 | |
---|
| 47 | if (rect(2) == 0. .and. rect(3) == 0.) then |
---|
| 48 | azm = 0. |
---|
| 49 | ! (arbitrary value, azimuth is not well defined since vector |
---|
| 50 | ! position is parallel |
---|
| 51 | ! to polar axis) |
---|
| 52 | else |
---|
| 53 | azm = atan2(rect(3),rect(2)) |
---|
| 54 | end if |
---|
| 55 | |
---|
| 56 | end subroutine rectsph |
---|
| 57 | |
---|
| 58 | !******************************************************************* |
---|
| 59 | |
---|
| 60 | real function sphbase(col,azm) |
---|
| 61 | |
---|
| 62 | ! This function returns the matrix of the spherical base: (radial |
---|
| 63 | ! vector, colatitude vector, azimuthal vector) in the cartesian |
---|
| 64 | ! vector base: (x, y, z). Vector x is assumed to be the polar |
---|
| 65 | ! vector. Colatitude "col" and azimuth "azm" are the angular |
---|
| 66 | ! spherical coordinates of the radial vector (azimuth is measured |
---|
| 67 | ! from vector "y"). Note that if col = 0 (radial vector parallel |
---|
| 68 | ! to x) then the choice of either the colatitude vector or the |
---|
| 69 | ! azimuthal vector is arbitrary. The choice made depends on the |
---|
| 70 | ! input value of azm. If azm is 0 then sphbase will return vector |
---|
| 71 | ! y as the colatitude vector ("sphbase" is then the identity |
---|
| 72 | ! matrix). Uses "sphrect", which converts spherical coordinates to |
---|
| 73 | ! rectangular coordinates. |
---|
| 74 | |
---|
| 75 | dimension sphbase(3,3) ! Transformation matrix (no dimension) |
---|
| 76 | |
---|
| 77 | real, intent(in):: col ! Colatitude (in radians) |
---|
| 78 | real, intent(in):: azm ! Azimuth (in radians) |
---|
| 79 | |
---|
| 80 | ! Local variable: |
---|
| 81 | real pi |
---|
| 82 | |
---|
| 83 | !----------------------------------------------------------------- |
---|
| 84 | |
---|
| 85 | pi = acos(-1.) |
---|
| 86 | |
---|
| 87 | ! Radial vector: |
---|
| 88 | sphbase(:,1) = sphrect(col, azm) |
---|
| 89 | |
---|
| 90 | ! Colatitude vector: |
---|
| 91 | sphbase(:,2) = sphrect(col + pi/2, azm) |
---|
| 92 | |
---|
| 93 | ! Azimuthal vector: |
---|
| 94 | sphbase(:,3) = (/0., - sin(azm), cos(azm)/) |
---|
| 95 | |
---|
| 96 | end function sphbase |
---|
| 97 | |
---|
| 98 | !***************************************************************** |
---|
| 99 | |
---|
| 100 | real function sphrect(col,azm,r) |
---|
| 101 | |
---|
| 102 | ! Converts spherical coordinates to rectangular coordinates. (Ox) |
---|
| 103 | ! is taken as the polar axis. Azimuth "azm" is the angle in the |
---|
| 104 | ! (Oyz) plane, measured from vector y. |
---|
| 105 | ! If r is not present then we assume r = 1. |
---|
| 106 | |
---|
| 107 | dimension sphrect(3) ! Rectangular coordinates = (x,y,z) |
---|
| 108 | |
---|
| 109 | real, intent(in):: col ! Colatitude (in radians) |
---|
| 110 | real, intent(in):: azm ! Azimuth (in radians) |
---|
| 111 | real, intent(in), optional:: r ! Radius |
---|
| 112 | |
---|
| 113 | !--------------------------------------------------------------- |
---|
| 114 | |
---|
| 115 | sphrect = (/cos(col), sin(col) * cos(azm), sin(col) * sin(azm)/) |
---|
| 116 | if (present(r)) sphrect = sphrect * r |
---|
| 117 | |
---|
| 118 | end function sphrect |
---|
| 119 | |
---|
| 120 | end module spherical |
---|