/**************************************************/
/* Program written by David Oesper on 16 Jan 2019 */
/**************************************************/
proc fcmp outlib=work.funcs.user;
   function deg2rad(deg);
      rad = deg*(constant('pi')/180);
      return(rad);
   endsub;
   function dms2rad(dms $);
      rad = 2*constant('pi')*((input(scan(dms,1,' '),3.) + input(scan(dms,2,' '),2.)/60 + input(scan(dms,3,' '),5.2)/3600)/360);
      return(rad);
   endsub;
   function hms2rad(hms $);
      rad = 2*constant('pi')*((input(scan(hms,1,' '),2.) + input(scan(hms,2,' '),2.)/60 + input(scan(hms,3,' '),5.2)/3600)/24);
      return(rad);
   endsub;
   function rad2deg(rad);
      deg = rad*(180/constant('pi'));
      return(deg);
   endsub;
run;

options cmplib=(work.funcs);

data evening_lq (drop=lst0);
   length event $42;
   event = 'Evening End of Astronomical Twilight';
   input @1 date  $6.
         @8 lst0 $11.
		 ;
   lst = hms2rad(lst0);
datalines;
Jan 27  3 12 08.12
Feb 26  5 46  1.42
Mar 27  8 17 50.60
Apr 26 11 02 40.58
May 26 13 52 58.16
Jun 25 16 16 59.10
Jul 24 17 42 40.95
Aug 23 18 42 48.89
Sep 21 19 37 22.52
Oct 21 20 43 21.52
Nov 19 22 07 11.88
Dec 18 23 58 59.76
;
run;

data morning_nm (drop=lst0);
   length event $42;
   event = 'Morning Beginning of Astronomical Twilight';
   input @1 date  $6.
         @8 lst0 $11.
		 ;
   lst = hms2rad(lst0);
datalines;
Jan  5 12 48 27.39
Feb  4 14 32 24.44
Mar  6 15 49 24.60
Apr  5 16 50 19.23
May  4 17 45 24.53
Jun  3 18 56 20.74
Jul  2 20 49 16.61
Jul 31 23 25 55.82
Aug 30  2 13 20.15
Sep 28  4 46 29.73
Oct 27  7 14 35.02
Nov 26  9 45 32.58
Dec 25 12 01 34.75
;
run;

data zodiacal;
   set evening_lq morning_nm;
   lat = dms2rad('+42 57 36.9');
   obl = deg2rad(23.4368); /* 2019.0 */
   ecliptic_horizon_angle = rad2deg(arcos(cos(obl)*sin(lat) - sin(obl)*cos(lat)*sin(lst)));
   if ecliptic_horizon_angle >= 60 then flag = 'BEST';
run;

proc print data=zodiacal;
   title 'Ecliptic Horizon Angle';
run;

proc fcmp outlib=work.funcs.user;
   deletefunc deg2rad;
   deletefunc dms2rad;
   deletefunc hms2rad;
   deletefunc rad2deg;
run;