File Coverage

blib/lib/GIS/Distance/ALT.pm
Criterion Covered Total %
statement 49 49 100.0
branch 2 4 50.0
condition n/a
subroutine 7 7 100.0
pod n/a
total 58 60 96.6


line stmt bran cond sub pod time code
1             package GIS::Distance::ALT;
2 1     1   21 use 5.008001;
  1         3  
3 1     1   6 use strictures 2;
  1         7  
  1         39  
4             our $VERSION = '0.19';
5              
6 1     1   719 use parent 'GIS::Distance::Formula';
  1         300  
  1         5  
7              
8 1     1   615 use Math::Trig qw( deg2rad acos pi );
  1         13621  
  1         91  
9 1     1   515 use GIS::Distance::Constants qw( :all );
  1         2  
  1         108  
10 1     1   8 use namespace::clean;
  1         2  
  1         4  
11              
12             my $DEG_RADS = pi / 180;
13              
14             # Sphere with equal meridian length
15             my $RM = 6367449.14582342;
16              
17             # WGS 84 Ellipsoid
18             my $A = 6378137;
19             my $B = 6356752.314245;
20             my $F = 1 / 298.257223563;
21              
22             sub _distance {
23 5     5   11 my ($lat1, $lon1, $lat2, $lon2) = @_;
24              
25 5         11 my $f = 0.5 * ($lat2 + $lat1) * $DEG_RADS;
26 5         10 my $g = 0.5 * ($lat2 - $lat1) * $DEG_RADS;
27 5         7 my $l = 0.5 * ($lon2 - $lon1) * $DEG_RADS;
28              
29 5         25 my $sf = sin($f); my $sg = sin($g); my $sl = sin($l);
  5         9  
  5         10  
30 5         6 my $s2f = $sf * $sf; my $s2g = $sg * $sg; my $s2l = $sl * $sl;
  5         8  
  5         8  
31 5         6 my $c2f = 1.0 - $s2f; my $c2g = 1.0 - $s2g; my $c2l = 1.0 - $s2l;
  5         9  
  5         6  
32              
33 5         9 my $s2 = $s2g * $c2l + $c2f * $s2l;
34 5         7 my $c2 = $c2g * $c2l + $s2f * $s2l;
35              
36 5         11 my ($s, $c, $omega, $rr, $aa, $bb, $pp, $qq, $d2, $qp, $eps1, $eps2);
37              
38 5 50       13 return 0 if $s2 == 0;
39 5 50       9 return pi * $RM * 0.001 if $c2 == 0;
40              
41 5         8 $s = sqrt($s2); $c = sqrt($c2);
  5         6  
42 5         12 $omega = atan2($s, $c);
43 5         7 $rr = $s * $c;
44 5         9 $aa = $s2g * $c2f / $s2 + $s2f * $c2g / $c2;
45 5         6 $bb = $s2g * $c2f / $s2 - $s2f * $c2g / $c2;
46 5         8 $pp = $rr / $omega;
47 5         8 $qq = $omega / $rr;
48 5         7 $d2 = $s2 - $c2;
49 5         10 $qp = $qq + 6 * $pp;
50 5         7 $eps1 = 0.5 * $F * (-$aa - 3 * $bb * $pp);
51 5         13 $eps2 = 0.25 * $F * $F * ((-$qp * $bb + (-3.75 + $d2 * ($qq + 3.75 * $pp)) *
52             $aa + 4. - $d2 * $qq) * $aa - (7.5 * $d2 * $bb * $pp - $qp) * $bb);
53              
54 5         9 my $d = 2 * $omega * $A * (1 + $eps1 + $eps2);
55 5         23 return $d * 0.001;
56             }
57              
58             1;
59             __END__