
function[sbsllat, sbsllon] = subsol(iyr, iday, ihr, imn, sec) 

% Find subsolar geographic latitude and longitude given the
% date and time (Universal Time).
%
%      This is based on formulas in Astronomical Almanac for the
%      year 1996, p.  C24. (U.S.  Government Printing Office,
%      1994).  According to the Almanac, results are good to at
%      least 0.01 degree latitude and 0.025 degree longitude
%      between years 1950 and 2050.  Accuracy for other years has
%      not been tested although the algorithm has been designed to
%      accept input dates from 1601 to 2100.  Every day is assumed
%      to have exactly 86400 seconds; thus leap seconds that
%      sometimes occur on June 30 and December 31 are ignored:
%      their effect is below the accuracy threshold of the algorithm.
%
%      961026 A. D. Richmond, NCAR
%
%      INPUTS:
%         IYR  = Year (e.g., 1994). IYR must be in the range: 1601 to 2100.
%         IDAY = Day number of year (e.g., IDAY = 32 for Feb 1)
%         IHR  = Hour of day    (e.g., 13 for 13:49)
%         IMN  = Minute of hour (e.g., 49 for 13:49)
%         SEC  = Second and fraction after the hour/minute.
%          Note:  While IYR is bounds tested, there is no constraint
%                 placed on values: IDAY,IHR,IMN,SEC; e.g., IHR=25 is
%                 properly interpreted.
%       RETURNS:
%         SBSLLAT = geographic latitude of subsolar point (degrees)
%         SBSLLON = geographic longitude of subsolar point (degrees,between -180 and +180)

% Constants:
	D2R = 0.0174532925199432957692369076847; 
        R2D = 57.2957795130823208767981548147;
	MSGUN = 6;

% Number of years from 2000 to IYR (negative if IYR < 2000):
      YR = iyr - 2000;

% NLEAP (final) = number of leap days from (2000 January 1) to (IYR January 1)
%                 (negative if IYR is before 1997)
      NLEAP = (iyr-1601)/4;
      NLEAP = NLEAP - 99;
      if iyr <= 1900 
	if iyr <= 1600 
	 fprintf('SUBSOLR INVALID BEFORE 1601: INPUT YEAR = %i\n', iyr);       
	else 
	NCENT = (iyr-1601)/100;
	NCENT = 3 - NCENT;
	NLEAP = NLEAP + NCENT;
	end 
      else if iyr >= 2101
	fprintf('SUBSOLR INVALID AFTER 2100:  INPUT YEAR = %i\n',iyr);
      else
      end
      end

% L0 = Mean longitude of Sun at 12 UT on January 1 of IYR:
%     L0 = 280.461 + .9856474*(365*(YR-NLEAP) + 366*NLEAP) 
%	   - (ARBITRARY INTEGER)*360.
%        = 280.461 + .9856474*(365*(YR-4*NLEAP) + (366+365*3)*NLEAP) 
%	   - (ARBITRARY INTEGER)*360.
%        = (280.461 - 360.) + (.9856474*365 - 360.)*(YR-4*NLEAP) 
%	   + (.9856474*(366+365*3) - 4*360.)*NLEAP,
%  where ARBITRARY INTEGER = YR+1.  This gives:
      L0 = -79.549 + (-.238699*(YR-4*NLEAP) + 3.08514E-2*NLEAP);

% G0 = Mean anomaly at 12 UT on January 1 of IYR:
%     G0 = 357.528 + .9856003*(365*(YR-NLEAP) + 366*NLEAP) 
%	   - (ARBITRARY INTEGER)*360.
%        = 357.528 + .9856003*(365*(YR-4*NLEAP) + (366+365*3)*NLEAP) 
%	   - (ARBITRARY INTEGER)*360.
%        = (357.528 - 360.) + (.9856003*365 - 360.)*(YR-4*NLEAP) 
%	   + (.9856003*(366+365*3) - 4*360.)*NLEAP,
%  where ARBITRARY INTEGER = YR+1.  This gives:
      G0 = -2.472 + (-.2558905*(YR-4*NLEAP) - 3.79617E-2*NLEAP);

% Universal time in seconds:
      UT = double(ihr*3600 + imn*60) + sec;

% Days (including fraction) since 12 UT on January 1 of IYR:
      DF = (UT/86400. - 1.5) + iday;

% Addition to Mean longitude of Sun since January 1 of IYR:
      LF = .9856474*DF;

% Addition to Mean anomaly since January 1 of IYR:
      GF = .9856003*DF;

% Mean longitude of Sun:
      L = L0 + LF;

% Mean anomaly:
      G = G0 + GF;
      GRAD = G*D2R;

% Ecliptic longitude:
      LAMBDA = L + 1.915*sin(GRAD) + .020*sin(2.*GRAD);
      LAMRAD = LAMBDA*D2R;
      SINLAM = sin(LAMRAD);

% Days (including fraction) since 12 UT on January 1 of 2000:
      N = DF + double(365*YR + NLEAP);

% Obliquity of ecliptic:
      EPSILON = 23.439 - 4.e-7*N;
      EPSRAD = EPSILON*D2R;

% Right ascension:
      ALPHA = atan2(cos(EPSRAD)*SINLAM,cos(LAMRAD))*R2D;

% Declination:
      DELTA = asin(sin(EPSRAD)*SINLAM)*R2D;

% Subsolar latitude:
      sbsllat = DELTA;

% Equation of time (degrees):
      ETDEG = L - ALPHA;
      NROT = round(ETDEG/360.);
      ETDEG = ETDEG - double(360*NROT);

% Apparent time (degrees):
      APTIME = UT/240. + ETDEG;
%          Earth rotates one degree every 240 s.

% Subsolar longitude:
      sbsllon = 180. - APTIME;
      NROT = round(sbsllon/360.);
      sbsllon = sbsllon - double(360*NROT);


