大地坐标换算

// Coords.h
#ifndef __HEADER_COORDS_GEO__
#define  __HEADER_COORDS_GEO__
namespace  Geo
... {

/**//**
* 直角坐标点
**/
struct UPnt 
...{
    /**//**
    * x 坐标值
    **/
    double x;

    /**//**
    * y 坐标值
    **/
    double y;
};

/**//**
* 大地坐标
**/
struct BLCoord
...{

    /**//**
    * 经度,单位为弧度
    **/
    double B;

    /**//**
    * 纬度,单位为弧度
    **/
    double L;
};
}

#endif // __HEADER_COORDS_GEO__

// Ellipsoid.h
#ifndef __HEADER_ELLIPSOID_GEO__
#define  __HEADER_ELLIPSOID_GEO__
#include  " Geo.h "
namespace  Geo
... {

/**//**
* 椭球参数
**/
struct  Ellipsoid
...{
    /**//**
    * 长半坐标轴长
    **/
    double a;

    /**//**
    * 短半坐标轴长
    **/
    double b;

    /**//**
    * 返回扁率
    **/
    double Getf()
    ...{
        return (a-b)/a;

    }

    /**//**
    * 返回第一偏心率
    **/
    double Gete1() const
    ...{
        double rate = b/a;
        return sqrt(1.0-rate*rate);

    }

    /**//**
    * 返回第二偏心率
    **/
    double Gete2() const
    ...{
        double rate = a/b;
        return sqrt(rate*rate-1.0);

    }
    
    /**//**
    * 计算卯酉圈曲率半径
    * @param [in] B 纬度(弧度)
    **/
    double GetN(double B) const
    ...{
        double e1     = Gete1();
        double e1sinB = e1*sin(B);

        return a/sqrt(1 - e1sinB*e1sinB);

    }

    /**//**
    * 计算子午圈曲率半径
    * @param [in] B 纬度(弧度)
    **/

    double GetR(double B) const
    ...{
        double e1     = Gete1();
        double e1sinB = e1*sin(B);
        double div    = sqrt(1 - e1sinB*e1sinB);

        return a*(1-e1*e1)/(div*div*div);

    }

};

}
#endif // __HEADER_ELLIPSOID_GEO__

// Geo.h
#ifndef  __HEADER_GEO__
#define   __HEADER_GEO__

#ifdef    GEO_EXPORTS
     #define  GEO_API _declspec(dllexport)
#else
     #define  GEO_API _declspec(dllimport)
#endif

#include  < math.h >
#include  < assert.h >
#include  < windows.h >
#endif // __HEADER_GEO__

// Prj.h
#ifndef __HEADER_PRJ_GEO__
#define  __HEADER_PRJ_GEO__
#include  " Ellipsoid.h "
#include  " Coords.h "
namespace  Geo
... {

/**//**
* 投影参数
**/
struct Prj
...{
    /**//**
    * 原点纬度
    **/
    double orgB;

    /**//**
    * 中央经度
    **/
    double centerL;

    /**//**
    * 东偏移
    **/ 
    double easteOffset;

    /**//**
    * 北偏移
    **/
    double northOffset;

    /**//**
    * 投影比例因子
    **/
    double k;
    
    /**//**
    * 保留参数
    **/
    double reserved0;
    double reserved1;
    double reserved2;
    
};

/**//**
* 高斯投影正算
* @param [in] ellipsoid  椭球参数
* @param [in] prj        投影参数
* @param [in] blCoords   经纬度坐标数组
* @param [in] count      经纬度坐标数组的长度
* @param [out] uPnts     传出计算的坐标点
**/
extern "C" GEO_API void __stdcall GaussProj(const Ellipsoid* ellipsoid, const Prj* prj, const BLCoord* blCoords, int count, UPnt* uPnts);

/**//**
* 高斯投影反算
* @param [in] ellipsoid  椭球参数
* @param [in] prj        投影参数
* @param [in] blCoords   坐标点数组
* @param [in] count      坐标点数组的长度
* @param [out] uPnts     传出计算的经纬度坐标
**/
extern "C" GEO_API void __stdcall GaussProjInvert(const Ellipsoid* ellipsoid, const Prj* prj, const UPnt* uPnts, int count, BLCoord* blCoords);
}
#endif // __HEADER_PRJ_GEO__

// Publish.h
#ifndef __HEADER_PUBLISH_GEO__
#define  __HEADER_PUBLISH_GEO__
#include  " Geo.h "
#include  " Ellipsoid.h "
#include  " Coords.h "
#include  " Prj.h "
#include  " Util.h "
#endif // __HEADER_PUBLISH_GEO__

// Util.h
#ifndef __HEADER_UTIL_GEO__
#define  __HEADER_UTIL_GEO__
#include  " Geo.h "
namespace  Geo
... {

/**//**
* 弧度转角度(如20.0940即20度,09分40秒)
* @param [in] rad   弧度
**/
extern "C" GEO_API double __stdcall RadToDegree(double rad);

/**//**
* 角度转弧度
* @param [in] degree 角度(格式同上)  
**/
extern "C" GEO_API double __stdcall DegreeToRad(double degree);
}
#endif // __HEADER_UTIL_GEO__

// Geo.cpp
#include  " Geo.h "
BOOL APIENTRY DllMain( HANDLE hModule, 
                       DWORD  ul_reason_for_call, 
                       LPVOID lpReserved
                     )
... {
    return TRUE;
}

