Commit | Line | Data |
---|---|---|

5aabfad6 PP |
1 | # |

2 | # Trigonometric functions, mostly inherited from Math::Complex. | |

d54bf66f | 3 | # -- Jarkko Hietaniemi, since April 1997 |

5cd24f17 | 4 | # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex) |

5aabfad6 PP |
5 | # |

6 | ||

7 | require Exporter; | |

8 | package Math::Trig; | |

9 | ||

17f410f9 | 10 | use 5.005_64; |

5aabfad6 PP |
11 | use strict; |

12 | ||

13 | use Math::Complex qw(:trig); | |

14 | ||

17f410f9 | 15 | our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS); |

5aabfad6 PP |
16 | |

17 | @ISA = qw(Exporter); | |

18 | ||

19 | $VERSION = 1.00; | |

20 | ||

ace5de91 GS |
21 | my @angcnv = qw(rad2deg rad2grad |

22 | deg2rad deg2grad | |

23 | grad2rad grad2deg); | |

5aabfad6 PP |
24 | |

25 | @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, | |

26 | @angcnv); | |

27 | ||

d54bf66f JH |
28 | my @rdlcnv = qw(cartesian_to_cylindrical |

29 | cartesian_to_spherical | |

30 | cylindrical_to_cartesian | |

31 | cylindrical_to_spherical | |

32 | spherical_to_cartesian | |

33 | spherical_to_cylindrical); | |

34 | ||

35 | @EXPORT_OK = (@rdlcnv, 'great_circle_distance'); | |

36 | ||

37 | %EXPORT_TAGS = ('radial' => [ @rdlcnv ]); | |

38 | ||

6570f784 GS |
39 | sub pi2 () { 2 * pi } # use constant generates warning |

40 | sub pip2 () { pi / 2 } # use constant generates warning | |

d54bf66f JH |
41 | use constant DR => pi2/360; |

42 | use constant RD => 360/pi2; | |

43 | use constant DG => 400/360; | |

44 | use constant GD => 360/400; | |

45 | use constant RG => 400/pi2; | |

46 | use constant GR => pi2/400; | |

5aabfad6 PP |
47 | |

48 | # | |

49 | # Truncating remainder. | |

50 | # | |

51 | ||

52 | sub remt ($$) { | |

53 | # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available. | |

54 | $_[0] - $_[1] * int($_[0] / $_[1]); | |

55 | } | |

56 | ||

57 | # | |

58 | # Angle conversions. | |

59 | # | |

60 | ||

ace5de91 | 61 | sub rad2deg ($) { remt(RD * $_[0], 360) } |

5aabfad6 | 62 | |

ace5de91 | 63 | sub deg2rad ($) { remt(DR * $_[0], pi2) } |

5aabfad6 | 64 | |

ace5de91 | 65 | sub grad2deg ($) { remt(GD * $_[0], 360) } |

5aabfad6 | 66 | |

ace5de91 | 67 | sub deg2grad ($) { remt(DG * $_[0], 400) } |

5aabfad6 | 68 | |

ace5de91 | 69 | sub rad2grad ($) { remt(RG * $_[0], 400) } |

5aabfad6 | 70 | |

ace5de91 | 71 | sub grad2rad ($) { remt(GR * $_[0], pi2) } |

5aabfad6 | 72 | |

d54bf66f JH |
73 | sub cartesian_to_spherical { |

74 | my ( $x, $y, $z ) = @_; | |

75 | ||

76 | my $rho = sqrt( $x * $x + $y * $y + $z * $z ); | |

77 | ||

78 | return ( $rho, | |

79 | atan2( $y, $x ), | |

80 | $rho ? acos( $z / $rho ) : 0 ); | |

81 | } | |

82 | ||

83 | sub spherical_to_cartesian { | |

84 | my ( $rho, $theta, $phi ) = @_; | |

85 | ||

86 | return ( $rho * cos( $theta ) * sin( $phi ), | |

87 | $rho * sin( $theta ) * sin( $phi ), | |

88 | $rho * cos( $phi ) ); | |

89 | } | |

90 | ||

91 | sub spherical_to_cylindrical { | |

92 | my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); | |

93 | ||

94 | return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); | |

95 | } | |

96 | ||

97 | sub cartesian_to_cylindrical { | |

98 | my ( $x, $y, $z ) = @_; | |

99 | ||

100 | return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); | |

101 | } | |

102 | ||

103 | sub cylindrical_to_cartesian { | |

104 | my ( $rho, $theta, $z ) = @_; | |

105 | ||

106 | return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); | |

107 | } | |

108 | ||

109 | sub cylindrical_to_spherical { | |

110 | return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); | |

111 | } | |

112 | ||

