function [F_solar,Ls] = calc_insolation(ecc,epsilon,Lsp,S0,lat,Ls);
%% HEADER: calc_insolation.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%Description: Calculates daily averaged solar insolation values for given
%              orbital elements, latitude, and Ls values
%
%%Inputs:
%        lat = vector of latitudes in degrees (optiona) [default is -90:90]]
%        ecc = orbital eccentricity
%        epsilon = obliquity  [radians]
%        Lsp = Ls of perihelion [radians]
%        S0 = Solar Constant [W/m^2]
%        Ls = vector of Ls values in degrees (optional) [defualt is to use
%        NOTE: If Orbital Elements are not input, J2000 values for
%              Saturn/Titan are used
%
%%Outputs:
%         F_solar = daily averaged solar flux in W/m^2
%         Ls = Ls values in degrees (same as input or default)
%
%%Calls: None
%
%%Required Files: None
%
%%Source:
%        Berger, 1978
%
%%Author: Alexander G. Hayes (hayes@gps.caltech.edu)
%
%%Change Log: 
%            09/13/08: Draft Code Escalated to V1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% BEGIN CODE: calc_insolation.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%lat defaults to a vector from the south to north pole in one degree step
if ~exist('lat','var');
    lat = [-90:90];
elseif length(lat)==0;
    lat = [-90:90];
end;
lat = lat.*pi./180;

%Orbital Paramters Default to Titan/Saturn at J2000
if (~exist('ecc','var') | ~exist('epsilon','var') | ...
        ~exist('Lsp','var') | ~exist('S0','var') );
    disp(['calc_insolation: Orbital Elements Not Input! Using Titan/Saturn at J2000!']);
    %J2000 = 2451545.5, Elements Calculated from JPL Horizons
    ecc = 0.0558249946541988124;
    epsilon = 0.476771591767290983 - 0.33.*pi./180;
    Lsp = 4.86583101394663888818;
    S0 = 14.9;
elseif (length(ecc) == 0) | (length(epsilon) == 0) | ...
        (length(Lsp) == 0) | (length(S0) == 0);
        disp(['calc_insolation: Orbital Elements Not Input! Using Titan/Saturn at J2000!']);
    %J2000 = 2451545.5, Elements Calculated from JPL Horizons
    ecc = 0.0558249946541988124;
    epsilon = 0.476771591767290983 -0.33.*pi./180;
    Lsp = 4.86583101394663888818;
    S0 = 14.9;
end;
    
%Check to make sure orbital elements are singleton values
if (length(ecc) > 1) | (length(epsilon) > 1) | ...
        (length(Lsp) > 1) | (length(S0) > 1);
    disp('calc_insolation: Only 1-Element Vectors of Orbital Parameters are Supported!');
    disp('calc_insolation: Using Index Position 1 on all input elements...');
    ecc = ecc(1);
    Lsp = Lsp(1);
    epsilon = epsilon(1);
    S0 = S0(1);
end;

%Ls defaults to a vector of 0 to 360
if ~exist('Ls','var');
    Ls = linspace(0,2.*pi,1000);
elseif length(Ls)==0;
    Ls = linspace(0,2.*pi,1000);
else;
    if max(Ls(~isnan(Ls))) > 2*pi | min(Ls(~isnan(Ls))) < -pi;
        Ls = Ls.*pi./180;
    else;
        Ls = Ls; %Units are now radians
    end;    
end;

%make a grid if we have both lat and Ls varying
if ( (length(lat) > 1) & (length(Ls) > 1) );
    lat2 = zeros(length(Ls),length(lat));
    Ls2 =  zeros(length(Ls),length(lat));
    for i=1:length(Ls);
        lat2(i,:) = lat;
    end;
    for i=1:length(lat);      
        Ls2(:,i) = Ls;
    end;
    lat = lat2;
    Ls = Ls2;   
    clear lat2 Ls2
end;

%calculate the declination angle of the sun
% From Berger, 1978
dec=asin(sin(epsilon).*sin(Ls)); 

%calculate the absolute value of the hour angle at sunrise / sunset
% From Berger, 1978
Ho=acos(-tan(lat).*tan(dec)); 

%correct for latitudes where there is no sunrise or sunset
% From Berger, 1978, Equations 8 & 9
ind = find( (abs(lat) >= (pi/2 - abs(dec))) & (lat.*dec > 0) );
Ho(ind) = pi;
ind = find( (abs(lat) >= (pi/2 - abs(dec))) & (lat.*dec <= 0) );
Ho(ind) = 0;

%calculate the normalized distance to the sun
% From Berger, 1978
p = (1-ecc.^2)./(1+ecc.*cos(Ls-Lsp));

%Calculate the Solar Insolation (W/m^2)
% From Berger 1978 eq (10)
F_solar = S0./pi./p.^2 .*( Ho.*sin(lat).*sin(dec) + cos(lat).*cos(dec).*sin(Ho) );

%% END CODE: calc_insolation.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%