急》白塞尔大地主题解算反算的问题(150分)

下面是我计算白塞尔大地主题解算的代码：
采用克拉索夫斯基椭球元素
正算结果完全正确，反算结果有很大的出入，大富翁们看看
//copyleft by wangfeng
unit Unit1;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, Grids, StdCtrls, math;

type
TForm1 = class(TForm)
GroupBox1: TGroupBox;
GroupBox2: TGroupBox;
StringGrid1: TStringGrid;
StringGrid2: TStringGrid;
StringGrid3: TStringGrid;
Edit1: TEdit;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
Button1: TButton;
Label5: TLabel;
Label6: TLabel;
Label7: TLabel;
StringGrid4: TStringGrid;
StringGrid5: TStringGrid;
StringGrid6: TStringGrid;
StringGrid7: TStringGrid;
Label8: TLabel;
Label9: TLabel;
Label10: TLabel;
Label11: TLabel;
Button2: TButton;
Label12: TLabel;
Label13: TLabel;
Label14: TLabel;
procedure Button1Click(Sender: TObject);
procedure Button2Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;


var
Form1: TForm1;
const
PI=3.14159265358979324;
e=0.0818133340169315;

implementation

{$R *.dfm}

function HuToJiao(h:do
uble): String;
var
s: string;
d,f: integer;
m,vd:do
uble;
begin

vd:=h*180/PI;
d:=Trunc(vd);
f:=Trunc((vd-d)*60);
m:=((vd-d)*60-f)*60;
s:=IntToStr(d)+'度 '+IntToStr(f)+'分 '+Format('%.5f秒',[m]);
result:=s;
end;



procedure TForm1.Button1Click(Sender: TObject);
var
B1,B2,L1,L2,A12,A21,A0,u1,u2,S:do
uble;
A,B,C,aa,bb:do
uble;
d1,d0,d,d2,t,r:do
uble;
begin

try
B1:=(StrToFloat(StringGrid1.cells[0,0])+StrToFloat(StringGrid1.cells[1,0])/60+StrToFloat(StringGrid1.cells[2,0])/3600)*PI/180;
L1:=(StrToFloat(StringGrid2.cells[0,0])+StrToFloat(StringGrid2.cells[1,0])/60+StrToFloat(StringGrid2.cells[2,0])/3600)*PI/180;
A12:=(StrToFloat(StringGrid3.cells[0,0])+StrToFloat(StringGrid3.cells[1,0])/60+StrToFloat(StringGrid3.cells[2,0])/3600)*PI/180;
S:=StrToFloat(Edit1.Text);
except
ShowMessage('认真一点，拜托！');
exit;
end;

