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 15  2 24 47.00
Feb 13  4 42 13.00 
Mar 15  7 03 14.00
Apr 13  9 21 38.00
May 12 11 44 39.00
Jun 10 14 03 35.00
Jul 10 16 03 21.00
Aug  8 17 31 29.00
Sep  6 18 45 25.00
Oct  6 20 01 39.00
Nov  5 21 29 59.00
Dec  5 23 20 02.00
;
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 21 13 36 52.00
Feb 20 15 16 26.00
Mar 21 16 35 50.00
Apr 20 17 51 42.00
May 19 19 12 01.00
Jun 18 20 58 02.00
Jul 17 23 07 59.00
Aug 16  1 34 35.00
Sep 15  3 57 54.00
Oct 14  6 12 26.00
Nov 13  8 32 21.00
Dec 12 10 47 58.00
;
run;

data zodiacal;
   set evening_lq morning_nm;
   lat = dms2rad('+32 16 10.1');
   obl = deg2rad(23.4363); /* 2023.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;