113 | sub great_circle_distance { | |

114 | my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_; | |

115 | ||

116 | $rho = 1 unless defined $rho; # Default to the unit sphere. | |

117 | ||

118 | my $lat0 = pip2 - $phi0; | |

119 | my $lat1 = pip2 - $phi1; | |

120 | ||

121 | return $rho * | |

122 | acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + | |

123 | sin( $lat0 ) * sin( $lat1 ) ); | |

124 | } | |

125 | ||

126 | =pod | |

127 | ||

5aabfad6 PP |
128 | =head1 NAME |

129 | ||

130 | Math::Trig - trigonometric functions | |

131 | ||

132 | =head1 SYNOPSIS | |

133 | ||

134 | use Math::Trig; | |

3cb6de81 | 135 | |

5aabfad6 PP |
136 | $x = tan(0.9); |

137 | $y = acos(3.7); | |

138 | $z = asin(2.4); | |

3cb6de81 | 139 | |

5aabfad6 PP |
140 | $halfpi = pi/2; |

141 | ||

ace5de91 | 142 | $rad = deg2rad(120); |

5aabfad6 PP |
143 | |

144 | =head1 DESCRIPTION | |

145 | ||

146 | C<Math::Trig> defines many trigonometric functions not defined by the | |

4ae80833 | 147 | core Perl which defines only the C<sin()> and C<cos()>. The constant |

5aabfad6 PP |
148 | B<pi> is also defined as are a few convenience functions for angle |

149 | conversions. | |

150 | ||

151 | =head1 TRIGONOMETRIC FUNCTIONS | |

152 | ||

153 | The tangent | |

154 | ||

d54bf66f JH |
155 | =over 4 |

156 | ||

157 | =item B<tan> | |

158 | ||

159 | =back | |

5aabfad6 PP |
160 | |

161 | The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot | |

162 | are aliases) | |

163 | ||

d54bf66f | 164 | B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan> |

5aabfad6 PP |
165 | |

166 | The arcus (also known as the inverse) functions of the sine, cosine, | |

167 | and tangent | |

168 | ||

d54bf66f | 169 | B<asin>, B<acos>, B<atan> |

5aabfad6 PP |
170 | |

171 | The principal value of the arc tangent of y/x | |

172 | ||

d54bf66f | 173 | B<atan2>(y, x) |

5aabfad6 PP |
174 | |

175 | The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc | |

176 | and acotan/acot are aliases) | |

177 | ||

d54bf66f | 178 | B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan> |

5aabfad6 PP |
179 | |

180 | The hyperbolic sine, cosine, and tangent | |

181 | ||

d54bf66f | 182 | B<sinh>, B<cosh>, B<tanh> |

5aabfad6 PP |
183 | |

184 | The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch | |

185 | and cotanh/coth are aliases) | |

186 | ||

d54bf66f | 187 | B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh> |

5aabfad6 PP |
188 | |

189 | The arcus (also known as the inverse) functions of the hyperbolic | |

190 | sine, cosine, and tangent | |

191 | ||

d54bf66f | 192 | B<asinh>, B<acosh>, B<atanh> |

5aabfad6 PP |
193 | |

194 | The arcus cofunctions of the hyperbolic sine, cosine, and tangent | |

195 | (acsch/acosech and acoth/acotanh are aliases) | |

196 | ||

d54bf66f | 197 | B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh> |

5aabfad6 PP |
198 | |

199 | The trigonometric constant B<pi> is also defined. | |

200 | ||

d54bf66f | 201 | $pi2 = 2 * B<pi>; |

5aabfad6 | 202 | |

5cd24f17 PP |
203 | =head2 ERRORS DUE TO DIVISION BY ZERO |

204 | ||

205 | The following functions | |

206 | ||

d54bf66f | 207 | acoth |

5cd24f17 | 208 | acsc |

5cd24f17 | 209 | acsch |

d54bf66f JH |
210 | asec |

211 | asech | |

212 | atanh | |

213 | cot | |

214 | coth | |

215 | csc | |

216 | csch | |

217 | sec | |

218 | sech | |

219 | tan | |

220 | tanh | |

5cd24f17 PP |
221 | |

222 | cannot be computed for all arguments because that would mean dividing | |

8c03c583 JH |
223 | by zero or taking logarithm of zero. These situations cause fatal |

224 | runtime errors looking like this | |

5cd24f17 PP |
225 | |

226 | cot(0): Division by zero. | |

227 | (Because in the definition of cot(0), the divisor sin(0) is 0) | |

228 | Died at ... | |

229 | ||

8c03c583 JH |
230 | or |

231 | ||

232 | atanh(-1): Logarithm of zero. | |

233 | Died at... | |

234 | ||

235 | For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, | |

236 | C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the | |

237 | C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the | |

