Koordinatentransformation GaussKrueger->Laengen/Breitengrad

[brechi]

Hallo,
ich versuch seit einiger Zeit GaussKrueger Koordinaten in andere Formate umzuwandeln (bzw. Umwandling kartesisch, geografisch etc.)

Als Codegrundlage verwende ich die Diplomarbeit von http://home.arcor.de/ernst_werner/diplom/d3.pdf
Zum Vergleich auf Richtigkeit nehme ich dabei die Umgewandelten Koodinaten von MapInfo.

Schon bei der Umwandlung von GaussKrueger in Geografische Koordinaten (Längen/Breitengrade) kommt es zu leichten Abweichungen, im Gegensatz zu dem was mir MapInfo umrechnet. Bei Rückumwandlung dieser in GK hab ich dann Abweichugnen von ca. 400Metern (was natürlich viel zu viel ist).

Der Code:

Delphi-Quellcode:

procedure GkToEll(X, Y, a, e: Extended; var fi, lambda: Extended);
var
  f, n, Grd, sigma, Bf, Nf, Mf, tf, b1, b2, b3, b4, b5, {b6, }est, eta: Extended;

  e2: Extended;
begin
  e2 := e*e;
  est := Sqrt(e2 / (1 - e2));

  f := 1 - Sqrt(1 - (e2));
  n := f / (2-f);
  Grd := (a / (1 + n)) * (1 + 1/4 * Power(n,2) + (1 / 64) * Power(n,4));
  sigma := x / Grd;

  Bf := sigma + (3 / 2) * (n - (9 / 16) * Power(n,3)) * Sin(2 * sigma) + (21 / 16) * Power(n,2) * Sin(4 * sigma) + (151 / 96) * Power(n,3) * Sin(6 * sigma);

  Nf := a / Sqrt(1 - e2 * Power(Sin(Bf),2));
  Mf := Power((a * (1 - e2)) / (1 - e2 * Power(Sin(Bf),2)), (3 / 2));
  tf := Tan(Bf);
  eta := est * Cos(Bf);

  b1 := 1 / (Nf * Cos(Bf));
  b2 := -tf / (2 * Mf * Nf);

  b3 := -(1 + 2 * Power(tf,2) + Power(eta,2)) / (6 * Power(nf,3) * Cos(Bf));

  b4 := (tf / (24 * Mf * Power(nf,3))) * (5 + 3 * Power(tf,2) + Power(eta,2) - 9 * Power(eta,2) * Power(tf,2) - 4 * Power(eta,4));

  b5 := (1 / (120 * Power(nf,5) * Cos(Bf))) * (28 * Power(tf,2) + 24 * Power(tf,4) + 6 * Power(eta,2) + 8 * Power(eta,2) * Power(tf,2));

  fi := Bf + b2 * Power(y,2) + b4 * Power(y,4) {+ b6 * Power(y,6)};
  lambda := b1 * y + b3 * Power(y,3) + b5 * Power(y,5);
end;

var
  Hochwert, Rechtswert: Extended;
  X, Y: Extended;
  bg, lg, hg: Extended;
  L: Integer;

  a, b: Extended;
  linEx, numEx1, numEx2: Extended;

  fi, lambda: Extended;
  alpha, gamma, beta, delta, epsilon, zeta: Extended;
begin
  // GaussKrueger DHDN (Bessel)
  Rechtswert := 3563200.64; // 9.95134518... Bessel laut MapInfo
  Hochwert := 5924051.07; // 53.44584275...  Bessel laut MapInfo
  L := 9; // 9° = GK3

  a := 6377397.155; // Bessel Halbachse groß
  b := 6356078.96282; // Bessel Halbachse klein

  linEx := sqrt(a * a - b * b);
  numEx1 := linEx / a;
  numEx2 := linEx / b;

  Y := Rechtswert - 500000 - L * 1000000 / 3;
  X := Hochwert;

  GkToEll(X, Y, a, numEx1, fi, lambda);

  Writeln(Format('Breitengrad: %12.5f, Laengengrad: %12.5f',[lambda * 180 / pi + L, fi * 180 / pi]));

  // Ausgabe: 9.95132 / 53.44965
end.

Ich finde den verdammten Fehler nicht. Andere Umwandlungen verwenden für GaussKrüger immer zu viele Konstanten. Im Grunde sollten aber die Halbachsen + Offset usw. reichen. Dann kann ich die z.b. für WGS84 und mit der 7-Parametertransformation später leicht austauschen und somit alle denkbaren Umwandlungen vornehmen.