[Koleksi Kode Sumber] (10) : Menghitung koordinat latitude dan longitude dari koordinat VH

/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
: :
: 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 :
: :
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
[/sourcecode]


[sourcecode language="cpp"]
#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);
}