program T15;
var
  LatA,LonA: real;
  LatB,LonB: real;
  LatC,LonC: real;
  LatD,LonD: real;

  A1,A2,A3: real;
  B1,B2,B3: real;
  C1,C2,C3: real;
  D1,D2,D3: real;

  e1,e2,e3: real;
  f1,f2: real;
  g1,g2,g3: real;

  v1,v2,v3: real;
  M11,M12,M13: real;
  M21,M22,M23: real;
  M31,M32,M33: real;

  n: real;
  r: real;
  w: real;
  CW,SW: real;

  P1,P2,P3: real;
  Q1,Q2,Q3: real;

  A1P,A2P,A3P,AW: real;
  B1P,B2P,B3P,BW: real;
  P1P,P2P,P3P,PW: real;
  Q1P,Q2P,Q3P,QW: real;

  A1PP,A2PP,A3PP: real;
  B1PP,B2PP,B3PP: real;
  P1PP,P2PP,P3PP: real;
  Q1PP,Q2PP,Q3PP: real;

  E1T,E2T,E3T: real;
  E1TT,E2TT,E3TT: real;

  LatP,LonP,HP: real;
  LatQ,LonQ,HQ: real;

  Inside: integer;

begin
   r := 6371000.8;

   LatA:=49;
   LonA:=8;
   LatB:=49;
   LonB:=10;
   LatC:=50;
   LonC:=9;
   LatD:=48;
   LonD:=9;

   Writeln('A: ',WGS84(LatA,LonA));
   Writeln('B: ',WGS84(LatB,LonB));
   Writeln('C: ',WGS84(LatC,LonC));
   Writeln('D: ',WGS84(LatD,LonD));

   SphereToXYZ(LatA,LonA,R,A1,A2,A3);
   SphereToXYZ(LatB,LonB,R,B1,B2,B3);
   SphereToXYZ(LatC,LonC,R,C1,C2,C3);
   SphereToXYZ(LatD,LonD,R,D1,D2,D3);

   e1:=a2*b3-a3*b2;
   e2:=a3*b1-a1*b3;
   e3:=a1*b2-a2*b1;

   f1:=e1*c1+e2*c2+e3*c3;
   f2:=e1*d1+e2*d2+e3*d3;

   n:=-(f2/f1);

   g1:=n*c1+d1;
   g2:=n*c2+d2;
   g3:=n*c3+d3;

   n := sqrt(g1*g1+g2*g2+g3*g3);
   g1:=g1/n;
   g2:=g2/n;
   g3:=g3/n;
   
   P1:=r*g1;
   P2:=r*g2;
   P3:=r*g3;

   Q1:=-P1;
   Q2:=-P2;
   Q3:=-P3;

(* auf XY abbilden *)
   n := sqrt(e1*e1+e2*e2+e3*e3);
   e1:=e1/n;
   e2:=e2/n;
   e3:=e3/n;

(* um z drehen *)
   w:=-arctan2(E2,E1);
   cw:=cos(w);
   sw:=sin(w);

   M11:=cw; M12:=-sw; M13:=0;
   M21:=sw; M22:=cw;  M23:=0;
   M31:=0;  M32:=0;   M33:=1;

   E1T:= M11*E1 + M12*E2 + M13* E3;
   E2T:= M21*E1 + M22*E2 + M23* E3;
   E3T:= M31*E1 + M32*E2 + M33* E3;

   A1P:= M11*A1 + M12*A2 + M13* A3;
   A2P:= M21*A1 + M22*A2 + M23* A3;
   A3P:= M31*A1 + M32*A2 + M33* A3;

   B1P:= M11*B1 + M12*B2 + M13* B3;
   B2P:= M21*B1 + M22*B2 + M23* B3;
   B3P:= M31*B1 + M32*B2 + M33* B3;

   P1P:= M11*P1 + M12*P2 + M13* P3;
   P2P:= M21*P1 + M22*P2 + M23* P3;
   P3P:= M31*P1 + M32*P2 + M33* P3;

   Q1P:= M11*Q1 + M12*Q2 + M13* Q3;
   Q2P:= M21*Q1 + M22*Q2 + M23* Q3;
   Q3P:= M31*Q1 + M32*Q2 + M33* Q3;

