www.pudn.com > test_gauss_proj.rar > GaussProj.cpp

//GaussProj.cpp 
#include "stdafx.h" 
#include  
#include "GaussProj.h" 
 
 
//////////////////////////////////////////// 
// Common functions 
//////////////////////////////////////////// 
double ConvertDms2Rad(double Dms) 
{ 
    double Degree; 
    double Miniute; 
    double Second; 
    double Rad; 
    int    Sign; 
     
    if(Dms >= 0) 
    { 
        Sign = 1; 
    } 
    else 
    { 
        Sign = -1; 
    } 
    Dms = fabs(Dms); 
    Degree = floor(Dms); 
    Miniute = floor(fmod(Dms * 100.0, 100.0)); 
    Second = fmod(Dms * 10000.0, 100.0); 
    Rad = Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * PI / 180.0; 
     
    return (Rad); 
} 
 
double ConvertRad2Dms(double Rad) 
{ 
    double Degree; 
    double Miniute; 
    double Second; 
    double Dms; 
    int    Sign; 
     
    if(Rad >= 0) 
    { 
        Sign = 1; 
    } 
    else 
    { 
        Sign = -1; 
    } 
    Rad = fabs(Rad * 180.0 / PI); 
    Degree = floor(Rad); 
    Miniute = floor(fmod(Rad * 60.0, 60.0)); 
    Second = fmod(Rad * 3600.0, 60.0); 
    Dms = Sign * (Degree + Miniute / 100.0 + Second / 10000.0); 
     
    return (Dms); 
} 
 
 
CProjection::CProjection() 
{ 
    vSetEllipsoidType(ELLIPSOID_TYPE_DEFAULT); 
     
    fgSetCurL0(0); 
     
    B = 0; 
    l = 0; 
     
    x = 0; 
    y = 0;   
} 
 
CProjection::~CProjection() 
{ 
 
} 
 
void CProjection::vSetEllipsoidType(ELLIPSOID_TYPE_T eType) 
{ 
	ELLIPSOID_TYPE_T eEllipsoidType; 
    switch(eType) 
    { 
    case ELLIPSOID_TYPE_IAG75: 
    case ELLIPSOID_TYPE_WGS84: 
        eEllipsoidType = eType; 
        break; 
    default: 
        eEllipsoidType = ELLIPSOID_TYPE_KRASSOVSKY; 
        break; 
    } 
     
    //set coefficents in stucture ELLIP_COEFF_T 
    vSetEllipCoeff(eEllipsoidType); 
} 
 
void CProjection::vSetEllipCoeff(ELLIPSOID_TYPE_T eEllipType) 
{ 
    switch(eEllipType) 
    { 
    case ELLIPSOID_TYPE_IAG75: 
        eEllipCoeff.a   = 6378140; 
        eEllipCoeff.A1  = 111133.0046793; 
        eEllipCoeff.A2  = -16038.52818; 
        eEllipCoeff.A3  = 16.83263; 
        eEllipCoeff.A4  = -0.02198; 
        eEllipCoeff.A5  = 0.00003; 
         
        eEllipCoeff.e2  = 0.006694384999588; 
        eEllipCoeff.e12 = 0.006739501819473; 
        break; 
    case ELLIPSOID_TYPE_WGS84: 
        eEllipCoeff.a   = 6378137; 
        eEllipCoeff.A1  = 111132.9525494; 
        eEllipCoeff.A2  = -16038.50840; 
        eEllipCoeff.A3  = 16.83260; 
        eEllipCoeff.A4  = -0.02198; 
        eEllipCoeff.A5  = 0.00003; 
         
        eEllipCoeff.e2  = 0.0066943799013; 
        eEllipCoeff.e12 = 0.00673949674227; 
        break; 
         
    default: 
        //look it as ELLIPSOID_TYPE_KRASSOVSKY 
        eEllipCoeff.a   = 6378245; 
        eEllipCoeff.A1  = 111134.8610828; 
        eEllipCoeff.A2  = -16036.48022; 
        eEllipCoeff.A3  = 16.82805; 
        eEllipCoeff.A4  = -0.02197; 
        eEllipCoeff.A5  = 0.00003; 
         
        eEllipCoeff.e2  = 0.006693421622966; 
        eEllipCoeff.e12 = 0.006738525414683; 
        break; 
    } 
} 
 