// Prj.cpp
#include  " Prj.h "
namespace  Geo
... {

extern "C" GEO_API void __stdcall GaussProj(const Ellipsoid* ellipsoid, const Prj* prj, const BLCoord* blCoords, int count, UPnt* uPnts)
...{
    double e1 = ellipsoid->Gete1();
    double e2 = ellipsoid->Gete2();

    double ee2 = e1*e1;
    double ee4 = ee2*ee2;
    double ee6 = ee2*ee4;

    double k0  = (1 - ee2/4 - 3*ee4/64 - 5*ee6/256);
    double k1  = (-3*ee2/8 - 3*ee4/32 - 45*ee6/1024);
    double k2  = (15*ee4/256 + 45*ee6/1024);
    double k3  = (-35*ee6/3072);
 
    double B, L, T, C, A, M, N;
    for ( int cnt =0; cnt < count; cnt++ )
    ...{
        B = blCoords[cnt].B;
        L = blCoords[cnt].L;
        
        double tanB   = tan(B);
        double cosB   = cos(B);
        double e2cosB = e2*cosB;
        double sin2B  = sin(2*B);
        double sin4B  = sin(4*B);
        double sin6B  = sin(6*B);
        double L0     = prj->centerL;

        T = tanB*tanB;
        C = e2cosB*e2cosB;
        A = (L-L0)*cosB;
        M = ellipsoid->a*(k0*B + k1*sin2B + k2*sin4B + k3*sin6B);
        N = ellipsoid->GetN(B);
        
        double AA2 = A*A;
        double AA4 = AA2*AA2;
        
        uPnts[cnt].x = prj->k*(M + N*tanB*(AA2/2 + (5 - T + 9*C + 4*C*C)*AA4/24) 
            + (61 - 58*T + T*T + 270*C - 330*T*C)*AA2*AA4/720);
        uPnts[cnt].y = prj->easteOffset + prj->k*N*(A + (1 - T + C)*AA2*A/6
            + (5 - 18*T + T*T + 14*C - 58*T*C)*AA4*A/120);


    }

}

extern "C" GEO_API void __stdcall GaussProjInvert(const Ellipsoid* ellipsoid, const Prj* prj, const UPnt* uPnts, int count, BLCoord* blCoords)
...{
    double e1 = ellipsoid->Gete1();
    double e2 = ellipsoid->Gete2();
    double ek = (1 - ellipsoid->b/ellipsoid->a)/(1 + ellipsoid->b/ellipsoid->a);

    double ee2 = e1*e1;
    double ee4 = ee2*ee2;
    double ee6 = ee2*ee4;

    double ek2 = ek*ek;
    double ek3 = ek2*ek;
    double ek4 = ek3*ek;

    double k0  = (3*ek/2 - 27*ek3/32);
    double k1  = (21*ek2/16 - 55*ek4/32);
    double k2  = (151*ek3/96);
    

    double XN, YE, Mf, Bf, Rf, Nf, Tf, Cf, D;
    for ( int cnt =0; cnt < count; cnt++ )
    ...{
        XN = uPnts->x;
        YE = uPnts->y;

        Mf = (XN - prj->northOffset)/prj->k;

        double rad = Mf/(ellipsoid->a*(1 - ee2/4 - 3*ee4/64 - 5*ee6/256));
        double sin2rad  = sin(2*rad);
        double sin4rad  = sin(4*rad);
        double sin6rad  = sin(6*rad);

        Bf = rad + k0*sin2rad + k1*sin4rad + k2*sin6rad;

        double div = sqrt(1 - ee2*sin(Bf)*sin(Bf));
        Rf = ellipsoid->a*(1 - ee2)/(div*div*div);
        Nf = ellipsoid->a/div;

        Tf  = tan(Bf);
        Tf *= Tf;

        Cf  = e2*cos(Bf);
        Cf *= Cf;

        D = (YE - prj->easteOffset)/(prj->k*Nf);
        double D2 = D*D;
        double D4 = D2*D2;
        double D6 = D4*D2;

        blCoords[cnt].B = Bf - (Nf*tan(Bf)/Rf)*(D2/2 - (5 + 3*Tf + Cf - 9*Cf*Tf)*D4/24 + (61 + 90*Tf + 45*Tf*Tf)*D6/720);
        blCoords[cnt].L = prj->centerL + (D -(1 + 2*Tf + Cf)*D2*D/6 + (5 + 28*Tf + 6*Cf + 8*Tf*Cf + 24*Tf*Tf)*D4*D/120)/cos(Bf);

    }

}

}

// Util.cpp
#include  " Util.h "
namespace  Geo
... {
const double GEO_PI = 3.14159265358979323846;
extern "C" GEO_API double __stdcall RadToDegree(double rad)
...{
    /**//*assert(rad);*/

    double res = 0.0;

    double right = (rad*180.0)/GEO_PI;
    int    left   = int(right);
    res += left;

    right  = right - left;
    right *= 60;
    left = int(right);
    res += left/100.0;
    
    right = right -left;
    right *= 60;
    left = int(right);
    res += left/10000.0;

    right = right -left;
    right *= 60;
    left = int(right);
    res += left/1000000.0;

  
    return res;


}

extern "C" GEO_API double __stdcall DegreeToRad(double degree)
...{
    /**//*assert(degree);*/

    int    left  = int(degree);
    double right = degree - left;
    double res = (left/180.0)*GEO_PI;
    

    right *= 100.0;
    left = int(right);
    res += (left/10800.0)*GEO_PI;

    right = right - left;
    right *= 100.0;
    left = int(right);
    res += (left/648000.0)*GEO_PI;

    right = right - left;
    right *= 100.0;
    left = int(right);
    res += (left/38880000.0)*GEO_PI;
    
    return res;


}

}