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;