lib/Math/Trig/Radial.pm

   1 package Math::Trig::Radial;
   2
   3 use strict;
   4 use vars qw(@ISA @EXPORT);
   5 @ISA = qw(Exporter);
   6
   7 @EXPORT =
   8     qw(
   9        cartesian_to_cylindrical
  10        cartesian_to_spherical
  11        cylindrical_to_cartesian
  12        cylindrical_to_spherical
  13        spherical_to_cartesian
  14        spherical_to_cylindrical
  15        great_circle_distance
  16       );
  17
  18 use Math::Trig;
  19
  20 sub pip2 { pi/2 }
  21
  22 =pod
  23
  24 =head1 NAME
  25
  26 Math::Trig::Radial - spherical and cylindrical trigonometry
  27
  28 =head1 SYNOPSIS
  29
  30     use Math::Trig::Radial;
  31
  32     ($rho, $theta, $z)     = cartesian_to_cylindrical($x, $y, $z);
  33     ($rho, $theta, $phi)   = cartesian_to_spherical($x, $y, $z);
  34     ($x, $y, $z)           = cylindrical_to_cartesian($rho, $theta, $z);
  35     ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
  36     ($x, $y, $z)           = spherical_to_cartesian($rho, $theta, $phi);
  37     ($rho_c, $theta, $z)   = spherical_to_cylindrical($rho_s, $theta, $phi);
  38
  39 =head1 DESCRIPTION
  40
  41 This module contains a few basic spherical and cylindrical
  42 trigonometric formulas.  B<All angles are in radians>, if needed
  43 use C<Math::Trig> angle unit conversions.
  44
  45 =head2 COORDINATE SYSTEMS
  46
  47 B<Cartesian> coordinates are the usual rectangular I<xyz>-coordinates.
  48
  49 Spherical coordinates are three-dimensional coordinates which define a
  50 point in three-dimensional space.  They are based on a sphere surface.
  51 The radius of the sphere is B<rho>, also known as the I<radial>
  52 coordinate.  The angle in the I<xy>-plane (around the I<z>-axis) is
  53 B<theta>, also known as the I<azimuthal> coordinate.  The angle from
  54 the I<z>-axis is B<phi>, also known as the I<polar> coordinate.  The
  55 `North Pole' is therefore I<0, 0, rho>, and the `Bay of Guinea' (think
  56 Africa) I<0, pi/2, rho>.
  57
  58 Cylindrical coordinates are three-dimensional coordinates which define
  59 a point in three-dimensional space.  They are based on a cylinder
  60 surface.  The radius of the cylinder is B<rho>, also known as the
  61 I<radial> coordinate.  The angle in the I<xy>-plane (around the
  62 I<z>-axis) is B<theta>, also known as the I<azimuthal> coordinate.
  63 The third coordinate is the I<z>.
  64
  65 =head2 CONVERSIONS
  66
  67 Conversions to and from spherical and cylindrical coordinates are
  68 available.  Please notice that the conversions are not necessarily
  69 reversible because of the equalities like I<pi> angles equals I<-pi>
  70 angles.
  71
  72 =over 4
  73
  74 =item cartesian_to_cylindrical
  75
  76         ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
  77
  78 =item cartesian_to_spherical
  79
  80         ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
  81
  82 =item cylindrical_to_cartesian
  83
  84         ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
  85
  86 =item cylindrical_to_spherical
  87
  88         ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
  89
  90 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
  91
  92 =item spherical_to_cartesian
  93
  94         ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
  95
  96 =item spherical_to_cylindrical
  97
  98         ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
  99
 100 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
 101
 102 =back
 103
 104 =head2 GREAT CIRCLE DISTANCE
 105
 106     $distance = great_circle_distance($x0, $y0, $z0, $x1, $y1, $z1 [, $rho]);
 107
 108 The I<great circle distance> is the shortest distance between two
 109 points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
 110 optional, it defaults to 1 (the unit sphere), therefore the distance
 111 defaults to radians.  The coordinates C<$x0> ... C<$z1> are in
 112 cartesian coordinates.
 113
 114 =head EXAMPLES
 115
 116 To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N
 117 139.8E) in kilometers:
 118
 119         use Math::Trig::Radial;
 120         use Math::Trig;
 121
 122         my @L = spherical_to_cartesian(1, map { deg2rad $_ } qw(51.3  -0.5));
 123         my @T = spherical_to_cartesian(1, map { deg2rad $_ } qw(35.7 139.8));
 124
 125         $km = great_circle_distance(@L, @T, 6378);
 126
 127 The answer may be off by up to 0.3% because of the irregular (slightly
 128 aspherical) form of the Earth.
 129
 130 =head2 AUTHOR
 131
 132 Jarkko Hietaniemi F<E<lt>jhi@iki.fiE<gt>>
 133
 134 =cut
 135
 136 sub cartesian_to_spherical {
 137     my ( $x, $y, $z ) = @_;
 138
 139     my $rho = sqrt( $x * $x + $y * $y + $z * $z );
 140
 141     return ( $rho,
 142              atan2( $y, $x ),
 143              $rho ? acos( $z / $rho ) : 0 );
 144 }
 145
 146 sub spherical_to_cartesian {
 147     my ( $rho, $theta, $phi ) = @_;
 148
 149     return ( $rho * cos( $theta ) * sin( $phi ),
 150              $rho * sin( $theta ) * sin( $phi ),
 151              $rho * cos( $phi   ) );
 152 }
 153
 154 sub spherical_to_cylindrical {
 155     my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
 156
 157     return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
 158 }
 159
 160 sub cartesian_to_cylindrical {
 161     my ( $x, $y, $z ) = @_;
 162
 163     return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
 164 }
 165
 166 sub cylindrical_to_cartesian {
 167     my ( $rho, $theta, $z ) = @_;
 168
 169     return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
 170 }
 171
 172 sub cylindrical_to_spherical {
 173     return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
 174 }
 175
 176 sub great_circle_distance {
 177     my ( $x0, $y0, $z0, $x1, $y1, $z1, $rho ) = @_;
 178
 179     $rho = 1 unless defined $rho; # Default to the unit sphere.
 180
 181     my ( $r0, $theta0, $phi0 ) = cartesian_to_spherical( $x0, $y0, $z0 );
 182     my ( $r1, $theta1, $phi1 ) = cartesian_to_spherical( $x1, $y1, $z1 );
 183
 184     my $lat0 = pip2 - $phi0;
 185     my $lat1 = pip2 - $phi1;
 186
 187     return $rho *
 188         acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
 189              sin( $lat0 ) * sin( $lat1 ) );
 190 }
 191
 192 1;
 193