PROGRAM PARABOLBAHN (IN,OUT);
CONST
  PI=3.14159266535;K=0.01720209895;BOG=0.01745329252;

TYPE
  DREIREL=ARRAY [1..3] OF REAL;

VAR
  R,A,B,C,X,Y,Z,JD:DREIREL;
  A12,B12,C12,A23,B23,C23,G,F1,L1,M,H1,H2,H3,H,F3,L3,F,L,ETA,THETA,
  S,RHO,R1,R3,R13,RN13,X1,Y1,Z1,X3,Y3,Z3,HST,HST1,X0,Y0,Z0,R0,V31,
  EKL,S1,V1,V3,Q,T,SV,CV,PX,PY,PZ,QX,QY,QZ,H31,H12,H23,IB,KNL,OM,D,
  RA,DEC:REAL;
  I,J1:INTEGER;

FUNCTION ARCSIN (X:REAL):REAL;
BEGIN
  ARCSIN:=ARCTAN (X/SQRT (-X*X+1));
END;

FUNCTION ARCCOS (X:REAL):REAL;
BEGIN
  ARCCOS:=-ARCSIN (X)+PI/2;
END;

FUNCTION TAN (X:REAL):REAL;
BEGIN
  TAN:=SIN(X)/COS(X);
END;

PROCEDUR XYZ (VAR D,M:REAL;J:INTEGER);
VAR
  S,N,YS,EN,D1,T,OM,EX,EL,MA,E,R,V:REAL;
BEGIN
  YS:=J+1900;N:=M+1;IF N<=3 THEN BEGIN YS:=J+1899;N:=N+12 END;
  JD[I]:=TRUNC (365.25*YS)+TRUNC (30.6001*N)+D+1720981.5;
  D1:=TRUNC((J-1)/4)+TRUNC(30.6*M+0.53/SQR (M-1.55)-32.3)+D-0.5;
  IF (J MOD 4=0) AND (M>=3) THEN D1:=D1+1;
  T:=(365*J+D1)/36525;OM:=BOG*(281.220833+1.719175*T+3.61E-4*T*T);
  EX:=0.01675104-4.18E-5*T;EKL:=BOG*(23.452294-0.0130125*T);
  MA:=BOG*(-1.524155-1.5E-4*T*T-0.25590255*J+0.98560027*D1);
  E:=MA;
  REPEAT EN:=E;E:=MA+EX*SIN (EN) UNTIL ABS (EN-E)<1E-7;
  R:=1-EX*COS (E);V:=2*ARCTAN (SQRT ((1+EX)/(1-EX))*TAN (E));
  EL:=V+OM;
  X[I]:=R*COS (EL);Y[I]:=R*SIN (EL)*COS (EKL);Z[I]:=R*SIN (EL)*SIN (EKL);
END;

PROCEDUR LOESUNG (VAR F:REAL);
VAR
  H,AQ1,AQ2,AQ3,P,Q,D:REAL;
BEGIN
  H:=-(A[2]*X[2]+B[2]*Y[2]+C[2]*Z[2]);
  AQ1:=1/3/R[2]*(4*F+5*H);AQ2:=2/3+1/SQR (R[2])*(2*F*H+SQR (G)/3/SQR (H));
  AQ3:=1/3/SQR (R[2])*(2*F+SQR (G/H/R[2])*H);
  P:=-1/3*SQR (AQ1)+AQ2;Q:=2/27*AQ1*AQ1*AQ1-1/3*AQ1*AQ2+AQ3;
  D:=SQR (P/3)*P/3+SQR (Q/2);
  IF D>0 THEN WRITELN("ES GIBT NUR EINE LOESUNG(MIT NACHRECHNUNG)!")
    ELSE BEGIN WRITELN("ES GIBT 3 LOESUNGEN!");HALT;END;
END;

