#! /usr/bin/perl

# Program Name:    ll2vh-openmap.pl

use strict;
use warnings;
use diagnostics;

# Get Latitude and Longitude
my $in1     = $ARGV[0];
my $in2     = $ARGV[1];

my $lat     = $in1;
my $lon     = $in2;

# Load appropriate packages
use Math::Trig;             # Trigonometric functions

# Declare constants

# Polynomial constants
my $k1      =  .99435487;
my $k2      =  .00336523;
my $k3      = -.00065596;
my $k4      =  .00005606;
my $k5      = -.00000188;

#  PI in various forms
my $pi      = 3.14159;
my $m_pi_2  = $pi/2.0;

# spherical coordinates of eastern reference point
#   $ex^2 + $ey^2 + EZ^2 = 1
my $ex      =  .40426992;
my $ey      =  .68210848;
my $ez      =  .60933887;

# spherical coordinates of western reference point
#   WX^2 + WY^2 + WZ^2 = 1
my $wx      =  .65517646;
my $wy      =  .37733790;
my $wz      =  .65449210;

# spherical coordinates of V-H coordinate system
#   PX^2 + PY^2 + PZ^2 = 1
my $px      = -0.555977821730048699;
my $py      = -0.345728488161089920;
my $pz      =  0.755883902605524030;

# Rotation by 76 degrees
my $rot     = deg2rad(76.597497064);
my $rotc    = cos($rot);
my $rots    = sin($rot);

my $radians = deg2rad($pi/180);
my $degrees = rad2deg(180/$pi);

# Orthogonal translation values
my $transv  = 6363.235;
my $transh  = 2250.700;

# Radius of Earth in sqrt(0.1)-mile units, minus 0.3 percent
my $radius  = 12481.103;

my $k9      = $radius*$rotc;
my $k10     = $radius*$rots;

# Variables

my ($e, $w, $x, $y, $z) = 0;
my ($ht, $vt, $h, $v) = 0;
my $cos_lat1 = 0;

# Main program

$lat = deg2rad($lat);
$lon = deg2rad($lon);

# Translate east by 52 degrees #
my $lon1 = $lon + deg2rad(52.0);

# Convert latitude to geocentric latitude using Horner's rule #
my $latsq = $lat*$lat;
my $lat1 = $lat*($k1 + ($k2 + ($k3 + ($k4 + $k5*$latsq)*$latsq)*$latsq
+)*$latsq);

# x, y, and z are the spherical coordinates corresponding to lat, lon.
+ #
$cos_lat1 = cos($lat1);
$x = $cos_lat1*sin($lon1*-1);
$y = $cos_lat1*cos($lon1*-1);
$z = sin($lat1);

# e and w are the cosine of the angular distance (radians) between our
+-
#   point and the east and west centers. #
$e = $ex*$x + $ey*$y + $ez*$z;
$w = $wx*$x + $wy*$y + $wz*$z;

if ($e > 1.0) {
    $e = 1.0;
}

if ($w > 1.0) {
    $w = 1.0;
}

$e = $m_pi_2 - atan($e/sqrt(1 - $e*$e));
$w = $m_pi_2 - atan($w/sqrt(1 - $w*$w));

# e and w are now in radians. #
$ht = ($e*$e - $w*$w + .16)/.8;
$vt = sqrt(abs($e*$e - $ht*$ht));

if (($px*$x + $py*$y + $pz*$z) < 0) {
    $vt = $vt * -1;
}

# rotate and translate to get final v and h. #
$v = $transv +  $k9*$ht - $k10*$vt;
$h = $transh + $k10*$ht +  $k9*$vt;

print "V = $v\tLatitude = $in1\nH = $h\tLongitude = $in2\n\n";
[download]

And here's my test program for your module...

#! /usr/bin/perl

# Program Name:    ll2vh.pl

my $lat= $ARGV[0];
my $lon= $ARGV[1];
use Geo::Coordinates::VandH::XS;
my( $v, $h ) = Geo::Coordinates::VandH::XS::toVH( $lat, $lon );
print "V = $v\tLatitude = $lat\nH = $h\tLongitude = $lon\n\n";
[download]

I tested both programs against the output of this test program for Geo::Coordinates::VandH.

#! /usr/bin/perl

# Program Name:    vh2ll.pl

my $v= $ARGV[0];
my $h= $ARGV[1];
use Geo::Coordinates::VandH;
$blah=new Geo::Coordinates::VandH;
($lat,$lon) = $blah->vh2ll($v,$h);
print "V = $v\tLatitude = $lat\nH = $h\tLongitude = $lon\n\n";
[download]

My test run looks like this

jcoxen@jcoxen vandh>./vh2ll.pl 7350 5909
V = 7350 Latitude = 40.423792190582
H = 5909 Longitude = 104.791265047498

jcoxen@jcoxen vandh>./ll2vh.pl 40.423792190582 104.791265047498
V = -9500.26569570891 Latitude = 40.423792190582
H = 17778.1669205222 Longitude = 104.791265047498

jcoxen@jcoxen vandh>./ll2vh-openmap.pl 40.423792190582 104.791265047498
V = -9500.25261315459 Latitude = 40.423792190582
H = 17778.1566122488 Longitude = 104.791265047498

Any ideas or suggestions? I'm still going over the code to see if I messed up somewhere but so far, I haven't found it. Geo::Coordinates::VandH::XS converts from V&H to Lat/Lon just fine. I just can't get it to go the other way.

BTW, my converted code is very raw. I haven't cleaned out unnecessary lines, prettyied anything up or even put in a big, bold comment stating how I stole re-used this code from OpenMap.

Jack