{$E+,N+}                  {alle Rechnungen im IEEE-Format(1+15+64 Bit) }
program tvsat;            {TurboPascal 5.5, 6 oder 7}

const D     : extended = 12756.27;           {Erddurchmesser am Äquator}
      qR    : extended = 40680606.1;         {Erdradius am Äquator(=6378.135 km)}
                                             {quadriert}
      sath  : extended = 42164.18;           {Höhe der geostationären Umlaufbahn}
                                             {über Erdmittelpunkt in km}
      qsath : extended = 1777818075.0724;    {Höhe quadriert}
      bm    : extended = 0.0174532925199433; {pi/180 = Bogenmaß-Koeffizient}

var   satpos,                                {Satpositionswinkel in Ostrichtung}
      c,                                     {Entfernung des Sat vom QTH in km}
      dl,                                    {Längendifferenz von Sat und QTH}
      qthl,                                  {QTH-Länge in Grad Ost}
      qthb,                                  {QTH-Breite in Grad Nord}
      elevation,                             {...}
      azimut : extended;                     {...}
      key : char;                            {Eingabedummy}

function arccos(wert:extended):extended;
{berechnet den Arcuscosinus im Bogenmaß}
begin
  arccos:=pi/2-arctan(wert/(sqrt(1-wert*wert)+1e-18))
end;

begin
    {Eingabe}
  writeln('TV-SAT Version 1.1 (C) DB1HZ 1990');
  writeln;
  write('Standort des TV-Satelliten in Grad Ost: ');
  readln(satpos);
  writeln;
  write('eigene Länge in Grad Ost: ');
  readln(qthl);
  writeln;
  write('eigene Breite in Grad Nord: ');
  readln(qthb);
  qthb:=qthb*bm;
    {Längendifferenz von Satellit und QTH}
  dl:=(qthl-satpos)*bm;
    {Entfernung QTH-Satellit in km über Seitencosinussatz}
  c:=sqrt(qR+qsath-D*sath*cos(dl)*cos(qthb));
    {Elevationswinkel des Satelliten am QTH vom Erdmittelpunkt in Grad}
    {über Seitencosinussatz}
  elevation:=arccos((qR+c*c-qsath)/(D*c))/bm-90;
    {Azimutwinkel des Satelliten am QTH zur Südrichtung in Grad}
    {über rechtwinkliges sphärisches Dreieck}
  azimut:=arctan(sin(dl)/(cos(dl)*sin(qthb)))/bm;
    {Nördl. und südl. Hemisphäre unterscheiden, da ARCTAN nicht eindeutig}
  if (qthb<0) then azimut:=360*frac(azimut/360+1) else azimut:=azimut+180;
    {Ausgabe}
  writeln;
  writeln('Azimut    = ',azimut:3:3,' Grad');
  writeln;
  writeln('Elevation = ',elevation:3:3,' Grad');
  writeln
end.