DZone Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

Snippets has posted 5883 posts at DZone. View Full User Profile

Calculate Latitude And Longitude From Rate Center V&H In C Library

03.26.2010
| 3706 views |
  • submit to reddit
        ZIPCodeWorld.com provides this function to calculate the latitude and longitude coordinates from Vertical and Horizontal (V&H) coordinates in C library. V&H's are used to identify locations and hence relative distances between network elements and between rate centers listed in AreaCodeWorld(TM) North American Area Code Database Subscription which includes NPA (area code), NXX (exchange), country, state, county, latitude, longitude, etc. 

#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);
}