大地坐标系转换为高斯坐标系的方法-C

1、高斯坐标系的原理，请参阅相关文件。

        我国的地形图采用高斯－克吕格平面直角坐标系。在该坐标系中，横轴：赤道，用Ｙ表示；纵轴：中央经线，用Ｘ表示；坐标原点：中央经线与赤道的交点，用Ｏ表示。赤道以南为负，以北为正；中央经线以东为正，以西为负。每一投影带采用各自独立的高斯-克吕格平面直角坐标系统，并规定y坐标值加500km以避免出现负值。

2、C语言的实现部分如下：（为了便于计算这里将中央经线固定，3/6度带原理自行查阅）

//笛卡尔坐标系
typedef struct tagCRDCARTESIAN
{
 double x;
 double y;
 double z;
}CRDCARTESIAN,*PCRDCARTESIAN;

//大地坐标系
typedef struct tagCRDGEODETIC
{
 double longitude; //经度
 double latitude; //纬度
 double height;    //大地高,可设为0
}CRDGEODETIC;
typedef CRDGEODETIC *PCRDGEODETIC;

#include <math.h>
#include <stdio.h>
/*++
把球面坐标转为平面坐标体系中去采用WGS84标准去转
输入：
[IN]PCRDCARTESIAN pcc 转出的平面坐标系
[IN]PCRDGEODETIC pcg 原始的球面经纬度坐标
返回
--*/
void GeodeticToCartesian (PCRDCARTESIAN pcc, PCRDGEODETIC pcg)
{
 double B; //纬度度数
 double L; //经度度数
 double L0; //中央经线度数
 double l; //L-L0
 double t; //tanB
 double m; //ltanB
 double N; //卯酉圈曲率半径
 double q2;
 double x; //高斯平面纵坐标
 double y; //高斯平面横坐标
 double s; //赤道至纬度B的经线弧长
 double f; //参考椭球体扁率
 double e2; //椭球第一偏心率
 double a; //参考椭球体长半轴
 //double b; //参考椭球体短半轴
 double a1;
 double a2;
 double a3;
 double a4;
 double b1;
 double b2;
 double b3;
 double b4;
 double c0;
 double c1;
 double c2;
 double c3;
 double LLL;
 int Datum=84; //投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面为84  中国2000基准面2000
 int prjno=0; //投影带号
 int zonewide=3;
 double IPI=0.0174532925199433333333; //3.1415926535898/180.0
 B=pcg->latitude ; //纬度
 L=pcg->longitude ; //经度
 
 //LLL = L;
 LLL = 120;//固定度带号
 if (zonewide==6)
 {
  prjno=(int)(LLL/zonewide)+1;
  L0=prjno*zonewide-3;
 }
 else
 {
  prjno=(int)((LLL-1.5)/3)+1;
  L0=prjno*3;
 }
 
 if(Datum==54)
 {
  a=6378245;
  f=1/298.3;
 }
 else if(Datum==84)
 {
  a=6378137;
  f=1/298.257223563;
 }
 else if(Datum==2000)
 {
  a=6378137;
  f=1/298.257222101;
 }

 L0=L0*IPI;
 L=L*IPI;
 B=B*IPI;
 
 e2=2*f-f*f;//(a*a-b*b)/(a*a);
 l=L-L0;
 t=tan(B);
 m=l * cos(B);
 N=a/sqrt(1-e2* sin(B) * sin(B));
 q2=e2/(1-e2)* cos(B)* cos(B);
 a1=1+(double)3/4*e2+(double)45/64*e2*e2+(double)175/256*e2*e2*e2+(double)11025/16384*e2*e2*e2*e2+(double)43659/65536*e2*e2*e2*e2*e2;
 a2=(double)3/4*e2+(double)15/16*e2*e2+(double)525/512*e2*e2*e2+(double)2205/2048*e2*e2*e2*e2+(double)72765/65536*e2*e2*e2*e2*e2;
 a3=(double)15/64*e2*e2+(double)105/256*e2*e2*e2+(double)2205/4096*e2*e2*e2*e2+(double)10359/16384*e2*e2*e2*e2*e2;
 a4=(double)35/512*e2*e2*e2+(double)315/2048*e2*e2*e2*e2+(double)31185/13072*e2*e2*e2*e2*e2;
 b1=a1*a*(1-e2);
 b2=(double)-1/2*a2*a*(1-e2);
 b3=(double)1/4*a3*a*(1-e2);
 b4=(double)-1/6*a4*a*(1-e2);
 c0=b1;
 c1=2*b2+4*b3+6*b4;
 c2=-(8*b3+32*b4);
 c3=32*b4;
 s=c0*B+cos(B)*(c1*sin(B)+c2*sin(B)*sin(B)*sin(B)+c3*sin(B)*sin(B)*sin(B)*sin(B)*sin(B));
 x=s+(double)1/2*N*t*m*m+(double)1/24*(5-t*t+9*q2+4*q2*q2)*N*t*m*m*m*m+(double)1/720*(61-58*t*t+t*t*t*t)*N*t*m*m*m*m*m*m;
 y=N*m+(double)1/6*(1-t*t+q2)*N*m*m*m+(double)1/120*(5-18*t*t+t*t*t*t-14*q2-58*q2*t*t)*N*m*m*m*m*m;
 y=y+1000000*prjno+500000;
 pcc->x=y;
 //pcc->y=y-38000000;
 pcc->y=x;//xy反向
 pcc->z=0;
 printf("safd : %f,%f\n" , pcc->x , pcc->y);
}

void gps_x_y_lo_la(double *x , double *y , double *lo , double *la , unsigned int type)
{
 CRDCARTESIAN pcc;
 CRDGEODETIC pcg;
 if(type == 0)
 {
  pcg.longitude = *lo;
  pcg.latitude = *la;
  GeodeticToCartesian(&pcc , &pcg); 
  *x = pcc.x;
  *y = pcc.y;
 }
}