function [cosz] = coszeta(time, lat, declmax, taudecl, taurot)
%cos(zeta) -- function calculates cosine of zenith angle as a function of time 
% for a given latitude on a planet in a circular orbit
% according to formula cos(zeta)=cos(lat)*cos(hour)*cos(decl)+sin(decl)*sin(lat)
% (and filtered for only positive / zero values)
% time = array of time values -- units unspecified, but taudecl, taurot
% must be same units

% lat = latitude, in degrees
% declmax = maximum declination, in degrees
% taudecl = period of revolution (e.g. length of year)
% taurot = period of planeraty rotation (e.g. length of day)

cosz=cosd(lat).*cosd(declmax*(cos(2*pi/(taudecl)*(time)))).*cos(2*pi/taurot*time)+...
     sind(lat).*sind(declmax*(cos(2*pi/(taudecl)*(time))));
 
cosz=cosz.*(cosz>=0); 

end