(*---------- Hauptprogramm ----------*)
BEGIN
  FOR I:=1 TO 3 DO BEGIN
    WRITE("DATUM",I:2);READLN;READ(D,M,J1);
  XYZ(D,M,J1);
  R[I]:=SQRT (SQR (X[I])+SQR (Y[I])+SQR (Z[I]));
  WRITE("DEC,RA",I:2);READLN;READ(DEC,RA);DEC:=BOG*DEC;RA:=BOG*15*RA;
  A[I]:=COS (DEC)*COS(RA);B[I]:=COS (DEC)*SIN (RA);C[I]:=SIN (DEC);
                   END;
  A12:=B[1]*C[2]-B[2]*C[1];B12:=C[1]*A[2]-C[2]*A[1];
  C12:=A[1]*B[2]-A[2]*B[1];A23:=B[2]*C[3]-B[3]*C[2];
  B23:=C[2]*A[3]-C[3]*A[2];C23:=A[2]*B[3]-A[3]*B[2];
  G:=SQRT (SQR (X[3]-X[1])+SQR (Y[3]-Y[1])+SQR (Z[3]-Z[1]));
  F1:=-(A[1]*X[1]+B[1]*Y[1]+C[1]*Z[1]);L1:=SQRT (SQR  (R[1])-SQR (F1));
  M:=(JD[3]-JD[2])/(JD[2]-JD[1]);
  M:=M*(A12*X[2]+B12*Y[2]+C12*Z[2])/(A23*X[2]+B23*Y[2]+C23*Z[2]);
  H1:=M*A[3]-A[1];H2:=M*B[3]-B[1];H3:=M*C[3]-C[1];
  H:=SQRT (SQR (H1)+SQR (H2)+SQR (H3));
  F3:=-(A[3]*X[3]+B[3]*Y[3]+C[3]*Z[3])/M;
  L3:=SQRT (SQR (R[3]/M)-SQR (F3));
  F:=(H1*(X[1]-X[3])+H2*(Y[1]-Y[3])+H3*(Z[1]-Z[3]))/SQR (H);
  L:=SQRT (SQR (G/H)-SQR (F));
  IF F<=0 THEN LOESUNG(F) ELSE WRITELN("ES GIBT NUR EINE LOESUNG!");
  RN13:=2;

  (*========== Iteration R13 ==========*)
  REPEAT R13:=RN13;
         ETA:=2*K*(JD[3]-JD[1])/EXP (R13*LN(1.5));
         THETA:=ARCSIN (ETA*1.060660172);
         S:=R13*2*SQRT (2)*SIN (THETA/3)*SQRT (COS (2/3*THETA));
         RHO:=SQRT (SQR (S/H)-SQR (L))-F;
         R1:=SQRT (SQR (RHO+F1)+SQR (L1));
         R3:=M*SQRT (SQR (RHO+F3)+SQR (L3));
         RN13:=R1+R3;
  UNTIL  ABS (RN13-R13)<1E-8;

  (*========== Ableitung der Bahnelemente ==========*)
  X[1]:=RHO*A[1]-X[1];Y[1]:=RHO*B[1]-Y[1];Z[1]:=RHO*C[1]-Z[1];
  RHO:=M*RHO;
  X[3]:=RHO*A[3]-X[3];Y[3]:=RHO*B[3]-Y[3];Z[3]:=RHO*C[3]-Z[3];
  HST:=X[1]*X[3]+Y[1]*Y[3]+Z[1]*Z[3];
  HST1:=HST/SQR (R1);
  X0:=X[3]-HST1*X[1];
  Y0:=Y[3]-HST1*Y[1];
  Z0:=Z[3]-HST1*Z[1];
  R0:=SQRT (SQR (X0)+SQR (Y0)+SQR (Z0));
  V31:=ARCTAN (R0*R1/HST);
  H1:=1/SQRT (R1);
  H2:=H1*COS (V31/2);
  H3:=1/SQRT (R3);
  S1:=(H2-H3)/SIN (V31/2);
  V1:=2*ARCTAN (S1/H1);
  V3:=V1+V31;
  Q:=R1*SQR (COS (V1/2));
  T:=JD[1]-SQRT (2)/K*EXP (Q*LN(3/2))*TAN (V1/2)*(1+1/3*SQRT (TAN (V1/2)));
  CV:=COS (V1);SV:=SIN (V1);

  (*========== Gauss'sche Konstanten ==========*)
  PX:=X[1]/R1*CV-X0/R0*SV;
  PY:=Y[1]/R1*CV-Y0/R0*SV;
  PZ:=Z[1]/R1*CV-Z0/R0*SV;
  QX:=X[1]/R1*SV+X0/R0*CV;
  QY:=Y[1]/R1*SV+Y0/R0*CV;
  QZ:=Z[1]/R1*SV+Z0/R0*CV;
    WRITELN("*=========================================*");
    WRITELN("DIE GAUSS'SCHEN KONSTANTEN DES OBJEKTS:");
    WRITELN("PX =",PX:16:6);
    WRITELN("PY =",PY:16:6);
    WRITELN("PZ =",PZ:16:6);
    WRITELN("QX =",QX:16:6);
    WRITELN("QY =",QY:16:6);
    WRITELN("QZ =",QZ:16:6);
  H31:=PZ*QX-PX*QZ;
  H12:=PX*QY-PY*QX;
  H23:=PY*QZ-PZ*QY;
  IB:=ARCCOS (H12*COS (EKL)-H31*SIN (EKL));
  KNL:=ARCSIN (H23/SIN (IB));
  OM:=ARCSIN (-(QX*COS (KNL)+QY*SIN (KNL)*COS (EKL)+QZ*SIN (KNL)*SIN (EKL)));

  (*========== Ausgabe ==========*)
  WRITELN("*=========================================*");
  WRITELN("DIE BAHNELEMENTE:");
    WRITELN("  T  =",T:20:5," (JD)");
    WRITELN("  Q  =",Q:20:5);
    WRITELN("  I  =",IB/BOG:20:5);
    WRITELN("  KNL=",KNL/BOG:20:5);
    WRITELN("  OMG=",OM/BOG:20:5);
  WRITELN("*==========================================*");
END.