首先，你要确定椭球参数：
C#代码
a=6378140;//西安80椭球IGA75
e2=0.006694384999588;
m0=a*(1-e2);
m2=3.0/2*e2*m0;
m4=5.0/4*e2*m2;
m6=7.0/6*e2*m4;
m8=9.0/8*e2*m6;
a0=m0+m2/2+(3.0/8.0)*m4+(5.0/16.0)*m6+(35.0/128.0)*m8;
a2=m2/2+m4/2+15.0/32*m6+7.0/16*m8;
a4=m4/8+3.0/16*m6+7.0/32*m8;
a6=m6/32+m8/16;
a8=m8/128;
xx=0;
yy=0;
_x=0;
_y=0;
BB=0;
LL=0;
下面才开始正题：
高斯正算：把经纬度坐标转换为平面坐标
C#代码
voidGaussPositive(doubleB,doubleL,doubleL0)
{
doubleX,t,N,h2,l,m,Bmiao,Lmiao;
intBdu,Bfen,Ldu,Lfen;
Bdu=(int)B;
Bfen=(int)(B*100)%100;
Bmiao=(B-Bdu-Bfen*0.01)*10000.0;
B=Bdu*PI/180.0+(Bfen/60.0)*PI/180.0+Bmiao/3600.0*PI/180.0;
Ldu=(int)L;
Lfen=(int)(L*100)%100;
Lmiao=(L-Ldu-Lfen*0.01)*10000.0;
L=Ldu*PI/180.0+(Lfen/60.0)*PI/180+Lmiao/3600.0*PI/180.0;
l=L-L0*PI/180;
X=a0*B-Math.Sin(B)*Math.Cos(B)*((a2-a4+a6)+(2*a4-16.0/3.0*a6)*Math.Sin(B)*Math.Sin(B)+16.0/3.0*a6*Math.Pow(Math.Sin(B),4))+a8/8.0*Math.Sin(8*B);
t=Math.Tan(B);
h2=e2/(1-e2)*Math.Cos(B)*Math.Cos(B);
N=a/Math.Sqrt(1-e2*Math.Sin(B)*Math.Sin(B));
m=Math.Cos(B)*l;
xx=X+N*t*((0.5+(1.0/24.0*(5-t*t+9*h2+4*h2*h2)+1.0/720.0*(61-58*t*t+Math.Pow(t,4))*m*m)*m*m)*m*m);
yy=N*((1+(1.0/6.0*(1-t*t+h2)+1.0/120.0*(5-18*t*t+Math.Pow(t,4)+14*h2-58*h2*t*t)*m*m)*m*m)*m);
yy=yy+500000;
}
高斯反算：把平面坐标转换为经纬度坐标：
C#代码
voidGaussNegative(doublex,doubley,doubleL0)
doubleBf,Vf,l,tf,hf2,Nf,Bmiao,Lmiao;
intBdu,Bfen,Ldu,Lfen;
y=y-500000;
Bf=hcfansuan(x);
Vf=Math.Sqrt(1+e2/(1-e2)*Math.Cos(Bf)*Math.Cos(Bf));
tf=Math.Tan(Bf);
hf2=e2/(1-e2)*Math.Cos(Bf)*Math.Cos(Bf);
Nf=a/Math.Sqrt(1-e2*Math.Sin(Bf)*Math.Sin(Bf));
BB=(Bf-0.5*Vf*Vf*tf*(Math.Pow(y/Nf,2)-1.0/12*(5+3*tf*tf+hf2-9*hf2*tf*tf)*Math.Pow(y/Nf,4)+1.0/360*(61+90*tf*tf+45*tf*tf)*Math.Pow(y/Nf,6)))*180/PI;
Bdu=(int)BB;
Bfen=(int)((BB-Bdu)*60);
Bmiao=((BB-Bdu)*60-Bfen)*60;
BB=Bdu+0.01*Bfen+0.0001*Bmiao;
l=1.0/Math.Cos(Bf)*(y/Nf-1.0/6.0*(1+2*tf*tf+hf2)*Math.Pow(y/Nf,3)+1.0/120.0*(5+28*tf*tf+24*Math.Pow(tf,4)+6*hf2+8*hf2*tf*tf)*Math.Pow(y/Nf,5))*180.0/PI;
LL=L0+l;
Ldu=(int)LL;
Lfen=(int)((LL-Ldu)*60);
Lmiao=((LL-Ldu)*60-Lfen)*60;
LL=Ldu+0.01*Lfen+0.0001*Lmiao;
里面涉及到的弧长反算：
C#代码
doublehcfansuan(doublepX)
{
doubleBf0=pX/a0;
doubleBf1,Bf2;
Bf1=Bf0;
Bf2=(pX-F(Bf1))/a0;
while((Bf2-Bf1)>1.0E-11)
{
Bf1=Bf2;
Bf2=(pX-F(Bf1))/a0;
}
returnBf1;
}
高斯换带就比较简单了：
C#代码
voidGaussZone(doublex,doubley,doubleL0,doubleL0new)
{
GaussNegative(x,y,L0);
GaussPositive(BB,LL,L0new);
}
