function [I]=inso2(kyear,lat,day,orbit)

% Usage:   [I]=inso2(kyear,lat,day)
%
% Purpose: Computes daily average insolation as function of day and
%   latitude at any point during past 5 million years.
%
% Input arguments:
%  kyear: Thousands of years before present (0 to 5000).
%  lat: Latitude in degrees (-90 to 90). 
%  day: calendar day (1-365.2422) is proportional to the amount of time that
%    has elapsed since the vernal equinox (always occuring at day=80). The
%    difference between calendar dates and solar longitude dates is related to
%    the fact that the Earth's elliptical orbit sweeps through angles more
%    slowly when the Earth is farther from the Sun (Kepler's second law).
%
% Output: Daily average solar radiation in W/m^2. 
%
% Required file: inso2.mat
%
% References: 
%  Berger A. and Loutre M.F. (1991). Insolation values for the climate of
%   the last 10 million years. Quaternary Science Reviews, 10(4), 297-317.
%  Berger A. (1978). Long-term variations of daily insolation and
%   Quaternary climatic changes. Journal of Atmospheric Science, 35(12),
%   2362-2367.
%
% Coded by Ian Eisenman, Harvard 2006


%Check input arguments
if nargin<3, disp('I = inso(kyear,lat,day)'); return; end

%Get orbital parameters
if nargin<4, 
  [ecc,epsilon,omega]=orbital_parameters(kyear); % function is below in this file
  obliquity=epsilon*180/pi;
  long_perh=omega*180/pi;
else,
  ecc=orbit(1);
  epsilon=orbit(2); 
  omega=orbit(3); 
end;

  %Calculate insolation 
lat=lat*pi/180; % latitude; phi in Berger 1978
  % calculate true longitude (lambda) from calendar day by computing
  % mean longitude, based on approximation from Berger 1978 section 3
delta_lambda_m=(day-80)*2*pi/365.2422;
beta=(1-ecc.^2).^(1/2);
lambda_m0=-2*( (1/2*ecc+1/8*ecc.^3).*(1+beta).*sin(-omega)-1/4*ecc.^2.*(1/2+beta).*sin(-2*omega)+1/8*ecc.^3.*(1/3+beta).*(sin(-3*omega)) );
lambda_m=lambda_m0+delta_lambda_m;
lambda=lambda_m+(2*ecc-1/4*ecc.^3).*sin(lambda_m-omega)+(5/4)*ecc.^2.*sin(2*(lambda_m-omega))+(13/12)*ecc.^3.*sin(3*(lambda_m-omega));
So=1365; % solar constant (W/m^2)
delta=asin(sin(epsilon).*sin(lambda)); % declination of the sun

for ctlat=1:length(lat),
  Ho=acos(-tan(lat(ctlat)).*tan(delta)); % hour angle at sunrise/sunset
  % no sunrise or no sunset: Berger 1978 eqn (8),(9)
  Ho( ( abs(lat(ctlat)) >= pi/2 - abs(delta) ) & ( lat(ctlat).*delta > 0 ) )=pi;
  Ho( ( abs(lat(ctlat)) >= pi/2 - abs(delta) ) & ( lat(ctlat).*delta <= 0 ) )=0;
  % Insolation: Berger 1978 eq (10)
  I(ctlat,:)=So/pi*(1+ecc.*cos(lambda-omega)).^2 ./ (1-ecc.^2).^2 .* ( Ho.*sin(lat(ctlat)).*sin(delta) + cos(lat(ctlat)).*cos(delta).*sin(Ho) );
end;

function [ecc,epsilon,omega] = orbital_parameters(kyear)
load inso2.mat;
kyear0=-m(:,1); % kyears before present for data (kyear0>=0);
ecc0=m(:,2); % eccentricity
  % add 180 degrees to omega (see lambda definition, Berger 1978 Appendix)
omega0=m(:,3)+180; % longitude of perihelion (precession angle)
omega0=unwrap(omega0*pi/180)*180/pi; % remove discontinuities (360 degree jumps)
epsilon0=m(:,4); % obliquity angle
ecc=interp1(kyear0,ecc0,kyear,'spline',NaN); % eccs means array of ecc values
omega=interp1(kyear0,omega0,kyear,'spline',NaN) * pi/180;
epsilon=interp1(kyear0,epsilon0,kyear,'spline',NaN) * pi/180;