*****************************************************;
* Program Name: Geostationary Satellite Declination *;
* Last Updated: 06/24/20 by David Oesper            *;
*****************************************************;

data decsatm; /* Geostationary Satellite Declination on the Observer's Meridian */
   /**************************************************/
   latobs = 42.96025; /* latitude of the observer    */
   /**************************************************/
   latobs2 = latobs*(constant('pi')/180);
   decsat = atan((6.611/sin(latobs2))-(1/tan(latobs2)))*(180/constant('pi')) - 90;
run;

proc print data=decsatm;
   title "Geostationary Satellite Declination on the Observer's Meridian";
run;

data decsatl; /* Geostationary Satellite Declination Using the Longitude of the Satellite */
   /**************************************************/
   lngsat = -75.2;     /* longitude of the satellite */
   lngobs = -90.14197; /* longitude of the observer  */
   latobs = 42.96025;  /* latitude of the observer   */
   /**************************************************/
   lngsat2 = lngsat*(constant('pi')/180);  
   lngobs2 = lngobs*(constant('pi')/180);
   lngrel2 = lngsat2 - lngobs2;
   latobs2 = latobs*(constant('pi')/180);
   h2 = atan(sin(lngrel2)/(cos(lngrel2) - 0.15126*cos(latobs2)));
   decsat2 = atan(-0.15126*sin(latobs2)*sin(h2)/sin(lngrel2));
   decsat = decsat2*(180/constant('pi'));
run;

proc print data=decsatl;
   title 'Geostationary Satellite Declination Using the Longitude of the Satellite';
run;

data decsath; /* Geostationary Satellite Declination Using the Hour Angle of the Satellite (here, HA = Right Ascension of the Satellite - Local Sidereal Time */
   /**************************************************/
   h = 1.12139;       /* hour angle of the satellite */
   latobs = 42.96025; /* latitude of the observer    */
   /**************************************************/
   h2 = h*(constant('pi')/12);
   latobs2 = latobs*(constant('pi')/180);
   sinl = tan(h2)*(cos(h2)-0.15126*cos(latobs2)); /* starting value */
   do until(sinl-(tan(h2)*(cos(sinl)-0.15126*cos(latobs2))) < 0.000005);
      sinl = tan(h2)*(cos(sinl)-0.15126*cos(latobs2));
   end;
   decsat2 = atan(-0.15126*sin(latobs2)*sin(h2)/sinl);
   decsat = decsat2*(180/constant('pi'));
run;

proc print data=decsath;
   title1 'Geostationary Satellite Declination Using the Hour Angle of the Satellite';
   title2 '(here, Hour Angle = Right Ascension of the Satellite - Local Sidereal Time)';
run;