#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#:::                                                                         :::
#:::  TITLE: Calculate Latitude and Longitude from Rate Center V&H in Perl   :::
#:::                                                                         :::
#:::  This function calculates the latitude and longitude coordinates        :::
#:::  from Vertical and Horizontal (V&H) coordinates. V&H's are used to      :::
#:::  identify locations and hence relative distances between network        :::
#:::  elements and between rate centers listed in AreaCodeWorld(tm) Gold     :::
#:::  Edition in http://www.zipcodeworld.com.                                :::
#:::                                                                         :::
#:::  Function Input Parameters:                                             :::
#:::    V = Vertical value from 0 to 10000                                   :::
#:::    H = Horizontal value from 0 to 10000                                 :::
#:::                                                                         :::
#:::  Function Output Parameters:                                            :::
#:::    Lat = Latitude from Vertical value                                   :::
#:::    Long = Longitude from Horizontal value                               :::
#:::                                                                         :::
#:::  North American Area Code NPA NXX database with V&H values is           :::
#:::  available at http://www.zipcodeworld.com. This sample code is          :::
#:::  provided to database subscribers "AS IS" without warranty of any kind. :::
#:::                                                                         :::
#:::  Email: sales@zipcodeworld.com                                          :::
#:::                                                                         :::
#:::  URL:   http://www.zipcodeworld.com                                     :::
#:::                                                                         :::
#:::           ZIPCodeWorld.com © All Rights Reserved 2002-2005              :::
#:::                                                                         :::
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
use Math::Trig;

#my ($lat, $lon) = vh2latlong(5079,1444);

sub vh2latlong {
	my $v = 0;
	my $h = 0;
	($v, $h) = @_;
	my $m_pi = 3.14159265358979323846;
	my $transv = 6363.235;
	my $transh = 2250.7;
	my $rotc = 0.23179040;
	my $rots = 0.97276575;
	my $radius = 12481.103;
	my $ex = 0.40426992;
	my $ey = 0.68210848;
	my $ez = 0.60933887;
	my $wx = 0.65517646;
	my $wy = 0.37733790;
	my $wz = 0.65449210;
	my $px = -0.555977821730048699;
	my $py = -0.345728488161089920;
	my $pz = 0.755883902605524030;
	my $gx = 0.216507961908834992;
	my $gy = -0.134633014879368199;
	my $a = 0.151646645621077297;
	my $q = -0.294355056616412800;
	my $q2 = 0.0866448993556515751;
	my @bi = ( 1.00567724920722457, -0.00344230425560210245, 0.000713971534527667990, -0.0000777240053499279217, 0.00000673180367053244284, -0.000000742595338885741395,  0.0000000905058919926194134 );
	my $x = 0;
	my $y = 0;
	my $z = 0;
	my $delta =0;
	
	my $t1 = ($v - $transv) / $radius;
	my $t2 = ($h - $transh) / $radius;
	my $vhat = $rotc * $t2 - $rots * $t1;
	my $hhat = $rots * $t2 + $rotc * $t1;
	my $e = cos(sqrt($vhat * $vhat + $hhat * $hhat));
	my $w = cos(sqrt($vhat * $vhat + ($hhat - 0.4) * ($hhat - 0.4)));
	my $fx = $ey * $w - $wy * $e;
	my $fy = $ex * $w - $wx * $e;
	my $b = $fx * $gx + $fy * $gy;
	my $c = $fx * $fx + $fy * $fy - $q2;
	my $disc = $b * $b - $a * $c;

	if ($disc == 0.0) {
		$z = $b / $a;
		$x = ($gx * $z - $fx) / $q;
		$y = ($fy - $gy * $z) / $q;
	} else {
		$delta = sqrt($disc); 
		$z = ($b + $delta)/$a;
		$x = ($gx * $z - $fx) / $q;
		$y = ($fy - $gy * $z) / $q;
		if ( $vhat * ( $px * $x + $py * $y + $pz * $z) < 0 ) {
    	$z = ($b - $delta) / $a;
    	$x = ($gx * $z - $fx) / $q;
    	$y = ($fy - $gy * $z) / $q;
    };
  }; 
	my $lat = asin($z);
	my $lat2 = $lat * $lat;
	my $earthlat = 0;
	for (my $i=6; $i>=0 ; $i--) {
		$earthlat = ($earthlat + $bi[$i]) * ($i? $lat2 : $lat);
	};
	$earthlat *= 180/$m_pi;

	my $lon = atan2($x, $y) * 180 / $m_pi;

	my $earthlon = $lon + 52.0000000000000000;
	return ($earthlat, $earthlon);
};