//the parameters B&l shoul be in dec. DD.MMSS format! 
//for example: 10.302530 means 10 degrees,30 minutes, 25.30 seconds 
bool CProjection::fgConvertBl2xy(double ConvB, double Convl) 
{ 
    double X, N, t, t2, p, p2, eta2; 
    double sinB, cosB; 
     
    //save current B&l 
    B = ConvertDms2Rad(ConvB); 
    l = ConvertDms2Rad(Convl);     
     
    X = eEllipCoeff.A1 * B * 180.0 / PI + eEllipCoeff.A2 * sin(2 * B) +  
		eEllipCoeff.A3 * sin(4 * B) + eEllipCoeff.A4 * sin(6 * B) +  
		eEllipCoeff.A5 * sin(8 * B); 
    sinB = sin(B); 
    cosB = cos(B); 
    t = tan(B); 
    t2 = t * t; 
    N = eEllipCoeff.a / sqrt(1 - eEllipCoeff.e2 * sinB * sinB); 
    p = cosB * (l - CurL0); 
    p2 = p * p; 
    eta2 = cosB * cosB * eEllipCoeff.e2 / (1 - eEllipCoeff.e2); 
    x = X + N * t * ((0.5 + ((5 - t2 + 9 * eta2 + 4 * eta2 * eta2) / 24.0 +  
		(61 - 58 * t2 + t2 * t2) * p2 / 720.0) * p2) * p2); 
    y = N * p * ( 1 + p2 * ( (1 - t2 + eta2) / 6.0 + p2 * ( 5 - 18 * t2 + t2 * t2  
		+ 14 * eta2 - 58 * eta2 * t2 ) / 120.0)); 
    y += 500000; 
 
    return (true); 
} 
 
//the parameters x&y should be in meter unit! 
bool CProjection::fgConvertxy2Bl(double Convx, double Convy) 
{ 
    double sinB, cosB, t, t2, N ,eta2, V, yN; 
    double preB0, B0; 
    double deta; 
     
    //save current x&y 
    x = Convx; 
    y = Convy; 
     
    y -= 500000; 
    B0 = x / eEllipCoeff.A1; 
     
    do 
    { 
        preB0 = B0; 
        B0 = B0 * PI / 180.0; 
        B0 = (x - (eEllipCoeff.A2 * sin(2 * B0) + eEllipCoeff.A3 * sin(4 * B0) +  
                   eEllipCoeff.A4 * sin(6 * B0) + eEllipCoeff.A5 * sin(8 * B0))) /  
				   eEllipCoeff.A1; 
        deta = fabs(B0 - preB0); 
    }while(deta > 0.000000001); 
     
    B0 = B0 * PI / 180.0; 
    B = ConvertRad2Dms(B0); 
    sinB = sin(B0); 
    cosB = cos(B0); 
    t = tan(B0); 
    t2 = t * t; 
    N = eEllipCoeff.a / sqrt(1 - eEllipCoeff.e2 * sinB * sinB); 
    eta2 = cosB * cosB * eEllipCoeff.e2 / (1 - eEllipCoeff.e2); 
    V = sqrt(1 + eta2); 
    yN = y / N; 
    B = B0 - (yN * yN - (5 + 3 * t2 + eta2 - 9 * eta2 * t2) * yN * yN * yN * yN /  
        12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0)  
        * V * V * t / 2; 
    l = CurL0 + (yN - (1 + 2 * t2 + eta2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24  
        * t2 * t2 + 6 * eta2 + 8 * eta2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB; 
     
    return (true); 
} 
 
bool CProjection::fgGetConvertedxy(double *pConvx, double *pConvy) 
{ 
    if(!pConvx || !pConvy) 
    { 
        return (false);    
    } 
     
    *pConvx = x; 
    *pConvy = y; 
     
    return (true); 
} 
 
bool CProjection::fgGetConvertedBl(double *pConvB, double *pConvl) 
{ 
    if(!pConvB || !pConvl) 
    { 
        return (false);    
    } 
     
    *pConvB = ConvertRad2Dms(B); 
    *pConvl = ConvertRad2Dms(l); 
     
    return (true); 
} 
 
bool CProjection::fgSetCurL0(double L0) 
{ 
    CurL0 = ConvertDms2Rad(L0); 
 
    return(true); 
}