gaussproj.cpp

这是一个高斯投影变换的类

//GaussProj.cpp
#include "stdafx.h"
#include <math.h>
#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);
}