如何计算高斯投影

废话不多说，直接上代码（中科星图-----合肥）

首先讲高斯正算：

void GaussProjectDirect(double a,double efang,double B,double L,double L0,double& x,double &y,double& R)//高斯投影正算克氏

{

    double b=aefangtob(a,efang);

    double e2=seconde(a,b);

    double W=sqrt(1-efang*sin(B)*sin(B));

    printf("W=%f",W);

    double N=a/W;printf("N=%f",N);

    double M=a*(1-efang)/pow(W,3);

    printf("M=%f",M);

    double t=tee(B);

    double eitef=eitefang(a,b,B);

    double l=L-L0;

    //主曲率半径计算

    double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;

    m0=a*(1-efang);    

    n0=a;

    m2=3.0/2.0*efang*m0;  

    n2=1.0/2.0*efang*n0;

    m4=5.0/4.0*efang*m2;  

    n4=3.0/4.0*efang*n2;

    m6=7.0/6.0*efang*m4;  

    n6=5.0/6.0*efang*n4;

    m8=9.0/8.0*efang*m6;  

    n8=7.0/8.0*efang*n6;

    //子午线曲率半径

    double a0,a2,a4,a6,a8;

    a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;

    a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;

    a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;

    a6=m6/32+m8/16;

    a8=m8/128;

    double X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8*sin(8*B);

    x=X+N/2*t*cos(B)*cos(B)*l*l+N/24*t*(5-t*t+9*eitef+4*pow(eitef,2))*pow(cos(B),4)*pow(l,4)+N/720*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6);

    y=N*cos(B)*l+N/6*(1-t*t+eitef)*pow(cos(B),3)*pow(l,3)+N/120*(5-18*t*t+pow(t,4)+14*eitef-58*eitef*t*t)*pow(cos(B),5)*pow(l,5);

        R=sqrt(M*N);

}