高斯坐标经纬度互相转换算法(Delphi)

　　这个程序是根据网上找到的VC代码改写而成的Delphi库单元，经验算，比较准确，支持西安80及北京54。
　　本人原创，转载请保留本人信息。http://wallimn.iteye.com
　　代码及使用方法如下（javaeye的博客代码高亮竟然不支持delphi，抗议）：

unit Translate;
{
经纬度坐标与高斯-克吕格投影坐标的互算。
时间：2009-05-11
博客：http://wallimn.iteye.com
转载请保留此信息
}
interface

uses Math;
type
TTranslate = class(TObject)
protected
a,f,e2,e12:double;
A1,A2,A3,A4:double;
private
L0:double; // 中央子午线经度
public
procedure BL2xy(B,L :double; var x,y :double);
procedure xy2BL(x,y :double; var B,L :double);
procedure SetL0(dL0:double);
end;

TTranslate_Krasovsky = class(TTranslate)
public
    constructor Create;
end;

TTranslate_IUGG1975=class(TTranslate)
public
    constructor Create;
end;

function Dms2Rad( Dms:double) : double ;
function Rad2Dms( Rad:double)  : double  ;

implementation
{
将度、分、秒形式转化成弧度
}
function Dms2Rad( Dms:double) : double ;
var
  Degree,Miniute,Second:Double;
  Rad:Double;
  Sign:Integer;
begin
if(Dms >= 0) then
Sign := 1
else
Sign := -1;
Dms := abs(Dms);
Degree := floor(Dms);
Miniute := floor(Dms * 100) mod 100;
Second := floor(Dms * 10000) mod 100;
Rad := Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * PI / 180.0;
result:= Rad;
end;
{
将弧度转换成度、分、秒表示形式；
转换的结果是度，
}
function Rad2Dms( Rad:double)  : double  ;
var
  Degree, Miniute, Second:double;
  Sign:integer;
begin
if(Rad >= 0)     then
Sign := 1
else
Sign := -1;
Rad := abs(Rad * 180.0 / PI);
Degree := floor(Rad);
Miniute := floor(Rad * 60) mod 60;
Second := floor(Rad * 3600) mod 60;
Result := Sign * (Degree + Miniute / 100.0 + Second / 10000.0);
  //Result:=Rad * 180/PI;
end;
{ TTranslate }
{
B,L　为以度为单位的纬度及经度
x,y　为转换结果，即投影坐标，其中y不带带号
时间：2009-05-11
博客：http://wallimn.iteye.com
}
procedure TTranslate.BL2xy(B, L: double; var x, y: double);
var
  XX, N, t, t2, m, m2, ng2:double;
  sinB, cosB:double;
begin
  B:= B*PI/180.0;
  L:= L*PI/180.0;
XX := A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 * B);
sinB := sin(B);
cosB := cos(B);
t := tan(B);
t2 := t * t;
N := a / sqrt(1 - e2 * sinB * sinB);
m := cosB * (L - L0);
m2 := m * m;
ng2 := cosB * cosB * e2 / (1 - e2);
//x,y的计算公式见孔祥元等主编武汉大学出版社2002年出版的《控制测量学》
x := XX + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 -
58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2);
y := N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t2
+ 14 * ng2 - 58 * ng2 * t2 ) / 120.0));
y := y + 500000;
end;

{
设置中央子午线的经度，以度为单位
}
procedure TTranslate.SetL0(dL0: double);
begin
  //L0:= Dms2Rad(dL0);
  L0:=dL0*PI/180.0;
end;
{
x,y　投影坐标，其中y不带带号
B,L　为转换结果，以度为单位的纬度及经度
时间：2009-05-11
博客：http://wallimn.iteye.com
}

procedure TTranslate.xy2BL(x, y: double; var B, L: double);
var
  sinB, cosB, t, t2, N ,ng2, V, yN:double;
  preB0, B0:double;
  eta:double;
begin
y := y- 500000;
B0 := x / A1;
repeat
begin
preB0 := B0;
B0 := B0 * PI / 180.0;
B0 := (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1;
eta := abs(B0 - preB0);
end
  until(eta <= 0.000000001);
B0 := B0 * PI / 180.0;
B := Rad2Dms(B0);
sinB := sin(B0);
cosB := cos(B0);
t := tan(B0);
t2 := t * t;
N := a / sqrt(1 - e2 * sinB * sinB);
ng2 := cosB * cosB * e2 / (1 - e2);
V := sqrt(1 + ng2);
yN := y / N;
B := B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN /
12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0)
* V * V * t / 2;
L := L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24
* t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB;

  //B:=Rad2Dms(B);
  //L:=Rad2Dms(L);
  B:=B*180.0/PI;
  L:=L*180.0/PI;
end;

{ TTranslate_Krasovsky }

constructor TTranslate_Krasovsky.Create;
begin
a := 6378245;
f := 298.3;
e2 := 1 - ((f - 1) / f) * ((f - 1) / f);
e12 := (f / (f - 1)) * (f / (f - 1)) - 1;
A1 := 111134.8611;
A2 := -16036.4803;
A3 := 16.8281;
A4 := -0.0220;
end;

{ TTranslate_IUGG1975 }

constructor TTranslate_IUGG1975.Create;
begin
a := 6378140;
f := 298.257;
e2 := 1 - ((f - 1) / f) * ((f - 1) / f);
e12 := (f / (f - 1)) * (f / (f - 1)) - 1;
A1 := 111133.0047;  //这几个A是什么意思？
A2 := -16038.5282;
A3 := 16.8326;
A4 := -0.0220;
end;
{
引用此库单元，具体使用方法如下：
procedure TFrmMain.Button1Click(Sender: TObject);
var
t:TTranslate;
L,B:Double;
begin
t :=TTranslate_IUGG1975.create;
t.SetL0(111);
t.xy2BL(strToFloat(edtX.text),strToFloat(edtY.text),L,B);
showmessage('L='+FloatToStr(L)+' B='+FloatToStr(B));
//运行结果：L=20,B=109.15
end;

}
end.

注：我的网络硬盘(http://wallimn.ys168.com)上放了这个转换工具，需要的请自行下载