This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math::Complex and Math::Trig updates (Re: [perl #37117] Math::Complex atan2 bug)
[perl5.git] / lib / Math / Trig.t
1 #!./perl 
2
3 #
4 # Regression tests for the Math::Trig package
5 #
6 # The tests here are quite modest as the Math::Complex tests exercise
7 # these interfaces quite vigorously.
8
9 # -- Jarkko Hietaniemi, April 1997
10
11 BEGIN {
12     if ($ENV{PERL_CORE}) {
13         chdir 't' if -d 't';
14         @INC = '../lib';
15     }
16 }
17
18 use Math::Trig 1.03;
19
20 my $pip2 = pi / 2;
21
22 use strict;
23
24 use vars qw($x $y $z);
25
26 my $eps = 1e-11;
27
28 if ($^O eq 'unicos') { # See lib/Math/Complex.pm and t/lib/complex.t.
29     $eps = 1e-10;
30 }
31
32 sub near ($$;$) {
33     my $e = defined $_[2] ? $_[2] : $eps;
34     print "# near? $_[0] $_[1] $e\n";
35     $_[1] ? (abs($_[0]/$_[1] - 1) < $e) : abs($_[0]) < $e;
36 }
37
38 print "1..49\n";
39
40 $x = 0.9;
41 print 'not ' unless (near(tan($x), sin($x) / cos($x)));
42 print "ok 1\n";
43
44 print 'not ' unless (near(sinh(2), 3.62686040784702));
45 print "ok 2\n";
46
47 print 'not ' unless (near(acsch(0.1), 2.99822295029797));
48 print "ok 3\n";
49
50 $x = asin(2);
51 print 'not ' unless (ref $x eq 'Math::Complex');
52 print "ok 4\n";
53
54 # avoid using Math::Complex here
55 $x =~ /^([^-]+)(-[^i]+)i$/;
56 ($y, $z) = ($1, $2);
57 print 'not ' unless (near($y,  1.5707963267949) and
58                      near($z, -1.31695789692482));
59 print "ok 5\n";
60
61 print 'not ' unless (near(deg2rad(90), pi/2));
62 print "ok 6\n";
63
64 print 'not ' unless (near(rad2deg(pi), 180));
65 print "ok 7\n";
66
67 use Math::Trig ':radial';
68
69 {
70     my ($r,$t,$z) = cartesian_to_cylindrical(1,1,1);
71
72     print 'not ' unless (near($r, sqrt(2)))     and
73                         (near($t, deg2rad(45))) and
74                         (near($z, 1));
75     print "ok 8\n";
76
77     ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z);
78
79     print 'not ' unless (near($x, 1)) and
80                         (near($y, 1)) and
81                         (near($z, 1));
82     print "ok 9\n";
83
84     ($r,$t,$z) = cartesian_to_cylindrical(1,1,0);
85
86     print 'not ' unless (near($r, sqrt(2)))     and
87                         (near($t, deg2rad(45))) and
88                         (near($z, 0));
89     print "ok 10\n";
90
91     ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z);
92
93     print 'not ' unless (near($x, 1)) and
94                         (near($y, 1)) and
95                         (near($z, 0));
96     print "ok 11\n";
97 }
98
99 {
100     my ($r,$t,$f) = cartesian_to_spherical(1,1,1);
101
102     print 'not ' unless (near($r, sqrt(3)))     and
103                         (near($t, deg2rad(45))) and
104                         (near($f, atan2(sqrt(2), 1)));
105     print "ok 12\n";
106
107     ($x,$y,$z) = spherical_to_cartesian($r, $t, $f);
108
109     print 'not ' unless (near($x, 1)) and
110                         (near($y, 1)) and
111                         (near($z, 1));
112     print "ok 13\n";
113
114     ($r,$t,$f) = cartesian_to_spherical(1,1,0);
115
116     print 'not ' unless (near($r, sqrt(2)))     and
117                         (near($t, deg2rad(45))) and
118                         (near($f, deg2rad(90)));
119     print "ok 14\n";
120
121     ($x,$y,$z) = spherical_to_cartesian($r, $t, $f);
122
123     print 'not ' unless (near($x, 1)) and
124                         (near($y, 1)) and
125                         (near($z, 0));
126     print "ok 15\n";
127 }
128
129 {
130     my ($r,$t,$z) = cylindrical_to_spherical(spherical_to_cylindrical(1,1,1));
131
132     print 'not ' unless (near($r, 1)) and
133                         (near($t, 1)) and
134                         (near($z, 1));
135     print "ok 16\n";
136
137     ($r,$t,$z) = spherical_to_cylindrical(cylindrical_to_spherical(1,1,1));
138
139     print 'not ' unless (near($r, 1)) and
140                         (near($t, 1)) and
141                         (near($z, 1));
142     print "ok 17\n";
143 }
144
145 {
146     use Math::Trig 'great_circle_distance';
147
148     print 'not '
149         unless (near(great_circle_distance(0, 0, 0, pi/2), pi/2));
150     print "ok 18\n";
151
152     print 'not '
153         unless (near(great_circle_distance(0, 0, pi, pi), pi));
154     print "ok 19\n";
155
156     # London to Tokyo.
157     my @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
158     my @T = (deg2rad(139.8),deg2rad(90 - 35.7));
159
160     my $km = great_circle_distance(@L, @T, 6378);
161
162     print 'not ' unless (near($km, 9605.26637021388));
163     print "ok 20\n";
164 }
165
166 {
167     my $R2D = 57.295779513082320876798154814169;
168
169     sub frac { $_[0] - int($_[0]) }
170
171     my $lotta_radians = deg2rad(1E+20, 1);
172     print "not " unless near($lotta_radians,  1E+20/$R2D);
173     print "ok 21\n";
174
175     my $negat_degrees = rad2deg(-1E20, 1);
176     print "not " unless near($negat_degrees, -1E+20*$R2D);
177     print "ok 22\n";
178
179     my $posit_degrees = rad2deg(-10000, 1);
180     print "not " unless near($posit_degrees, -10000*$R2D);
181     print "ok 23\n";
182 }
183
184 {
185     use Math::Trig 'great_circle_direction';
186
187     print 'not '
188         unless (near(great_circle_direction(0, 0, 0, pi/2), pi));
189     print "ok 24\n";
190
191 # Retired test: Relies on atan2(0, 0), which is not portable.
192 #    print 'not '
193 #       unless (near(great_circle_direction(0, 0, pi, pi), -pi()/2));
194     print "ok 25\n";
195
196     my @London  = (deg2rad(  -0.167), deg2rad(90 - 51.3));
197     my @Tokyo   = (deg2rad( 139.5),   deg2rad(90 - 35.7));
198     my @Berlin  = (deg2rad ( 13.417), deg2rad(90 - 52.533));
199     my @Paris   = (deg2rad (  2.333), deg2rad(90 - 48.867));
200
201     print 'not '
202         unless (near(rad2deg(great_circle_direction(@London, @Tokyo)),
203                      31.791945393073));
204     print "ok 26\n";
205
206     print 'not '
207         unless (near(rad2deg(great_circle_direction(@Tokyo, @London)),
208                      336.069766430326));
209     print "ok 27\n";
210
211     print 'not '
212         unless (near(rad2deg(great_circle_direction(@Berlin, @Paris)),
213                      246.800348034667));
214     
215     print "ok 28\n";
216
217     print 'not '
218         unless (near(rad2deg(great_circle_direction(@Paris, @Berlin)),
219                      58.2079877553156));
220     print "ok 29\n";
221
222     use Math::Trig 'great_circle_bearing';
223
224     print 'not '
225         unless (near(rad2deg(great_circle_bearing(@Paris, @Berlin)),
226                      58.2079877553156));
227     print "ok 30\n";
228
229     use Math::Trig 'great_circle_waypoint';
230     use Math::Trig 'great_circle_midpoint';
231
232     my ($lon, $lat);
233
234     ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.0);
235
236     print 'not ' unless (near($lon, $London[0]));
237     print "ok 31\n";
238
239     print 'not ' unless (near($lat, $pip2 - $London[1]));
240     print "ok 32\n";
241
242     ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 1.0);
243
244     print 'not ' unless (near($lon, $Tokyo[0]));
245     print "ok 33\n";
246
247     print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
248     print "ok 34\n";
249
250     ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.5);
251
252     print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
253     print "ok 35\n";
254
255     print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
256     print "ok 36\n";
257
258     ($lon, $lat) = great_circle_midpoint(@London, @Tokyo);
259
260     print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
261     print "ok 37\n";
262
263     print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
264     print "ok 38\n";
265
266     ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.25);
267
268     print 'not ' unless (near($lon, 0.516073562850837)); # 29.5688 E
269     print "ok 39\n";
270
271     print 'not ' unless (near($lat, 1.170565013391510)); # 67.0684 N
272     print "ok 40\n";
273     ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.75);
274
275     print 'not ' unless (near($lon, 2.17494903805952)); # 124.6154 E
276     print "ok 41\n";
277
278     print 'not ' unless (near($lat, 0.952987032741305)); # 54.6021 N
279     print "ok 42\n";
280
281     use Math::Trig 'great_circle_destination';
282
283     my $dir1 = great_circle_direction(@London, @Tokyo);
284     my $dst1 = great_circle_distance(@London,  @Tokyo);
285
286     ($lon, $lat) = great_circle_destination(@London, $dir1, $dst1);
287
288     print 'not ' unless (near($lon, $Tokyo[0]));
289     print "ok 43\n";
290
291     print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
292     print "ok 44\n";
293
294     my $dir2 = great_circle_direction(@Tokyo, @London);
295     my $dst2 = great_circle_distance(@Tokyo,  @London);
296
297     ($lon, $lat) = great_circle_destination(@Tokyo, $dir2, $dst2);
298
299     print 'not ' unless (near($lon, $London[0]));
300     print "ok 45\n";
301
302     print 'not ' unless (near($lat, $pip2 - $London[1]));
303     print "ok 46\n";
304
305     my $dir3 = (great_circle_destination(@London, $dir1, $dst1))[2];
306
307     print 'not ' unless (near($dir3, 2.69379263839118)); # about 154.343 deg
308     print "ok 47\n";
309
310     my $dir4 = (great_circle_destination(@Tokyo,  $dir2, $dst2))[2];
311
312     print 'not ' unless (near($dir4, 3.6993902625701)); # about 211.959 deg
313     print "ok 48\n";
314
315     print 'not ' unless (near($dst1, $dst2));
316     print "ok 49\n";
317 }
318
319 # eof