2010. 6. 25. 17:56
갑자기 ASP에서 좌표변환 유틸이 필요하다는 요청에 dll 로 적업하게 됨. 4년만에 Visual studio 를 사용하는 느낌은 너무 생소했다 ㅋ
뭐 이런저런 알고리즘이 있는거 같기는 하지만.. 솔직히 이해는 못하고. 기존에 작업된 자바코드를 토대로 그냥 C++로 변환(??) 했음. 알고리즘을 물어본다면.. 죄송하다는 말밖에는 ㅡㅡ;;
Coordinate.h
Coordinate.cpp
Dll 사용방법
뭐 이런저런 알고리즘이 있는거 같기는 하지만.. 솔직히 이해는 못하고. 기존에 작업된 자바코드를 토대로 그냥 C++로 변환(??) 했음. 알고리즘을 물어본다면.. 죄송하다는 말밖에는 ㅡㅡ;;
Coordinate.h
#include
typedef struct{
double x; //x좌표
double y; //y좌표
} DPoint;
typedef struct{
double lon;
double lat;
double h;
} LonAndLat;
typedef struct{
double x;
double y;
double h;
} TM;
typedef struct{
int ind;
double r_major, r_minor;
double scale_factor;
double lon_center, lat_origin;
double false_northing, false_easting;
double e0, e1, e2, e3;
double e, es, esp, ml0;
} Adjustment;
class CCoordinate {
public:
DPoint convertToKATEC(double dx, double dy);
TM convLLToTM(LonAndLat inputL, int coordinates, int ellipse);
LonAndLat convTMToLL(TM inputT, int inputCoordinates, int targetEllipse);
LonAndLat datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse);
Adjustment setAdjustment(int ellipse, int coordinates);
TM LLtoTM(LonAndLat L, Adjustment A);
LonAndLat TMToLL(TM T, Adjustment A);
LonAndLat convertToWGS(int x, int y);
double mlfn(double e0, double e1, double e2, double e3, double phi);
double e0fn(double x);
double e1fn(double x);
double e2fn(double x);
double e3fn(double x);
double asinz(double con);
double dabs(double x);
BSTR m_x, m_y;
}
typedef struct{
double x; //x좌표
double y; //y좌표
} DPoint;
typedef struct{
double lon;
double lat;
double h;
} LonAndLat;
typedef struct{
double x;
double y;
double h;
} TM;
typedef struct{
int ind;
double r_major, r_minor;
double scale_factor;
double lon_center, lat_origin;
double false_northing, false_easting;
double e0, e1, e2, e3;
double e, es, esp, ml0;
} Adjustment;
class CCoordinate {
public:
DPoint convertToKATEC(double dx, double dy);
TM convLLToTM(LonAndLat inputL, int coordinates, int ellipse);
LonAndLat convTMToLL(TM inputT, int inputCoordinates, int targetEllipse);
LonAndLat datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse);
Adjustment setAdjustment(int ellipse, int coordinates);
TM LLtoTM(LonAndLat L, Adjustment A);
LonAndLat TMToLL(TM T, Adjustment A);
LonAndLat convertToWGS(int x, int y);
double mlfn(double e0, double e1, double e2, double e3, double phi);
double e0fn(double x);
double e1fn(double x);
double e2fn(double x);
double e3fn(double x);
double asinz(double con);
double dabs(double x);
BSTR m_x, m_y;
}
Coordinate.cpp
double DFACT = 36000.000000;
int KATEC = 3;
int WGS84 = 1;
double PI = 3.14159265358979323846;
double minor[]={6356078.96325f, 6356752.3142f};
double major[]={6377397.155f, 6378137.0f};
double dX_W2B = 128;
double dY_W2B = -481;
double dZ_W2B = -664;
int Bessel = 0; //TM은 Bessel 만 사용..
double SF[] = {1, 1, 1, 0.9999, 0.9996, 0.9996};
double FE[] = {200000, 200000, 200000, 400000, 500000, 500000};
double FN[] = {500000, 500000, 500000, 600000, 0.0, 0.0};
double LatCen[] = {0.663225115757845, 0.663225115757845, 0.663225115757845, 0.663225115757845,0, 0};
double LonCen[] = {2.18171200985643, 2.21661859489632, 2.2515251799362, 2.23402144255274, 2.25147473507269, 2.14675497995303};
double EPSLN = 0.0000000001;
DPoint CCoordinate::convertToKATEC(double dx, double dy)
{
DPoint point;
LonAndLat ll;
memset(&point, 0x00, sizeof(DPoint));
memset(&ll, 0x00, sizeof(LonAndLat));
point.x = (double)dx*DFACT;
point.y = (double)dy*DFACT;
double dX = point.x/DFACT;
double dY = point.y/DFACT;
ll.lon = dX;
ll.lat = dY;
ll.h = 0;
TM tm = convLLToTM(ll, KATEC, WGS84);
point.x = (int)tm.x;
point.y = (int)tm.y;
return point;
}
LonAndLat CCoordinate::convertToWGS(int x, int y)
{
// 좌표계 초기화
LonAndLat ll;
TM tm;
memset(&tm, 0x00, sizeof(TM));
memset(&ll, 0x00, sizeof(LonAndLat));
tm.x = (double)x;
tm.y = (double)y;
// 좌표 변환
ll = convTMToLL(tm, KATEC, WGS84);
return ll;
}
LonAndLat CCoordinate::convTMToLL(TM inputT, int inputCoordinates, int targetEllipse) {
TM tm;
memset(&tm, 0x00, sizeof(TM));
tm.x = inputT.x;
tm.y = inputT.y;
// TM좌표 -> 경위도로 변환
//Adjustment A = new Adjustment();
Adjustment A = setAdjustment(Bessel, inputCoordinates);
LonAndLat inputL = TMToLL(tm, A);
//System.out.println("좌표변환후 :"+inputL.lon +" "+inputL.lat);
//System.out.println(inputL.lon*180.0f/Math.PI+" "+inputL.lat*180.0f/Math.PI);
// 타원체를 변경한다.
LonAndLat outputL = datumTrans(inputL, Bessel, targetEllipse);
outputL.lon = outputL.lon*180.0f/PI;
outputL.lat = outputL.lat*180.0f/PI;
return outputL;
}
TM CCoordinate::convLLToTM(LonAndLat inputL, int coordinates, int ellipse) {
LonAndLat LL;
memset(&LL, 0x00, sizeof(LonAndLat));
LL.lon = inputL.lon;
LL.lat = inputL.lat;
LL.h = 0;
//convert degree to radian
LL.lon = inputL.lon * PI / 180.0f;
LL.lat = inputL.lat * PI / 180.0f;
// 타원체를 먼저 변경한다.
LonAndLat outputL = datumTrans(LL, ellipse, Bessel);
//System.out.println("타원체 변경후:"+outputL.lon*180/Math.PI+" "+outputL.lat*180/Math.PI);
//Adjustment A = new Adjustment();
// 경위도 좌표를 TM좌표로 변경한다.
Adjustment A = setAdjustment(Bessel, coordinates);
TM T = LLtoTM(outputL, A);
T.h = outputL.h;
return T;
}
LonAndLat CCoordinate::TMToLL(TM T, Adjustment A) {
//System.out.println("A.e0"+A.e0);
LonAndLat L;
memset(&L, 0x00, sizeof(LonAndLat));
double con, phi, delta_phi, sin_phi, cos_phi, tan_phi;
double c, cs, t, ts, n, r, d, ds, f, h, g, temp;
long max_iter = 6;
long i;
if (A.ind != 0) {
f = exp(T.x / (A.r_major * A.scale_factor));
g = 0.5 * (f - 1.0f / f);
temp = A.lat_origin + T.y / (A.r_major * A.scale_factor);
h = cos(temp);
con = sqrt((1.0f - h * h) / (1.0f + g * g));
L.lat= asinz(con);
if (temp < 0)
L.lat = -L.lat;
if ((g == 0) && (h == 0))
L.lon = A.lon_center;
else
L.lon = atan(g / h) + A.lon_center;
}
// TM to LL inverse equations from here
T.x = T.x - A.false_easting;
T.y = T.y - A.false_northing;
con = (A.ml0 + T.y / A.scale_factor) / A.r_major;
phi = con;
i = 0;
while (true) {
delta_phi = ((con + A.e1 * sin(2.0f * phi) - A.e2 * sin(4.0f * phi)
+ A.e3 * sin(6.0f * phi)) / A.e0) - phi;
phi = phi + delta_phi;
if (dabs(delta_phi) <= EPSLN)
break;
if (i >= max_iter) {
//ShowMessage("Latitude failed to convert...");
// System.out.println("Latitude failed to convert...");
}
i++;
}
if (fabs(phi) < (PI / 2.0f)) {
sin_phi = sin(phi);
cos_phi = cos(phi);
tan_phi = tan(phi);
c = A.esp * (cos_phi * cos_phi);
cs = c * c;
t = tan_phi * tan_phi;
ts = t * t;
con = 1.0f - A.es * sin_phi * sin_phi;
n = A.r_major / sqrt(con);
r = n * (1.0f - A.es) / con;
d = T.x / (n * A.scale_factor);
ds = d * d;
L.lat = phi - (n * tan_phi * ds / r) * (0.5f - ds / 24.0f *
(5.0f + 3.0f * t + 10.0f * c - 4.0f * cs - 9.0f * A.esp - ds / 30.0f *
(61.0f + 90.0f * t + 298.0f * c + 45.0f * ts - 252.0f * A.esp - 3.0f * cs)));
L.lon = A.lon_center + (d *
(1.0f - ds / 6.0f *
(1.0f + 2.0f * t + c - ds / 20.0f *
(5.0f - 2.0f * c + 28.0f * t - 3.0f * cs + 8.0f * A.esp + 24.0f * ts))) / cos_phi);
// *lat = 127.2866484932;
// *lon = 37.4402108468;
} else {
L.lat = PI/2.0f * sin(T.y);
L.lon = A.lon_center;
}
return L;
}
LonAndLat CCoordinate::datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse) {
LonAndLat outputL;
memset(&outputL, 0x00, sizeof(LonAndLat));
double input_a = major[inputEllipse];
double input_b = minor[inputEllipse];
double output_a = major[outputEllipse];
double output_b = minor[outputEllipse];
double delta_a, delta_f, delta_Phi, delta_Lamda, delta_H;
double Rm, Rn;
double temp, es_temp;
int gap = inputEllipse - outputEllipse;
double delta_X = gap*dX_W2B;
double delta_Y = gap*dY_W2B;
double delta_Z = gap*dZ_W2B;
temp = input_b / input_a;
es_temp = 1.0f - temp*temp;
delta_a = output_a - input_a;
delta_f = input_b / input_a - output_b / output_a;
Rm =(input_a * (1.0f - es_temp)) / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (3.0f/2.0f));
Rn = input_a / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (1.0f/ 2.0f));
delta_Phi = ((((-delta_X * sin(inputL.lat) * cos(inputL.lon) - delta_Y * sin(inputL.lat) * sin(inputL.lon))
+ delta_Z * cos(inputL.lat)) + delta_a * Rn * es_temp * sin(inputL.lat) * cos(inputL.lat) / input_a)
+ delta_f * (Rm / temp + Rn * temp) * sin(inputL.lat) * cos(inputL.lat)) / (Rm + inputL.h);
//System.out.println(delta_Phi);
delta_Lamda = ((-delta_X) * sin(inputL.lon) + delta_Y * cos(inputL.lon)) / ((Rn + inputL.h) * cos(inputL.lat));
delta_H = delta_X * cos(inputL.lat) * cos(inputL.lon) + delta_Y * cos(inputL.lat) * sin(inputL.lon)
+ delta_Z * sin(inputL.lat) - delta_a * input_a / Rn + delta_f * temp * Rn * sin(inputL.lat) * sin(inputL.lat);
//System.out.println(delta_Lamda);
outputL.lat = inputL.lat + delta_Phi;
outputL.lon = inputL.lon + delta_Lamda;
outputL.h = inputL.h + delta_H;
return outputL;
}
Adjustment CCoordinate::setAdjustment(int ellipse, int coordinates) {
Adjustment A;
memset(&A, 0x00, sizeof(Adjustment));
A.r_major = major[ellipse];
A.r_minor = minor[ellipse];
A.scale_factor = SF[coordinates];
A.lon_center = LonCen[coordinates];
A.lat_origin = LatCen[coordinates];
A.false_northing = FN[coordinates];
A.false_easting = FE[coordinates];
double temp = A.r_minor / A.r_major;
A.es = 1.0f - temp*temp;
A.e = sqrt(A.es);
A.e0 = e0fn(A.es);
A.e1 = e1fn(A.es);
A.e2 = e2fn(A.es);
A.e3 = e3fn(A.es);
A.ml0 = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, A.lat_origin);
A.esp = A.es / (1.0f - A.es);
if(A.es < 0.00001f)
A.ind = 1;
else
A.ind = 0;
return A;
}
TM CCoordinate::LLtoTM(LonAndLat L, Adjustment A)
{
TM T;
memset(&T, 0x00, sizeof(TM));
double delta_lon;
double sin_phi, cos_phi;
double al, als, b=0, c, t, tq;
double con, n, ml;
// LL to TM forward equations from here;
delta_lon = L.lon - A.lon_center;
sin_phi = sin(L.lat);
cos_phi = cos(L.lat);
if (A.ind != 0) {
b = cos_phi * sin(delta_lon);
if ((fabs(fabs(b) - 1.0f)) < 0.0000000001f) {
//ShowMessage("The point is going to the unlimited value...\n");
//System.out.println("The point is going to the unlimited value...");
}
} else {
T.x = 0.5f * A.r_major * A.scale_factor * log((1.0f + b) / 1.0f- b);
con = acos(cos_phi * cos(delta_lon) / sqrt(1.0f - b * b));
if (L.lat < 0) {
con = -con;
T.y = A.r_major * A.scale_factor * (con - A.lat_origin);
}
}
al = cos_phi * delta_lon;
als = al * al;
c = A.esp * cos_phi * cos_phi;
tq = tan(L.lat);
t = tq * tq;
con = 1.0f - A.es * sin_phi * sin_phi;
n = A.r_major / sqrt(con);
ml = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, L.lat);
T.x = A.scale_factor * n * al * (1.0f + als / 6.0f *
(1.0f - t + c + als / 20.0f *
(5.0f - 18.0f * t + t * t + 72.0f * c - 58.0f * A.esp)
)
) + A.false_easting;
T.y = A.scale_factor * (ml - A.ml0 + n * tq *
(als *
(0.5f + als / 24.0f *
(5.0f - t + 9.0f * c + 4.0f * c * c + als / 30.0f *
(61.0f - 58.0f * t + t * t + 600.0f * c - 330.0f * A.esp)
)))) + A.false_northing;
return T;
}
double CCoordinate::mlfn(double e0, double e1, double e2, double e3, double phi)
{
return (e0 * phi - e1 * sin(2.0f * phi) + e2 * sin(4.0f * phi) - e3 * sin(6.0f * phi));
}
double CCoordinate::e0fn(double x)
{
return (1.0f - 0.25f*x*(1.0f + x / 16.0f * (3+1.25f*x)));
}
double CCoordinate::e1fn(double x)
{
return(0.375f * x * (1. + 0.25f*x*(1.0f + 0.46875f*x)));
}
double CCoordinate::e2fn(double x)
{
return(0.05859375f*x*x*(1.0f+0.75f*x));
}
double CCoordinate::e3fn(double x)
{
return (x * x * x * (35.0f / 3072.0f));
}
double CCoordinate::asinz(double con)
{
if (fabs(con) > 1.0f)
if (con > 1.0f)
con = 1.0f;
else
con = -1.0f;
return asin(con);
}
double CCoordinate::dabs(double x)
{
if(x>=0) return x;
else return -x;
}
int KATEC = 3;
int WGS84 = 1;
double PI = 3.14159265358979323846;
double minor[]={6356078.96325f, 6356752.3142f};
double major[]={6377397.155f, 6378137.0f};
double dX_W2B = 128;
double dY_W2B = -481;
double dZ_W2B = -664;
int Bessel = 0; //TM은 Bessel 만 사용..
double SF[] = {1, 1, 1, 0.9999, 0.9996, 0.9996};
double FE[] = {200000, 200000, 200000, 400000, 500000, 500000};
double FN[] = {500000, 500000, 500000, 600000, 0.0, 0.0};
double LatCen[] = {0.663225115757845, 0.663225115757845, 0.663225115757845, 0.663225115757845,0, 0};
double LonCen[] = {2.18171200985643, 2.21661859489632, 2.2515251799362, 2.23402144255274, 2.25147473507269, 2.14675497995303};
double EPSLN = 0.0000000001;
DPoint CCoordinate::convertToKATEC(double dx, double dy)
{
DPoint point;
LonAndLat ll;
memset(&point, 0x00, sizeof(DPoint));
memset(&ll, 0x00, sizeof(LonAndLat));
point.x = (double)dx*DFACT;
point.y = (double)dy*DFACT;
double dX = point.x/DFACT;
double dY = point.y/DFACT;
ll.lon = dX;
ll.lat = dY;
ll.h = 0;
TM tm = convLLToTM(ll, KATEC, WGS84);
point.x = (int)tm.x;
point.y = (int)tm.y;
return point;
}
LonAndLat CCoordinate::convertToWGS(int x, int y)
{
// 좌표계 초기화
LonAndLat ll;
TM tm;
memset(&tm, 0x00, sizeof(TM));
memset(&ll, 0x00, sizeof(LonAndLat));
tm.x = (double)x;
tm.y = (double)y;
// 좌표 변환
ll = convTMToLL(tm, KATEC, WGS84);
return ll;
}
LonAndLat CCoordinate::convTMToLL(TM inputT, int inputCoordinates, int targetEllipse) {
TM tm;
memset(&tm, 0x00, sizeof(TM));
tm.x = inputT.x;
tm.y = inputT.y;
// TM좌표 -> 경위도로 변환
//Adjustment A = new Adjustment();
Adjustment A = setAdjustment(Bessel, inputCoordinates);
LonAndLat inputL = TMToLL(tm, A);
//System.out.println("좌표변환후 :"+inputL.lon +" "+inputL.lat);
//System.out.println(inputL.lon*180.0f/Math.PI+" "+inputL.lat*180.0f/Math.PI);
// 타원체를 변경한다.
LonAndLat outputL = datumTrans(inputL, Bessel, targetEllipse);
outputL.lon = outputL.lon*180.0f/PI;
outputL.lat = outputL.lat*180.0f/PI;
return outputL;
}
TM CCoordinate::convLLToTM(LonAndLat inputL, int coordinates, int ellipse) {
LonAndLat LL;
memset(&LL, 0x00, sizeof(LonAndLat));
LL.lon = inputL.lon;
LL.lat = inputL.lat;
LL.h = 0;
//convert degree to radian
LL.lon = inputL.lon * PI / 180.0f;
LL.lat = inputL.lat * PI / 180.0f;
// 타원체를 먼저 변경한다.
LonAndLat outputL = datumTrans(LL, ellipse, Bessel);
//System.out.println("타원체 변경후:"+outputL.lon*180/Math.PI+" "+outputL.lat*180/Math.PI);
//Adjustment A = new Adjustment();
// 경위도 좌표를 TM좌표로 변경한다.
Adjustment A = setAdjustment(Bessel, coordinates);
TM T = LLtoTM(outputL, A);
T.h = outputL.h;
return T;
}
LonAndLat CCoordinate::TMToLL(TM T, Adjustment A) {
//System.out.println("A.e0"+A.e0);
LonAndLat L;
memset(&L, 0x00, sizeof(LonAndLat));
double con, phi, delta_phi, sin_phi, cos_phi, tan_phi;
double c, cs, t, ts, n, r, d, ds, f, h, g, temp;
long max_iter = 6;
long i;
if (A.ind != 0) {
f = exp(T.x / (A.r_major * A.scale_factor));
g = 0.5 * (f - 1.0f / f);
temp = A.lat_origin + T.y / (A.r_major * A.scale_factor);
h = cos(temp);
con = sqrt((1.0f - h * h) / (1.0f + g * g));
L.lat= asinz(con);
if (temp < 0)
L.lat = -L.lat;
if ((g == 0) && (h == 0))
L.lon = A.lon_center;
else
L.lon = atan(g / h) + A.lon_center;
}
// TM to LL inverse equations from here
T.x = T.x - A.false_easting;
T.y = T.y - A.false_northing;
con = (A.ml0 + T.y / A.scale_factor) / A.r_major;
phi = con;
i = 0;
while (true) {
delta_phi = ((con + A.e1 * sin(2.0f * phi) - A.e2 * sin(4.0f * phi)
+ A.e3 * sin(6.0f * phi)) / A.e0) - phi;
phi = phi + delta_phi;
if (dabs(delta_phi) <= EPSLN)
break;
if (i >= max_iter) {
//ShowMessage("Latitude failed to convert...");
// System.out.println("Latitude failed to convert...");
}
i++;
}
if (fabs(phi) < (PI / 2.0f)) {
sin_phi = sin(phi);
cos_phi = cos(phi);
tan_phi = tan(phi);
c = A.esp * (cos_phi * cos_phi);
cs = c * c;
t = tan_phi * tan_phi;
ts = t * t;
con = 1.0f - A.es * sin_phi * sin_phi;
n = A.r_major / sqrt(con);
r = n * (1.0f - A.es) / con;
d = T.x / (n * A.scale_factor);
ds = d * d;
L.lat = phi - (n * tan_phi * ds / r) * (0.5f - ds / 24.0f *
(5.0f + 3.0f * t + 10.0f * c - 4.0f * cs - 9.0f * A.esp - ds / 30.0f *
(61.0f + 90.0f * t + 298.0f * c + 45.0f * ts - 252.0f * A.esp - 3.0f * cs)));
L.lon = A.lon_center + (d *
(1.0f - ds / 6.0f *
(1.0f + 2.0f * t + c - ds / 20.0f *
(5.0f - 2.0f * c + 28.0f * t - 3.0f * cs + 8.0f * A.esp + 24.0f * ts))) / cos_phi);
// *lat = 127.2866484932;
// *lon = 37.4402108468;
} else {
L.lat = PI/2.0f * sin(T.y);
L.lon = A.lon_center;
}
return L;
}
LonAndLat CCoordinate::datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse) {
LonAndLat outputL;
memset(&outputL, 0x00, sizeof(LonAndLat));
double input_a = major[inputEllipse];
double input_b = minor[inputEllipse];
double output_a = major[outputEllipse];
double output_b = minor[outputEllipse];
double delta_a, delta_f, delta_Phi, delta_Lamda, delta_H;
double Rm, Rn;
double temp, es_temp;
int gap = inputEllipse - outputEllipse;
double delta_X = gap*dX_W2B;
double delta_Y = gap*dY_W2B;
double delta_Z = gap*dZ_W2B;
temp = input_b / input_a;
es_temp = 1.0f - temp*temp;
delta_a = output_a - input_a;
delta_f = input_b / input_a - output_b / output_a;
Rm =(input_a * (1.0f - es_temp)) / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (3.0f/2.0f));
Rn = input_a / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (1.0f/ 2.0f));
delta_Phi = ((((-delta_X * sin(inputL.lat) * cos(inputL.lon) - delta_Y * sin(inputL.lat) * sin(inputL.lon))
+ delta_Z * cos(inputL.lat)) + delta_a * Rn * es_temp * sin(inputL.lat) * cos(inputL.lat) / input_a)
+ delta_f * (Rm / temp + Rn * temp) * sin(inputL.lat) * cos(inputL.lat)) / (Rm + inputL.h);
//System.out.println(delta_Phi);
delta_Lamda = ((-delta_X) * sin(inputL.lon) + delta_Y * cos(inputL.lon)) / ((Rn + inputL.h) * cos(inputL.lat));
delta_H = delta_X * cos(inputL.lat) * cos(inputL.lon) + delta_Y * cos(inputL.lat) * sin(inputL.lon)
+ delta_Z * sin(inputL.lat) - delta_a * input_a / Rn + delta_f * temp * Rn * sin(inputL.lat) * sin(inputL.lat);
//System.out.println(delta_Lamda);
outputL.lat = inputL.lat + delta_Phi;
outputL.lon = inputL.lon + delta_Lamda;
outputL.h = inputL.h + delta_H;
return outputL;
}
Adjustment CCoordinate::setAdjustment(int ellipse, int coordinates) {
Adjustment A;
memset(&A, 0x00, sizeof(Adjustment));
A.r_major = major[ellipse];
A.r_minor = minor[ellipse];
A.scale_factor = SF[coordinates];
A.lon_center = LonCen[coordinates];
A.lat_origin = LatCen[coordinates];
A.false_northing = FN[coordinates];
A.false_easting = FE[coordinates];
double temp = A.r_minor / A.r_major;
A.es = 1.0f - temp*temp;
A.e = sqrt(A.es);
A.e0 = e0fn(A.es);
A.e1 = e1fn(A.es);
A.e2 = e2fn(A.es);
A.e3 = e3fn(A.es);
A.ml0 = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, A.lat_origin);
A.esp = A.es / (1.0f - A.es);
if(A.es < 0.00001f)
A.ind = 1;
else
A.ind = 0;
return A;
}
TM CCoordinate::LLtoTM(LonAndLat L, Adjustment A)
{
TM T;
memset(&T, 0x00, sizeof(TM));
double delta_lon;
double sin_phi, cos_phi;
double al, als, b=0, c, t, tq;
double con, n, ml;
// LL to TM forward equations from here;
delta_lon = L.lon - A.lon_center;
sin_phi = sin(L.lat);
cos_phi = cos(L.lat);
if (A.ind != 0) {
b = cos_phi * sin(delta_lon);
if ((fabs(fabs(b) - 1.0f)) < 0.0000000001f) {
//ShowMessage("The point is going to the unlimited value...\n");
//System.out.println("The point is going to the unlimited value...");
}
} else {
T.x = 0.5f * A.r_major * A.scale_factor * log((1.0f + b) / 1.0f- b);
con = acos(cos_phi * cos(delta_lon) / sqrt(1.0f - b * b));
if (L.lat < 0) {
con = -con;
T.y = A.r_major * A.scale_factor * (con - A.lat_origin);
}
}
al = cos_phi * delta_lon;
als = al * al;
c = A.esp * cos_phi * cos_phi;
tq = tan(L.lat);
t = tq * tq;
con = 1.0f - A.es * sin_phi * sin_phi;
n = A.r_major / sqrt(con);
ml = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, L.lat);
T.x = A.scale_factor * n * al * (1.0f + als / 6.0f *
(1.0f - t + c + als / 20.0f *
(5.0f - 18.0f * t + t * t + 72.0f * c - 58.0f * A.esp)
)
) + A.false_easting;
T.y = A.scale_factor * (ml - A.ml0 + n * tq *
(als *
(0.5f + als / 24.0f *
(5.0f - t + 9.0f * c + 4.0f * c * c + als / 30.0f *
(61.0f - 58.0f * t + t * t + 600.0f * c - 330.0f * A.esp)
)))) + A.false_northing;
return T;
}
double CCoordinate::mlfn(double e0, double e1, double e2, double e3, double phi)
{
return (e0 * phi - e1 * sin(2.0f * phi) + e2 * sin(4.0f * phi) - e3 * sin(6.0f * phi));
}
double CCoordinate::e0fn(double x)
{
return (1.0f - 0.25f*x*(1.0f + x / 16.0f * (3+1.25f*x)));
}
double CCoordinate::e1fn(double x)
{
return(0.375f * x * (1. + 0.25f*x*(1.0f + 0.46875f*x)));
}
double CCoordinate::e2fn(double x)
{
return(0.05859375f*x*x*(1.0f+0.75f*x));
}
double CCoordinate::e3fn(double x)
{
return (x * x * x * (35.0f / 3072.0f));
}
double CCoordinate::asinz(double con)
{
if (fabs(con) > 1.0f)
if (con > 1.0f)
con = 1.0f;
else
con = -1.0f;
return asin(con);
}
double CCoordinate::dabs(double x)
{
if(x>=0) return x;
else return -x;
}
Dll 사용방법
<%
Set trans = Server.CreateObject("Interactive.Coordinate")
// WGS 타입을 KATEC 타입으로 변환.
trans.WGS2KATEC 126.83530, 37.62508
Response.write trans.X
Response.write " : "
Response.write trans.Y
// KATEC 타입을 WGS 타입으로 변화
trans.KATEC2WGS 297386, 558722
Response.write " "
Response.write trans.X
Response.write " : "
Response.write trans.Y
%>
Set trans = Server.CreateObject("Interactive.Coordinate")
// WGS 타입을 KATEC 타입으로 변환.
trans.WGS2KATEC 126.83530, 37.62508
Response.write trans.X
Response.write " : "
Response.write trans.Y
// KATEC 타입을 WGS 타입으로 변화
trans.KATEC2WGS 297386, 558722
Response.write " "
Response.write trans.X
Response.write " : "
Response.write trans.Y
%>