program moon;

(*   Program to calculate azimuth and elevation of the moon.

   May 1974

   Translated from FORTRAN to BASIC *Mitchell (WB4BWK)

   Re-written into PASCAL 30-May-80 *Mitchell (WB4BWK)

   This program is based upon Lance Collister's work (WA1JXN/WA3GPL) who
   acknowledged the vital help of Thomas Ake of the Warner-Swasey Observatory.

   WB4BWK < K3IXD < W6PO < WA1JXN/WA3GPL.

   *)

CONST
  pi = 3.14159265;

TYPE
  monthtypes=(zero,jan,feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec,nextyear);

  char3 = PACKED ARRAY[1..3] OF char;

VAR
  beginninghour, d6, d7, d8, d9, day, daycount, endinghour, h8, nrdays,
  o1, o2, o3, o4, time, timeinc, year, z : integer;

  done : boolean;

  month : monthtypes;

  monthdays : ARRAY[monthtypes] OF integer;

  monthname : ARRAY[monthtypes] OF char3;

  call : PACKED ARRAY[1..10] OF char;

  ans, monthin : char3;

  a, a1, a2, c1, c2, d1, degtoradians, e, e2, e3, elmax, factor1,
  f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, g, h, i1,
  j, j1, j2, k0, k1, k2, k2k2, k3, k4, k4k4, k5, k6, l, l1, l7, l8, latmin,
  latitude, longmin, longtitude, r1, radianstodeg, t, t1, t5, twopi, x : real;


  FUNCTION leapyear(year:integer) : boolean;
   BEGIN
     leapyear := (year MOD 4 = 0)
		 AND ((year MOD 100 <> 0) OR (year MOD 400 = 0))
   END;

  FUNCTION atangent(y,x:real) : real;
  VAR
    arctangent : real;
   BEGIN
     arctangent := arctan(y/x);
     IF x < 0.0
      THEN arctangent := arctangent + pi
      ELSE
      IF y < 0.0
       THEN arctangent := arctangent + twopi;
     atangent := arctangent
   END;

  PROCEDURE newday;
   BEGIN
     writeln;
     writeln;
     time := beginninghour;
     writeln('Position of the moon on ',
	     day:2, ' ', monthname[month], ' ', year:4);
     writeln;
     writeln(' GMT             Az         El           GHA         Dec');
     writeln;
     i1 := 2.0;
     IF month < mar
      THEN
       BEGIN
	 IF (year-1853) DIV 4 < 11
	  THEN c1 := 0.0
	  ELSE c1 := -1.0;
	 j1 := 365.0*(year-1853)
		+ day + 30.0*(ord(month)+9) + trunc((ord(month)+10)/2.0);
	 j2 := trunc((year-1853)/4.0) + 1.0 + c1
       END
      ELSE
       BEGIN
	 IF (year-1852) DIV 4 < 11
	  THEN c1 := 0.0
	  ELSE c1 := -1.0;
	 IF (month = sep) OR (month = nov)
	  THEN c2 := 1.0
	  ELSE c2 := 0.0;
	 j1 := 365.0*(year-1852)
		+ day + 30.0*(ord(month)-3) + trunc((ord(month)-2)/2.0);
	 j2 := trunc((year-1852)/4.0) + c1 + c2
       END;

      (*      Julian date - 2397547.5 for 0 hr GMT *)

     j  := j1 + j2;
     t1 := j - 17472.5;
   END;
  (*NEWday*)

  PROCEDURE timecheck;
   BEGIN
     d9 := time DIV 100 * 60 + time MOD 100;
     d7 := d9 - d6;
     d8 := d7 - timeinc;
   END;
  (*TIMEcheck*)

  PROCEDURE init;
   BEGIN

     done := false;
     daycount := 1;

     monthname[jan] := 'Jan';
     monthname[feb] := 'Feb';
     monthname[mar] := 'Mar';
     monthname[apr] := 'Apr';
     monthname[may] := 'May';
     monthname[jun] := 'Jun';
     monthname[jul] := 'Jul';
     monthname[aug] := 'Aug';
     monthname[sep] := 'Sep';
     monthname[oct] := 'Oct';
     monthname[nov] := 'Nov';
     monthname[dec] := 'Dec';

     monthdays[jan] := 31;
     monthdays[mar] := 31;
     monthdays[apr] := 30;
     monthdays[may] := 31;
     monthdays[jun] := 30;
     monthdays[jul] := 31;
     monthdays[aug] := 31;
     monthdays[sep] := 30;
     monthdays[oct] := 31;
     monthdays[nov] := 30;
     monthdays[dec] := 31;

     writeln;
     writeln;
     writeln('Moon  V7.0 (Pascal)');
     writeln;

     twopi := 2.0 * pi;
     radianstodeg := 180.0/pi;
     degtoradians := pi/180.0;
     f1  := 0.658  * degtoradians;
     f2  := 6.289  * degtoradians;
     f3  := 1.274  * degtoradians;
     f4  := 0.186  * degtoradians;
     f5  := 0.214  * degtoradians;
     f6  := 0.114  * degtoradians;
     f7  := 0.059  * degtoradians;
     f8  := 0.057  * degtoradians;
     f9  := 0.6593 * degtoradians;
     f10 := 6.2303 * degtoradians;
     f11 := 1.2720 * degtoradians;
     f12 := 5.144  * degtoradians;
     f13 := 0.146  * degtoradians;
     factor1 := 24.0 * 1.002738;
     elmax := 100.0;

     call := 'WB4BWK    ';
     latitude := 34.0;
     latmin := 19.7833;
     longtitude := 82.0;
     longmin := 22.6333;
     timeinc := 30;
     beginninghour := 0;
     endinghour := 2400;


     write('Use preset values?  ');
     read(ans);
     IF (ans[1]<>'Y') AND (ans[1]<>'y')
      THEN
       BEGIN
	 writeln;
	 write('Do you want preset moon rise and moon set only?  ');
	 readln;
	 read(ans);
	 IF (ans[1]='Y') OR (ans[1]='y')
	  THEN
	   BEGIN
	     timeinc := 1;
	     elmax := 1.0
	   END
	  ELSE
	   BEGIN
	     writeln;
	     write('Your call?  ');
	     readln;
	     read(call);
	     writeln;
	     write('Your latitude and longitude (degrees, minutes)?  ');
	     read(latitude, latmin, longtitude, longmin);
	     writeln;
	     write('Time increments?  ');
	     read(timeinc);
	     writeln;
	     write('Beginning hour and ending hour each day?  ');
	     read(beginninghour, endinghour);
	     writeln;
	     write('Limit elevation?  ');
	     readln;
	     read(ans);
	     IF (ans[1]='Y') OR (ans[1]='y')
	      THEN
	       BEGIN
		 writeln;
		 write('What maximum elevation?  ');
		 read(elmax)
	       END;
	   END;
       END;

     writeln;

      REPEAT
       write('Beginning month, day, year (21-Jan-80 = Jan 21 80)?  ');
       readln;
       read(monthin, day, year);
       IF ord(monthin[1]) > ord('S')
	THEN monthin[1] := chr( ord(monthin[1]) - 32 );
       IF ord(monthin[2]) < ord('a')
	THEN monthin[2] := chr( ord(monthin[2]) + 32 );
       IF ord(monthin[3]) < ord('a')
	THEN monthin[3] := chr( ord(monthin[3]) + 32 );
       month := jan;
       WHILE (month < dec) AND (monthname[month] <> monthin) DO
		month := succ(month)
      UNTIL monthname[month]=monthin;
     IF year < 100
      THEN year := year + 1900;
     IF leapyear(year)
      THEN monthdays[feb] := 29
      ELSE monthdays[feb] := 28;
     writeln;
     write('Calculations for how many days?  ');
     read(nrdays);
     writeln;
   END; (*Init*)

 BEGIN (*Main*)
   init;
   writeln;
   writeln;
   writeln('Earth-Moon-Earth data for ', call);
   writeln;
   writeln('Location:  ', latitude:3:0, ' degrees ',
	   latmin:6:4, ' minutes latitude');
   writeln('           ', longtitude:3:0, ' degrees ',
	   longmin:6:4, ' minutes longtitude');
   latitude   := (latitude   + latmin/60.0) * degtoradians;
   longtitude := (longtitude + longmin/60.0) * degtoradians;
   elmax := elmax * degtoradians;
   d6 := endinghour DIV 100 * 60 + endinghour MOD 100;
   newday;
   timecheck;

    REPEAT

     (*      Calculations of lat and long of the moon: *)

     t := (time MOD 100)/1440.0 + (time DIV 100)/24.0;
     t5 := t1 + t;

     k0 := 0.751213   + 0.036601102 * t5;
     k1 := (k0 - trunc(k0)) * twopi;

     k0 := 0.822513  + 0.0362916457 * t5;
     k2 := (k0 - trunc(k0)) * twopi;

     k0 := 0.995766 + 0.00273777852 * t5;
     k3 := (k0 - trunc(k0)) * twopi;

     k0 := 0.974271  + 0.0338631922 * t5;
     k4 := (k0 - trunc(k0)) * twopi;

     k0 := 0.0312525 + 0.0367481957 * t5;
     k5 := (k0 - trunc(k0)) * twopi;

     k2k2 := 2.0 * k2;
     k4k4 := 2.0 * k4;

     l8  :=  k1 + f1 * sin(k4k4) + f2 * sin(k2) - f3 * sin(k2 - k4k4)
		- f4 * sin(k3) + f5 * sin(k2k2) - f6 * sin(2.0 * k5)
		- f7 * sin(k2k2 - k4k4) - f8 * sin(k2 + k3 - k4k4);

     k6  :=  k5 + f9 * sin(k4k4) + f10 * sin(k2) - f11 * sin(k2 - k4k4);

     l7  := f12 * sin(k6) - f13 * sin(k5 - k4k4);

     (*      Calculations of right ascention (A = R1) and declination (D1) *)

     d1 := cos(l7) * sin(l8) * 0.397821 + sin(l7) * 0.917463;
     d1 := arctan(d1/(sqrt(1.0 - d1 * d1)));
     a2 := cos(l7) * cos(l8)/cos(d1);
     a1 := (cos(l7) * sin(l8) * 0.917463 - sin(l7) * 0.397821)/cos(d1);
     a := atangent(a1,a2);

     r1 := a;
     l1 := 0.065709822 * t1;
     l  := t * factor1 + 6.646055 + (l1-trunc(l1/24.0) * 24.0);
     l  := (l - trunc(l/24.0) * 24.0);

     (*      Calculations of greenwich hour angle, G, from sideral time. *)

     g := (l/24.0) * twopi - r1;
     IF g < 0.0
      THEN g := g + twopi
      ELSE
      IF g >= twopi
       THEN g := g - twopi;

     (*      Calculations of local hour angle, H, from GHA. *)

     h := longtitude - g;

     (*      Calculations of elevation, E, of object sin of el *)

     e3 := cos(latitude) * cos(h) * cos(d1) + sin(d1) * sin(latitude);
     e2 := sqrt(1.0 - e3 * e3);
     e := arctan(e3/e2);
     IF (e >= 0.0) AND ( e <= elmax)
      THEN
       BEGIN

	 (*      Calculation of azimuth, A, of object *)

	 a2 := sin(d1)/(cos(latitude) * cos(e));
	 a2 := a2-(sin(latitude)/cos(latitude)) * (sin(e)/cos(e));
	 a1 := sin(latitude) * sin(d1) + cos(latitude) * cos(d1) * cos(h);
	 a1 := (sin(h) * cos(d1))/(sqrt(1 - a1 * a1));
	 a  := atangent(a1,a2);

	 IF t-i1 > timeinc*2/1440.0 THEN writeln;

	 a  :=  a * radianstodeg;
	 e  :=  e * radianstodeg;
	 g  :=  g * radianstodeg;
	 d1 := d1 * radianstodeg;

	 h8 := time;
	 o1 := h8 DIV 1000;
	 h8 := h8 MOD 1000;
	 o2 := h8 DIV 100;
	 h8 := h8 MOD 100;
	 o3 := h8 DIV 10;
	 o4 := h8 MOD 10;
	 write(chr(o1+48), chr(o2+48), ':', chr(o3+48), chr(o4+48), '  ');
	 writeln(a:12:1, e:11:1, g:14:1, d1:12:1);
	 i1 := t
       END;

     time := time + timeinc;
     z := (time MOD 100) DIV 60;
     IF z > 0
      THEN time := (time DIV 100)*100 + z * 100;
     timecheck;
     IF (d7 > 0) AND (d8 >= 0)
      THEN
       BEGIN
	 day := day + 1;
	 daycount := daycount + 1;
	 IF daycount > nrdays
	  THEN done := true
	  ELSE
	   BEGIN
	     IF day>monthdays[month]
	      THEN
	       BEGIN
		 month := succ(month);
		 day := 1
	       END;
	     IF month = nextyear
	      THEN
	       BEGIN
		 month := jan;
		 year := year + 1;
		 IF leapyear(year)
		  THEN
		  monthdays[feb] := 29
		  ELSE
		  monthdays[feb] := 28
	       END;
	     newday;
	   END
       END
      ELSE
      IF d7 > 0
       THEN time := endinghour;
    UNTIL done
 END.