(* drehen um y  *)
   w:=arctan2(E3T,E1T)-pi/2;
   cw:=cos(w);
   sw:=sin(w);

   M11:=cw;  M12:=0;  M13:=sw;
   M21:=0;   M22:=1;  M23:=0;
   M31:=-sw; M32:=0;  M33:=cw;

   E1TT:= M11*E1T + M12*E2T + M13* E3T;
   E2TT:= M21*E1T + M22*E2T + M23* E3T;
   E3TT:= M31*E1T + M32*E2T + M33* E3T;

   A1PP:= M11*A1P + M12*A2P + M13* A3P;
   A2PP:= M21*A1P + M22*A2P + M23* A3P;
   A3PP:= M31*A1P + M32*A2P + M33* A3P;

   B1PP:= M11*B1P + M12*B2P + M13* B3P;
   B2PP:= M21*B1P + M22*B2P + M23* B3P;
   B3PP:= M31*B1P + M32*B2P + M33* B3P;

   P1PP:= M11*P1P + M12*P2P + M13* P3P;
   P2PP:= M21*P1P + M22*P2P + M23* P3P;
   P3PP:= M31*P1P + M32*P2P + M33* P3P;

   Q1PP:= M11*Q1P + M12*Q2P + M13* Q3P;
   Q2PP:= M21*Q1P + M22*Q2P + M23* Q3P;
   Q3PP:= M31*Q1P + M32*Q2P + M33* Q3P;

   writeln;
   writeln('E: ',E1 :12:6,E2 :12:6,E3 :12:6);
   writeln;
   writeln('A: ',A1:12:2,A2:12:2,A3:12:2);
   writeln('E: ',B1:12:2,B2:12:2,B3:12:2);
   writeln('P: ',P1:12:2,P2:12:2,P3:12:2);
   writeln('Q: ',Q1:12:2,Q2:12:2,Q3:12:2);
   writeln;
   writeln('E: ',E1T:12:6,E2T:12:6,E3T:12:6);
   writeln;
   writeln('A: ',A1P:12:2,A2P:12:2,A3P:12:2);
   writeln('B: ',B1P:12:2,B2P:12:2,B3P:12:2);
   writeln('P: ',P1P:12:2,P2P:12:2,P3P:12:2);
   writeln('Q: ',Q1P:12:2,Q2P:12:2,Q3P:12:2);

   writeln;
   writeln('E: ',E1TT:12:6,E2TT:12:6,E3TT:12:6);
   writeln;
   writeln('A: ',A1PP:12:2,A2PP:12:2,A3PP:12:2);
   writeln('B: ',B1PP:12:2,B2PP:12:2,B3PP:12:2);
   writeln('P: ',P1PP:12:2,P2PP:12:2,P3PP:12:2);
   writeln('Q: ',Q1PP:12:2,Q2PP:12:2,Q3PP:12:2);

   AW:=arctan2(A2PP,A1PP);
   BW:=arctan2(B2PP,B1PP);
   PW:=arctan2(P2PP,P1PP);
   QW:=arctan2(Q2PP,Q1PP);

   if AW>BW then begin
      BW:=BW+2*PI
   end;
   if AW>PW then begin
      PW:=PW+2*PI
   end;
   if AW>QW then begin
      QW:=QW+2*PI
   end;

   writeln;
   writeln('A: ',AW/Pi*180:7:2);
   writeln('B: ',BW/Pi*180:7:2);
   writeln('P: ',PW/Pi*180:7:2);
   writeln('Q: ',QW/Pi*180:7:2);
   writeln;

   if (AW<=PW) and (PW<=BW) then begin
      writeln('P liegt zwischen A und B')
   end;
   if (AW<=QW) and (QW<=BW) then begin
      writeln('Q liegt zwischen A und B')
   end;

   XYZtoSphere(P1,P2,P3,R,LatP,LonP);
   XYZtoSphere(Q1,Q2,Q3,R,LatQ,LonQ);

   writeln;
   Writeln('P: ',WGS84(LatP,LonP));
   Writeln('Q: ',WGS84(LatQ,LonQ));

end.