/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
:                                                                       :
:  TITLE: Calculate Latitude and Longitude from Rate Center V&H in C    :
:                                                                       :
:  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-2010             :
:                                                                       :
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/


#include <math.h>

/* constants contributed by vh2ll.c */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define TRANSV 6363.235
#define TRANSH 2250.7

#define ROTC 0.23179040
#define ROTS 0.97276575

#define RADIUS 12481.103

#define EX 0.40426992
#define EY 0.68210848
#define EZ 0.60933887

#define WX 0.65517646
#define WY 0.37733790
#define WZ 0.65449210

#define PX -0.555977821730048699
#define PY -0.345728488161089920
#define PZ  0.755883902605524030

#define GX  0.216507961908834992
#define GY -0.134633014879368199

#define A  0.151646645621077297

#define Q -0.294355056616412800
#define Q2  0.0866448993556515751

static void
vh2latlong(v, h)
double v, h;
{
	int i=0, latdeg=0, latmin=0, londeg=0, lonmin=0, latsec=0, lonsec=0;
	double t1=0.0, t2=0.0, vhat=0.0, hhat=0.0, fx=0.0, fy=0.0;
	double e=0.0, w=0.0;
	double b=0.0, c=0.0, disc=0.0, z=0.0, x=0.0, y=0.0, delta=0.0, lat=0.0, lat2=0.0, lon=0.0;
	double earthlat=0.0, earthlon=0.0;
	static double bi[7] = {
		1.00567724920722457,
		-0.00344230425560210245,
		0.000713971534527667990,
		-0.0000777240053499279217,
		0.00000673180367053244284,
		-0.000000742595338885741395,
		0.0000000905058919926194134
	};

	t1 = (v - TRANSV) / RADIUS;
	t2 = (h - TRANSH) / RADIUS;
	vhat = ROTC*t2 - ROTS*t1;
	hhat = ROTS*t2 + ROTC*t1;
	e = cos(sqrt(vhat*vhat + hhat*hhat));
	w = cos(sqrt(vhat*vhat + (hhat-0.4)*(hhat-0.4)));
	fx = EY*w - WY*e;
	fy = EX*w - WX*e;
	b = fx*GX + fy*GY;
	c = fx*fx + fy*fy - Q2;
	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;
		}
	}
	lat = asin(z);

	lat2 = lat*lat;
	earthlat = 0.0;
	for (i=6; i>=0; i--) {
		earthlat = (earthlat + bi[i]) * (i? lat2: lat);
	}
	earthlat *= 180/M_PI;
	lon = atan2(x, y) * 180/M_PI;
	earthlon = lon + 52;
	printf("%lf %lf\n", earthlat, earthlon);
}