238 | C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the | |

239 | C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * | |

240 | pi>, where I<k> is any integer. | |

5cd24f17 PP |
241 | |

242 | =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS | |

5aabfad6 PP |
243 | |

244 | Please note that some of the trigonometric functions can break out | |

245 | from the B<real axis> into the B<complex plane>. For example | |

246 | C<asin(2)> has no definition for plain real numbers but it has | |

247 | definition for complex numbers. | |

248 | ||

249 | In Perl terms this means that supplying the usual Perl numbers (also | |

250 | known as scalars, please see L<perldata>) as input for the | |

251 | trigonometric functions might produce as output results that no more | |

252 | are simple real numbers: instead they are complex numbers. | |

253 | ||

254 | The C<Math::Trig> handles this by using the C<Math::Complex> package | |

255 | which knows how to handle complex numbers, please see L<Math::Complex> | |

256 | for more information. In practice you need not to worry about getting | |

257 | complex numbers as results because the C<Math::Complex> takes care of | |

258 | details like for example how to display complex numbers. For example: | |

259 | ||

260 | print asin(2), "\n"; | |

3cb6de81 | 261 | |

5aabfad6 PP |
262 | should produce something like this (take or leave few last decimals): |

263 | ||

264 | 1.5707963267949-1.31695789692482i | |

265 | ||

5cd24f17 PP |
266 | That is, a complex number with the real part of approximately C<1.571> |

267 | and the imaginary part of approximately C<-1.317>. | |

5aabfad6 | 268 | |

d54bf66f | 269 | =head1 PLANE ANGLE CONVERSIONS |

5aabfad6 PP |
270 | |

271 | (Plane, 2-dimensional) angles may be converted with the following functions. | |

272 | ||

ace5de91 GS |
273 | $radians = deg2rad($degrees); |

274 | $radians = grad2rad($gradians); | |

3cb6de81 | 275 | |

ace5de91 GS |
276 | $degrees = rad2deg($radians); |

277 | $degrees = grad2deg($gradians); | |

3cb6de81 | 278 | |

ace5de91 GS |
279 | $gradians = deg2grad($degrees); |

280 | $gradians = rad2grad($radians); | |

5aabfad6 | 281 | |

5cd24f17 | 282 | The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians. |

5aabfad6 | 283 | |

d54bf66f JH |
284 | =head1 RADIAL COORDINATE CONVERSIONS |

285 | ||

286 | B<Radial coordinate systems> are the B<spherical> and the B<cylindrical> | |

287 | systems, explained shortly in more detail. | |

288 | ||

289 | You can import radial coordinate conversion functions by using the | |

290 | C<:radial> tag: | |

291 | ||

292 | use Math::Trig ':radial'; | |

293 | ||

294 | ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); | |

295 | ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); | |

296 | ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); | |

297 | ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); | |

298 | ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); | |

299 | ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); | |

300 | ||

301 | B<All angles are in radians>. | |

302 | ||

303 | =head2 COORDINATE SYSTEMS | |

304 | ||

305 | B<Cartesian> coordinates are the usual rectangular I<(x, y, | |

306 | z)>-coordinates. | |

307 | ||

308 | Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional | |

309 | coordinates which define a point in three-dimensional space. They are | |

310 | based on a sphere surface. The radius of the sphere is B<rho>, also | |

311 | known as the I<radial> coordinate. The angle in the I<xy>-plane | |

312 | (around the I<z>-axis) is B<theta>, also known as the I<azimuthal> | |

313 | coordinate. The angle from the I<z>-axis is B<phi>, also known as the | |

314 | I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and | |

