(*---------------------------------------------------------------------------*)
(*                                 SUNSET                                    *)
(*               Auf- und Untergangszeiten von Sonne und Mond                *)
(*                  Für Lazarus mit Free Pascal 24.9.2018                    *)
(*---------------------------------------------------------------------------*)

PROGRAM sunset2(INPUT,OUTPUT);


  VAR   ABOVE,RISE,SETT                      : BOOLEAN;
        DAY,MONTH,YEAR, I,IOBJ,NZ            : INTEGER;
        LAMBDA,ZONE,PHI,SPHI,CPHI            : REAL;
        TSTART,DATE,HOUR,HH,UTRISE,UTSET     : REAL;
        Y_MINUS,Y_0,Y_PLUS,ZERO1,ZERO2,XE,YE : REAL;
        SINH0                                : ARRAY[1..3] OF REAL;

(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\MATH\SN.PAS         *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\MATH\CS.PAS         *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\MATH\QUAD.PAS       *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\TIME\MJD.PAS        *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\TIME\CALDAT.PAS     *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\TIME\LMST.PAS       *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\SUN\MINISUN.PAS     *)
(*$I C:\USERS\SM64\DESKTOP\ASTROPAS\MOON\MINIMOON.PAS   *)


(*---------------------------------------------------------------------------*)
(* SIN_ALT: sin(Horizonthoehe)                                               *)
(*         IOBJ:  1=Mond, 2=Sonne                                            *)
(*---------------------------------------------------------------------------*)
FUNCTION SIN_ALT(IOBJ:INTEGER;MJD0,HOUR,LAMBDA,CPHI,SPHI:REAL): REAL;
  VAR MJD,T,RA,DEC,TAU: REAL;
  BEGIN
    MJD := MJD0 + HOUR/24.0;
    T   := (MJD-51544.5)/36525.0;
    IF (IOBJ=1)
      THEN  MINI_MOON(T,RA,DEC)
      ELSE  MINI_SUN (T,RA,DEC);
    TAU := 15.0 * (LMST(MJD,LAMBDA) - RA);
    SIN_ALT  := SPHI*SN(DEC) + CPHI*CS(DEC)*CS(TAU);
  END;
(*---------------------------------------------------------------------------*)
PROCEDURE WHM(UT:REAL);
  VAR H,M:INTEGER;
  BEGIN
    UT := TRUNC(UT*60.0+0.5)/60.0; (* Rundung auf 1 min *)
    H:=TRUNC(UT); M:=TRUNC(60.0*(UT-H)+0.5); WRITE(H:5,':',M:2,' ');
  END;
(*---------------------------------------------------------------------------*)
PROCEDURE GETINP(VAR DATE,LAMBDA,PHI,ZONE: REAL);
  VAR D,M,Y: INTEGER;
  BEGIN
    WRITELN;
    WRITELN('      SUNSET: Auf-/Untergangszeiten von Sonne und Mond');
    WRITELN('                      Aktuelle Version 23.09.2018     ');
    WRITELN;
    WRITELN;
    WRITE  (' Startdatum (TT MM JJJJ)                ... '); READLN(D,M,Y);
    WRITELN;
    WRITE  (' Beobachtungsort:  Laenge (oestl. neg.) ... '); READLN(LAMBDA);
    WRITE  ('                   Breite               ... '); READLN(PHI);
    WRITE  ('                   Zonenzeit-UT (in h)  ... '); READLN(ZONE);
    WRITELN;
    WRITELN;
    WRITELN('    Datum           Mond             Sonne            Daemmerung');
    WRITELN;
    WRITELN('               Auf-/Untergang    Auf-/Untergang      Ende/Anfang');
    WRITELN;
    ZONE := ZONE /24.0;
    DATE := MJD(D,M,Y,0)-ZONE;
  END;
(*---------------------------------------------------------------------------*)

BEGIN (* SUNSET *)

  SINH0[1] := SN ( +8.0/60.0); (* Mondaufgang      bei h= +8'     *)
  SINH0[2] := SN (-50.0/60.0); (* Sonnenaufgang    bei h=-50'     *)
  SINH0[3] := SN (   -12.0  ); (* naut. Daemmerung bei h=-12 Grad *)

  GETINP (TSTART,LAMBDA,PHI,ZONE);

  SPHI := SN(PHI);  CPHI := CS(PHI);

  FOR I:=0 TO 9 DO  (* Schleife ueber die Tage des Suchzeitraumes *)

     BEGIN

       DATE := TSTART + I;
       CALDAT(DATE+ZONE,DAY,MONTH,YEAR,HH);
       WRITE(DAY:3,'.',MONTH:2,'.',YEAR:4,' ');  (* aktuellen Tag anzeigen *)

       FOR IOBJ := 1 TO 3 DO

         BEGIN

           HOUR := 1.0;

           Y_MINUS := SIN_ALT(IOBJ,DATE,HOUR-1.0,LAMBDA,CPHI,SPHI)-SINH0[IOBJ];

           ABOVE := (Y_MINUS>0.0); RISE := FALSE; SETT := FALSE;

           (* Schleife ueber die Suchintervalle von [0h-2h] bis [22h-24h]  *)
           REPEAT

             Y_0   :=SIN_ALT(IOBJ,DATE,HOUR,LAMBDA,CPHI,SPHI)-SINH0[IOBJ];
             Y_PLUS:=SIN_ALT(IOBJ,DATE,HOUR+1.0,LAMBDA,CPHI,SPHI)-SINH0[IOBJ];

             (* Parabel durch die drei Werte Y_MINUS,Y_0,Y_PLUS legen *)
             QUAD ( Y_MINUS,Y_0,Y_PLUS, XE,YE, ZERO1,ZERO2, NZ );

             CASE (NZ) OF
               0: ;
               1: IF (Y_MINUS<0.0)
                    THEN BEGIN UTRISE:=HOUR+ZERO1;  RISE:=TRUE; END
                    ELSE BEGIN UTSET :=HOUR+ZERO1;  SETT:=TRUE; END;
               2: BEGIN
                    IF (YE<0.0)
                      THEN BEGIN UTRISE:=HOUR+ZERO2; UTSET:=HOUR+ZERO1; END
                      ELSE BEGIN UTRISE:=HOUR+ZERO1; UTSET:=HOUR+ZERO2; END;
                    RISE:=TRUE; SETT:=TRUE;
                  END;
              END;

              Y_MINUS := Y_PLUS;     (* Tabelle weitergeben fuer *)
              HOUR := HOUR + 2.0;    (* naechstes Suchintervall  *)

           UNTIL ( (HOUR=25.0) OR (RISE AND SETT) );

           IF (RISE OR SETT) (* Ausgabe *)
             THEN
               BEGIN
                 IF RISE THEN WHM(UTRISE) ELSE WRITE('----- ':9);
                 IF SETT THEN WHM(UTSET)  ELSE WRITE('----- ':9);
               END
             ELSE
               BEGIN
                 IF ABOVE
                   THEN CASE IOBJ OF
                          1,2: WRITE ('   immer sichtbar ');
                          3:   WRITE ('     immer hell   ');
                         END
                   ELSE CASE IOBJ OF
                          1,2: WRITE ('  immer unsichtbar');
                          3:   WRITE ('    immer dunkel  ');
                         END;
               END;

         END;

       WRITELN;

     END; (* Ende der Schleife ueber die Tage des Suchzeitraumes *)


  WRITELN;  WRITE  ('  alle Zeiten in Zonenzeit ( = UT ');
  IF ZONE>=0 THEN WRITE('+'); WRITELN(ZONE*24.0:5:1,'h )');

END.