u1:=Arctan(Sqrt(1-e*e)*tan(B1));
A0:=ArcSin(Cos(u1)*Sin(A12));
d1:=Arctan(tan(u1)/Cos(A12));
A:=6356863.020+(10708.949-13.474*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
B:=(5354.469-8.978*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
C:=2.238*Cos(A0)*Cos(A0)*Cos(A0)*Cos(A0)+0.006;
aa:=(33523299-(28189-70*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0))*0.0000000001;
bb:=(0.2907-0.0010*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
//Showmessage(FloatToStr(aa)+#13+FloatToStr(bb));
d0:=(S-(B+C*Cos(2*d1))*Sin(2*d1))/A;
d:=d0+((B+5*C*Cos(2*(d0+d1)))*Sin(2*(d0+d1)))/A;
d2:=d1+d;
t:=Sin(A0)*(aa*d+bb*(Sin(2*d2)-Sin(2*d1)));
u2:=ArcSin(Sin(u1)*Cos(d)+Cos(u1)*Cos(A12)*Sin(d));
B2:=Arctan(Tan(u2)/sqrt(1-e*e));
r:=Arctan(Sin(d)*Sin(A12)/(Cos(u1)*Cos(d)-Sin(u1)*Sin(d)*Cos(A12)));
L2:=L1+r-t;
A21:=Arctan(Cos(u1)*Sin(A12)/(Cos(u1)*Cos(d)*Cos(A12)-Sin(u1)*Sin(d)));
Label5.Caption:='B2: '+HutoJiao(B2);
Label6.Caption:='L2: '+HuToJiao(L2);
Label7.Caption:='A21: '+HuToJiao(A21+3.14159265358979324);
end;




procedure TForm1.Button2Click(Sender: TObject);
var
B1,B2,L1,L2, L,u1,u2,A1,A21,A0,S:do
uble;
r,t,t1,d,d1,d2:do
uble;
A,B,C,aa,bb:do
uble;
begin

try
B1:=(StrToFloat(StringGrid4.cells[0,0])+StrToFloat(StringGrid4.cells[1,0])/60+StrToFloat(StringGrid4.cells[2,0])/3600)*PI/180;
L1:=(StrToFloat(StringGrid5.cells[0,0])+StrToFloat(StringGrid5.cells[1,0])/60+StrToFloat(StringGrid5.cells[2,0])/3600)*PI/180;
B2:=(StrToFloat(StringGrid6.cells[0,0])+StrToFloat(StringGrid6.cells[1,0])/60+StrToFloat(StringGrid6.cells[2,0])/3600)*PI/180;
L2:=(StrToFloat(StringGrid7.cells[0,0])+StrToFloat(StringGrid7.cells[1,0])/60+StrToFloat(StringGrid7.cells[2,0])/3600)*PI/180;
except
ShowMessage('认真一点，拜托！');
exit;
end;

Showmessage(FloatToStr(B1)+#13+FloatToStr(L1)+#13+FloatToStr(B2)+#13+FloatToStr(L2));
u1:=ArcSin(Sqrt(1-e*e)*Sin(B1)/Sqrt(1-e*e*Sin(B1)*Sin(B1)));
u2:=ArcSin(Sqrt(1-e*e)*Sin(B2)/Sqrt(1-e*e*Sin(B2)*Sin(B2)));
L:=L2-L1;
r:=L;
t:=1;
t1:=0;
//Initialize
while (Abs(t1-t)>0.00000000000001)do

begin

A1:=Arctan(Sin(r)*Cos(u2)/(Cos(u1)*Sin(u2)-Sin(u1)*Cos(u2)*Cos(r)));
d:=Arctan((Cos(u2)*Sin(r)*Sin(A1)+Sin(A1)*(Cos(u1)*Sin(u2)-Sin(u1)*Cos(u2)*Cos(r)))/(Sin(u1)*Sin(u2)+Cos(u1)*Cos(u2)*Cos(r)));
d1:=Arctan(Tan(u1)*Tan(A1));
d2:=d1+d;
A0:=ArcSin(Cos(u1)*Sin(A1));
aa:=(33523299-(28189-70*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0))*0.0000000001;
//aa:=691.46768-(0.58143-0.00144*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
bb:=(0.2907-0.0010*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
t:=t1;
t1:=Cos(u1)*Sin(A1)*(aa*d+bb*(Sin(2*d2)-Sin(2*d1)));
ShowMessage(FloatToStr(d)+#13+FloatToStr(t1));
r:=L+t1;
//ShowMessage(FloatToStr(r));
end;

A:=6356863.020+(10708.949-13.474*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
B:=(5354.469-8.978*Cos(A0)*Cos(A0))*Cos(A0)*Cos(A0);
C:=2.238*Cos(A0)*Cos(A0)*Cos(A0)*Cos(A0)+0.006;

S:=A*d+Sin(2*d1)*(B+C*Cos(2*d1))-Sin(2*d2)*(B+C*Cos(2*d2));
A21:=Arctan(Sin(r)*Cos(u1)/(Cos(u1)*Sin(u2)*Cos(r)-Sin(u1)*Cos(u2)));
Label12.Caption:='S: '+FloatToStr(S);
Label13.Caption:='A12: '+HuToJiao(A1);
Label14.Caption:='A21: '+HuToJiao(A21);
end;


end.