315 | the `Bay of Guinea' (think of the missing big chunk of Africa) I<0, | |

4b0d1da8 JH |
316 | pi/2, rho>. In geographical terms I<phi> is latitude (northward |

317 | positive, southward negative) and I<theta> is longitude (eastward | |

318 | positive, westward negative). | |

d54bf66f | 319 | |

4b0d1da8 | 320 | B<BEWARE>: some texts define I<theta> and I<phi> the other way round, |

d54bf66f JH |
321 | some texts define the I<phi> to start from the horizontal plane, some |

322 | texts use I<r> in place of I<rho>. | |

323 | ||

324 | Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional | |

325 | coordinates which define a point in three-dimensional space. They are | |

326 | based on a cylinder surface. The radius of the cylinder is B<rho>, | |

327 | also known as the I<radial> coordinate. The angle in the I<xy>-plane | |

328 | (around the I<z>-axis) is B<theta>, also known as the I<azimuthal> | |

329 | coordinate. The third coordinate is the I<z>, pointing up from the | |

330 | B<theta>-plane. | |

331 | ||

332 | =head2 3-D ANGLE CONVERSIONS | |

333 | ||

334 | Conversions to and from spherical and cylindrical coordinates are | |

335 | available. Please notice that the conversions are not necessarily | |

336 | reversible because of the equalities like I<pi> angles being equal to | |

337 | I<-pi> angles. | |

338 | ||

339 | =over 4 | |

340 | ||

341 | =item cartesian_to_cylindrical | |

342 | ||

343 | ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); | |

344 | ||

345 | =item cartesian_to_spherical | |

346 | ||

347 | ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); | |

348 | ||

349 | =item cylindrical_to_cartesian | |

350 | ||

351 | ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); | |

352 | ||

353 | =item cylindrical_to_spherical | |

354 | ||

355 | ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); | |

356 | ||

357 | Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>. | |

358 | ||

359 | =item spherical_to_cartesian | |

360 | ||

361 | ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); | |

362 | ||

363 | =item spherical_to_cylindrical | |

364 | ||

365 | ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); | |

366 | ||

367 | Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. | |

368 | ||

369 | =back | |

370 | ||

371 | =head1 GREAT CIRCLE DISTANCES | |

372 | ||

373 | You can compute spherical distances, called B<great circle distances>, | |

374 | by importing the C<great_circle_distance> function: | |

375 | ||

376 | use Math::Trig 'great_circle_distance' | |

377 | ||

4b0d1da8 | 378 | $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); |

d54bf66f JH |
379 | |

380 | The I<great circle distance> is the shortest distance between two | |

381 | points on a sphere. The distance is in C<$rho> units. The C<$rho> is | |

382 | optional, it defaults to 1 (the unit sphere), therefore the distance | |

383 | defaults to radians. | |

384 | ||

4b0d1da8 JH |
385 | If you think geographically the I<theta> are longitudes: zero at the |

386 | Greenwhich meridian, eastward positive, westward negative--and the | |

2d06e7d7 | 387 | I<phi> are latitudes: zero at the North Pole, northward positive, |

4b0d1da8 | 388 | southward negative. B<NOTE>: this formula thinks in mathematics, not |

2d06e7d7 JH |
389 | geographically: the I<phi> zero is at the North Pole, not at the |

390 | Equator on the west coast of Africa (Bay of Guinea). You need to | |

391 | subtract your geographical coordinates from I<pi/2> (also known as 90 | |

392 | degrees). | |

4b0d1da8 JH |
393 | |

394 | $distance = great_circle_distance($lon0, pi/2 - $lat0, | |

395 | $lon1, pi/2 - $lat1, $rho); | |

396 | ||

51301382 | 397 | =head1 EXAMPLES |

d54bf66f JH |
398 | |

399 | To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N | |

400 | 139.8E) in kilometers: | |

401 | ||

402 | use Math::Trig qw(great_circle_distance deg2rad); | |

403 | ||

404 | # Notice the 90 - latitude: phi zero is at the North Pole. | |

405 | @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); | |

406 | @T = (deg2rad(139.8),deg2rad(90 - 35.7)); | |

407 | ||

408 | $km = great_circle_distance(@L, @T, 6378); | |

409 | ||

4b0d1da8 | 410 | The answer may be off by few percentages because of the irregular |

41bd693c JH |
411 | (slightly aspherical) form of the Earth. The used formula |

412 | ||

413 | lat0 = 90 degrees - phi0 | |

414 | lat1 = 90 degrees - phi1 | |

415 | d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) + | |

416 | sin(lat0) * sin(lat1)) | |

417 | ||

418 | is also somewhat unreliable for small distances (for locations | |

419 | separated less than about five degrees) because it uses arc cosine | |

420 | which is rather ill-conditioned for values close to zero. | |

d54bf66f | 421 | |

5cd24f17 | 422 | =head1 BUGS |

5aabfad6 | 423 | |

5cd24f17 PP |
424 | Saying C<use Math::Trig;> exports many mathematical routines in the |

425 | caller environment and even overrides some (C<sin>, C<cos>). This is | |

426 | construed as a feature by the Authors, actually... ;-) | |

5aabfad6 | 427 | |

5cd24f17 PP |
428 | The code is not optimized for speed, especially because we use |

429 | C<Math::Complex> and thus go quite near complex numbers while doing | |

430 | the computations even when the arguments are not. This, however, | |

431 | cannot be completely avoided if we want things like C<asin(2)> to give | |

432 | an answer instead of giving a fatal runtime error. | |

5aabfad6 | 433 | |

5cd24f17 | 434 | =head1 AUTHORS |

5aabfad6 | 435 | |

ace5de91 | 436 | Jarkko Hietaniemi <F<jhi@iki.fi>> and |

6e238990 | 437 | Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>. |

5aabfad6 PP |
438 | |

439 | =cut | |

440 | ||

441 | # eof |