function [out,data,ver] = skyline(az,el,strike,dip)

% skyline.m
%
% Syntax: [correction_factor,horizon,version] = skyline(azimuths,elevations,strike,dip);
%
% This function calculates the shielding factor for cosmogenic isotope 
% production on a surface which is either dipping, shielded by high 
% topography, or both. 
%
% The first two arguments are vectors of similar 
% length containing the azimuth (0-360) and elevation (0-90) or points 
% on the horizon. This presumes that in the field the horizon was 
% approximated by a series of points with straight lines between them. 
% Note that this is incompatible with field procedures in which 
% the horizon is approximated by the average rise angle in a series 
% of equal sectors. This latter procedure is inappropriate anyway 
% because the relationship between rise angle and shielding is nonlinear.
%
% The second two arguments are optional and are scalars describing the 
% strike and dip of the surface sampled. Use the convention that 
% the dip is down to the right of the strike direction, i.e. strike = 0,
% dip = 45 means that the surface is dipping 45 down to the E. 
%
% If you just want the shielding for a dipping surface enter
% empty vectors ([]) for the first two arguments.
%
% Note that this function first compares the "horizon" created by a
% dipping surface (second 2 arguments) with the far-field horizon 
% (first 2 arguments) and calculates a composite horizon before 
% determining the shielding factor. Thus the shielding due to a 
% dipping surface is never "double-counted" with a distant horizon.
%
% Uses highly simplistic numerical integration -- approximates
% horizon by a series of 1-degree sectors and uses an integration
% formula to calculate the shielding for each of these sectors,
% then sums them. 
% 
% The output argument correction_factor is the ratio of nuclide production 
% at the shielded site to that at an unshielded site at the same location.
%
% The output argument 'horizon' is a vector containing the horizon angle in degrees
% at 1-degree increments. Useful for plotting purposes. 
% 
% The output argument 'ver' is a string giving the version number of this function. 
%
% Written by Greg Balco -- UW Cosmogenic Nuclide Lab
% balcs@u.washington.edu
% First version, Feb. 2001
% Revised March, 2006
% Part of the CRONUS-Earth online calculators: 
%      http://hess.ess.washington.edu/math
% Further revised April 2018 to clean up some of the error handling. 
%
% Copyright 2001-2007, University of Washington
% Subsequently Greg Balco
% All rights reserved
% Developed in part with funding from the National Science Foundation.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License, version 2,
% as published by the Free Software Foundation (www.fsf.org).
%
% Note error fix in angular increment definition 20220629. 
 

% what's the version

ver = '2.0';

% Sort missing arguments

if nargin < 4; dip = 0; end;
if nargin < 3; strike = 0; end;

if isempty(dip) || isempty(strike); strike = 0; dip = 0; end;
% So now if nothing was submitted for strike and/or dip, there should be
% zero placeholders present. 

% set up angular framework -- 1-degree increments...

angles = [(pi./180):pi./180:(2.*pi)];

% Finally, check to make sure something was actually submitted
if dip == 0 && isempty(az) && isempty(el);
    % Case nothing to calculate. Return zeros. 
    out = 1; data = zeros(size(angles)); 
    return;
end;

% Proceed with error checks

if max(max(az)) > 360;
	error('skyline.m: azimuths 0 to 360 please');
elseif max(max(el)) > 90;
	error('skyline.m: elevations 0 to 90 please');
elseif length(az) ~= length(el);
	error('skyline.m: vectors the same length please');
elseif strike < 0 || strike > 360;
    error('skyline.m: strike must be 0-360');
elseif dip < 0 || dip > 90;
    error('skyline.m: dip must be 0-90');
end;

% Step 1. Determine if dipping surface needed; create horizon
% from such surface.

if dip > 0;

	% convert to radians
	
	strikeR = (strike/360) * (2*pi);
	dipR = (dip/360) * (2*pi);
	
	% calculate "horizon"
	% updip direction = strike direction - pi/2 by convention
	% a relative horiz angle to dip direction
	
	a = angles - (strikeR - (pi/2)); 
	
	horiz1 = atan(tan(dipR) .* cos(a));
	
	% remove entries less than 0, i.e. below horizon
	
	horiz1(find(horiz1 < 0)) = zeros(size(find(horiz1 < 0)));
	
else
    % Case zero dip
	horiz1 = zeros(size(angles));
	
end;

% Step 2. Interpolate between supplied horizon points. Note linear
% interpolation in degree space might not be right, but you'd have to make
% lots of other assumptions about the orientation of the topography to do
% it any better way. 

if ~(isempty(az) && isempty(el));
	
    % Azimuth, elv data present

	% convert to radians
	azR = (az/360) * (2*pi);
	elR = (el/360) * (2*pi);

	% sort in ascending order
	[azR,si] = sort(azR);
	elR = elR(si);

	% pad for interpolation;
	azR2(2:(length(azR)+1)) = azR;
	elR2(2:length(elR)+1) = elR;
	azR2(1) = azR(length(azR)) - (2*pi); 
	elR2(1) = elR(length(elR));
	azR2(length(azR)+2) = azR(1) + (2*pi); 
	elR2(length(azR)+2) = elR(1);

	% interpolate;

	horiz2 = interp1(azR2,elR2,angles);

	% flip?

	horiz2 = horiz2;
	
else
    % Case no azimuth, elevation submitted
	horiz2 = zeros(size(angles));
	
end;

% Step 3. Make composite horizon. 

horiz = max([horiz1;horiz2]);

% Step 4. Integrate shielding function.
% Here we do the simplest possible thing - integrate in 
% 1-degree increments and add them up.

% sector angle = 1 degree;

B = (1/360)*(2*pi);

% integration formula

S = (B/(2*pi)) .* (sin(horiz).^3.3);

% above formula gives percent shielding; subtract from 1 to give corr. factor --

out = 1 - sum(S);

data = horiz.*180/pi;


	
	
	
