File Coverage

blib/lib/GIS/Distance/Vincenty.pm
Criterion Covered Total %
statement 50 50 100.0
branch 1 2 50.0
condition 3 6 50.0
subroutine 6 6 100.0
pod n/a
total 60 64 93.7


line stmt bran cond sub pod time code
1             package GIS::Distance::Vincenty;
2 2     2   38 use 5.008001;
  2         7  
3 2     2   16 use strictures 2;
  2         15  
  2         87  
4             our $VERSION = '0.18';
5              
6 2     2   860 use parent 'GIS::Distance::Formula';
  2         302  
  2         11  
7              
8 2     2   676 use Math::Trig qw( deg2rad pi tan atan asin );
  2         13349  
  2         182  
9 2     2   17 use namespace::clean;
  2         4  
  2         15  
10              
11             sub _distance {
12 6     6   14 my ($lat1, $lon1, $lat2, $lon2) = @_;
13              
14 6 50 33     23 return 0 if (($lon1==$lon2) and ($lat1==$lat2));
15              
16 6         16 $lon1 = deg2rad($lon1);
17 6         63 $lat1 = deg2rad($lat1);
18 6         47 $lon2 = deg2rad($lon2);
19 6         98 $lat2 = deg2rad($lat2);
20              
21 6         44 my($a,$b,$f) = (6378137,6356752.3142,1/298.257223563);
22 6         9 my $l = $lon2 - $lon1;
23 6         20 my $u1 = atan((1-$f) * tan($lat1));
24 6         164 my $u2 = atan((1-$f) * tan($lat2));
25 6         90 my $sin_u1 = sin($u1); my $cos_u1 = cos($u1);
  6         12  
26 6         9 my $sin_u2 = sin($u2); my $cos_u2 = cos($u2);
  6         12  
27 6         9 my $lambda = $l;
28 6         7 my $lambda_pi = 2 * pi;
29 6         11 my $iter_limit = 20;
30 6         12 my($cos_sq_alpha,$sin_sigma,$cos2sigma_m,$cos_sigma,$sigma) = (0,0,0,0,0);
31 6   66     30 while( abs($lambda-$lambda_pi) > 1e-12 && --$iter_limit>0 ){
32 26         40 my $sin_lambda = sin($lambda); my $cos_lambda = cos($lambda);
  26         37  
33 26         55 $sin_sigma = sqrt(($cos_u2*$sin_lambda) * ($cos_u2*$sin_lambda) +
34             ($cos_u1*$sin_u2-$sin_u1*$cos_u2*$cos_lambda) * ($cos_u1*$sin_u2-$sin_u1*$cos_u2*$cos_lambda));
35 26         40 $cos_sigma = $sin_u1*$sin_u2 + $cos_u1*$cos_u2*$cos_lambda;
36 26         38 $sigma = atan2($sin_sigma, $cos_sigma);
37 26         57 my $alpha = asin($cos_u1 * $cos_u2 * $sin_lambda / $sin_sigma);
38 26         137 $cos_sq_alpha = cos($alpha) * cos($alpha);
39 26         44 $cos2sigma_m = $cos_sigma - 2*$sin_u1*$sin_u2/$cos_sq_alpha;
40 26         63 my $cc = $f/16*$cos_sq_alpha*(4+$f*(4-3*$cos_sq_alpha));
41 26         34 $lambda_pi = $lambda;
42 26         97 $lambda = $l + (1-$cc) * $f * sin($alpha) *
43             ($sigma + $cc*$sin_sigma*($cos2sigma_m+$cc*$cos_sigma*(-1+2*$cos2sigma_m*$cos2sigma_m)));
44             }
45 6         16 my $usq = $cos_sq_alpha*($a*$a-$b*$b)/($b*$b);
46 6         19 my $aa = 1 + $usq/16384*(4096+$usq*(-768+$usq*(320-175*$usq)));
47 6         28 my $bb = $usq/1024 * (256+$usq*(-128+$usq*(74-47*$usq)));
48 6         20 my $delta_sigma = $bb*$sin_sigma*($cos2sigma_m+$bb/4*($cos_sigma*(-1+2*$cos2sigma_m*$cos2sigma_m)-
49             $bb/6*$cos2sigma_m*(-3+4*$sin_sigma*$sin_sigma)*(-3+4*$cos2sigma_m*$cos2sigma_m)));
50 6         10 my $c = $b*$aa*($sigma-$delta_sigma);
51              
52 6         29 return $c / 1000;
53             }
54              
55             1;
56             __END__