Geographische in Gauß-Krüger-Koordinaten umrechnen

type Twinkel = record
  grad : integer;
  min : integer;
  sek : integer;
end;

{==============================================================================}
{ Umrechnung von Geographischen Koordinaten in Gauss-Krueger-Koordinaten       }
{ Formel: Grossmann,W., Geodätische Abbildungen, 1964, Seite 151               }
{ Parameter: geo.Breite (Grad.Min.Sek) in Altgrad  : Twinkel                   }
{            geo.Laenge (Grad.Min.Sek) in Altgrad  : Twinkel                   }
{            Zielsystemnummer (Meridiankennziffer) : longint                   }
{            Rechtswert (X) im Zielsystem          : double                    }
{            Hochwert (Y) im Zielsystem            : double                    }
{==============================================================================}
procedure GeoGk(br,la:Twinkel;sy:Longint;var x,y:double);
const
  {26}
  rho = 180 / pi;
var
  brDezimal,laDezimal,rm,e2,c,bf,g,co,g2,g1,t,dl,fa,grad,min,sek :extended;
begin
  {25}
  e2 := 0.0067192188;
  {27}
  c := 6398786.849;
  {in Dezimal}
  {Breite}
  brDezimal := br.grad + br.min / 60 + br.sek / 3600;
  {Laenge}
  laDezimal := la.grad + la.min / 60 + la.sek / 3600;
  {64}
  bf := brDezimal / rho;
  {65}
  g := 111120.61962 * brDezimal
       -15988.63853 * sin(2*bf)
       +16.72995 * sin(4*bf)
       -0.02178 * sin(6*bf)
       +0.00003 * sin(8*bf);
  {70}
  co := cos(bf);
  {71}
  g2 := e2 * (co * co);
  {72}
  g1 := c / sqrt(1+g2);
  {73}
  t := sin(bf) / cos(bf); {=tan(t)}
  {74}
  dl := laDezimal - sy * 3;
  {77}
  fa := co * dl / rho;
  {78}
  y := g
       + fa * fa * t * g1 / 2
       + fa * fa * fa * fa * t * g1 * (5 - t * t + 9 * g2) / 24;
  {81}
  rm := fa * g1
        + fa * fa * fa * g1 * (1 - t * t + g2) / 6
        + fa * fa * fa * fa * fa * g1 * (5 - 18 * t * t * t * t * t * t) / 120;
  {84}
  x := rm + sy * 1000000 + 500000;
end;

{==============================================================================}
{ Umrechnung von Geographischen Koordinaten in Gauss-Krueger-Koordinaten       }
{ Formel: Grossmann,W., Geodätische Abbildungen, 1964, Seite 153               }
{ Parameter: Rechtswert                            : double                    }
{            Hochwert                              : double                    }
{            geo.Breite (Grad.Min.Sek) in Altgrad  : Twinkel                   }
{            geo.Laenge (Grad.Min.Sek) in Altgrad  : Twinkel                   }
{==============================================================================}
procedure GkGeo(rw,hw:double;var br,la:Twinkel);
const
  {26}
  rho= 180 / pi;
var
  rm,e2,c,bI,bII,bf,co,g2,g1,t,fa,dl,gb,gl,grad,min,sek :extended;
  mKen :integer;
begin
  {25}
  e2 := 0.0067192188;
  {27}
  c := 6398786.849;
  {32}
  mKen := trunc(rw / 1000000);
  {33}
  rm := rw - mKen * 1000000 - 500000;
  {34}
  bI := hw / 10000855.7646;
  {35}
  bII := bI * bI;
  {36}
  bf := 325632.08677 *bI *((((((0.00000562025
        * bII - 0.00004363980)
        * bII + 0.00022976983)
        * bII - 0.00113566119)
        * bII + 0.00424914906)
        * bII - 0.00831729565)
        * bII + 1);
  {43}
  bf := bf / 3600 / rho;
  {44}
  co := cos(bf);
  {45}
  g2 := e2 * (co * co);
  {46}
  g1 := c / sqrt(1 + g2);
  {47}
  t := sin(bf) / cos(bf); {=tan(bf)}
  {51}
  fa := rm / g1;
  {52}
  gb := bf
        - fa * fa * t * (1 + g2) / 2
        + fa * fa * fa * fa * t * (5 + 3 * t * t + 6 * g2 - 6 * g2 * t * t) / 24;
  {55}
  gb := gb * rho;
  {56}
  dl := fa
        - fa * fa * fa * (1 + 2 * t * t + g2) / 6
        + fa * fa * fa * fa * fa * (1 + 28 * t * t + 24 * t * t * t * t) / 120;
  {59}
  gl := dl *rho / co + mKen * 3;
  {in grad.min.sek}
  {Breite}
  grad:=int(gb);
   sek:=60*(gb-grad);
   min:=int(sek);
   sek:=60*(sek-min);
  br.grad := trunc(grad);
  br.min  := trunc(min);
  br.sek  := trunc(sek);
  {Laenge}
  grad:=int(gl);
   sek:=60*(gl-grad);
   min:=int(sek);
   sek:=60*(sek-min);
  la.grad := trunc(grad);
  la.min  := trunc(min);
  la.sek  := trunc(sek);
end;