C语言版经纬度与高斯投影相互转换函数

C语言版经纬度与高斯投影相互转换函数，实现了不同坐标系之间转换

#include “math.h“
#include “coordinate.h“
#include “stdlib.h“


// 高斯投影转经纬度
void GaussToLL（double X double Y double *longitude double *latitudeenum Coordinate coordinate） 

{
	int ProjNo; int ZoneWide; ////带宽 
	double longitude1latitude1 longitude0latitude0 X0Y0 xvalyval;
	double e1e2fa ee NN TC M DRufai iPI;
	iPI = 0.0174532925199433; ////3.1415926535898/180.0; 

	switch （coordinate） 
	{
	case BJ54:
		a=6378245.0; f=1.0/298.3; //54年北京坐标系参数 
		break;
	case WGS84:
		a=6378137.0; f=1/298.257223563; //WGS84
		break;
	case XA80:
		a=6378140.0; f=1/298.257; //80年西安坐标系参数 
		break; 
	}

	//a=6378137.0; f=1/298.257223563; //WGS84
	//a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数 
	////a=6378140.0; f=1/298.257; //80年西安坐标系参数 

	ZoneWide = 6; ////6度带宽 
	ProjNo = （int）（X/1000000L） ; //查找带号
	longitude0 = （ProjNo-1） * ZoneWide + ZoneWide / 2; 
	longitude0 = longitude0 * iPI ; //中央经线
	X0 = ProjNo*1000000L+500000L; 
	Y0 = 0; 
	xval = X-X0; yval = Y-Y0; //带内大地坐标
	e2 = 2*f-f*f;
	e1 = （1.0-sqrt（1-e2））/（1.0+sqrt（1-e2））;
	ee = e2/（1-e2）;
	M = yval;
	u = M/（a*（1-e2/4-3*e2*e2/64-5*e2*e2*e2/256））;
	fai = u+（3*e1/2-27*e1*e1*e1/32）*sin（2*u）+（21*e1*e1/16-55*e1*e1*e1*e1/32）*sin（
		4*u）
		+（151*e1*e1*e1/96）*sin（6*u）+（1097*e1*e1*e1*e1/512）*sin（8*u）;
	C = ee*cos（fai）*cos（fai）;
	T = tan（fai）*tan（fai）;
	NN = a/sqrt（1.0-e2*sin（fai）*sin（fai））;
	R = a*（1-e2）/sqrt（（1-e2*sin（fai）*sin（fai））*（1-e2*sin（fai）*sin（fai））*（1-e2*sin
		（fai）*sin（fai）））;
	D = xval/NN;
	//计算经度（Longitude） 纬度（Latitude）
	longitude1 = longitude0+（D-（1+2*T+C）*D*D*D/6+（5-2*C+28*T-3*C*C+8*ee+24*T*T）*D
		*D*D*D*D/120）/cos（fai）;
	latitude1 = fai -（NN*tan（fai）/R）*（D*D/2-（5+3*T+10*C-4*C*C-9*ee）*D*D*D*D/24
		+（61+90*T+298*C+45*T*T-256*ee-3*C*C）*D*D*D*D*D*D/720）; 
	//转换为度 DD
	*longitude = longitude1 / iPI; 
	*latitude = latitude1 / iPI;
}


//经纬度转高斯投影
void LLToGauss（double longitude double latitude double *X double *Yenum Coordinate coordinate） 
{
	int ProjNo=0; int ZoneWide; ////带宽 
	double longitude1latitude1 longitude0latitude0 xval X0Y0yval;
	double af e2ee NN TCA M iPI;

	iPI = 0.0174532925199433; ////3.1415926535898/180.0; 

	ZoneWide = 6; ////6度带宽 
	
	switch （coordinate） 
	{
	case BJ54:
	    a=6378245.0; f=1.0/298.3; //54年北京坐标系参数 
		break;
	case WGS84:
		a=6378137.0; f=1/298.257223563; //WGS84
		break;
	case XA80:
	    a=6378140.0; f=1/298.257; //80年西安坐标系参数 
		break; 
	}

	ProjNo = （int）（longitude / ZoneWide） ; 
	
	//longitude0 = 0;
	longitude0 = ProjNo * ZoneWide + ZoneWide / 2;    //中央子午线

	longitude0 = longitude0 * iPI ;                   //中央子午线转换为弧度

	latitude0=0; 

	longitude1 = longitude * iPI ; //经度转换为弧度
	latitude1 = latitude * iPI ; //纬度转换为弧度

	e2=2*f-f*f;
	ee=e2*（1.0-e2）;
	NN=a/sqrt（1.0-e2*sin（latitude1）*sin（latitude1））;
	T=tan（latitude1）*tan（latitude1）;
	C=ee*cos（latitude1）*cos（latitude1）;
	A=（longitude1-longitude0）*cos（latitude1）;
	M=a*（（1-e2/4-3*e2*e2/64-5*e2*e2*

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3780  2017-09-27 08:49  coordinate.c

     文件        256  2017-09-27 08:46  coordinate.h

----------- ---------  ---------- -----  ----

                 4036                    2