function [sw,sz,saz,dir,sky] = SolarRadiance(alat,alon,jd, tutc,watvap, p,k1, k2, l)

%SOLARRADIANCE --	Radince, zenith, and azimuth

%=====================================================================

%[sw,sz,saz,dir,sky] = SolarRadiance(alat,alon,jd, tutc,watvap, p, k1, k2, l)

%

%%% COMPUTE THE SOLAR RADIANCE

%	sw = total solar flux through a horizontal surface

% 	sz = zenith declination (deg)

%	saz = solar azimuth (deg)

%  dir = direct beam irradiance on a surface normal to the beam.

%  sky = diffuse component on a horizontal surface

%INPUT -

% alat, alon = latitude and longitude in degrees

%  jday is the year day (utc)

%  tutc = f.p. hour of the day

%  watvap = integrated water vapor (g/m^2) -- use 5 for TWP, default = 3.

%  p = pressure, default = 1013 hPa

%  k1, k2 = aerosol absorbtion coefs, default = .3

% l = ozone thickness, default = 0.3;

%

%Uses function "solflux" to compute the direct beam and diffuse radiation.

%assumptions:

% transmittance from aerosols Ta = 0;

% sfc pressure = 1023 hPa

% ozone = 0.3 (see solflux.m comments);

%

%The ephem function came from P. Minnett

%Reynolds, 971029

%=====================================================================



% TESTING

%alat = 0;  alon = 162;  jd = 180;%julian(1997,2,5);

%tutc = [0:.5:24]';

%watvap = 5;

%fprintf('Testing solarradiance\n');



%fprintf('SolarRadiance: nargin = %d\n',nargin);

if nargin < 9

	l = 0.3;

end

if nargin < 8

	k2 = 0.3;

end

if nargin < 7

	k1 = 0.3;

end

if nargin < 6

	p = 1013;

end

if nargin < 5

	watvap = 3;

end

if nargin < 4,

	error('too few params');

end



%% constants

d2r = pi / 180;  r2d = 1. ./ d2r;

s0 = 1367;

s1 = -105;   	s2 = 596.2;   	s3 = 4.3;   	s4 = -12.7;

c1 = -429.3;   	c2 = -2.1;   	c3 = 19.3;



%% coordinates in radians

rlat = alat .* d2r;

rlon = alon .* d2r;

colatr = (90 - alat) * d2r;



%=================

%time check, if we are making a computation at a single time,

% we must do it for that time plus a small increment so,

%in the end we can adjust the azimuth

%==================

an = jd - 0.5 + tutc ./ 24;

if length(an) == 1

	tutc = [tutc; tutc + 1e-6];

	jd = [jd; jd];

	an = jd - 0.5 + tutc ./ 24;

	tlength = 1;

else

	tlength = length(an);

end





al = 280.460 + 0.9856474 .* an;

al = rem(al,360);



ix = find(al < 0 );

if length(ix) > 0

	al(ix) = al(ix) + 360;

end



%% CALCULATE MEAN ANOMOLY

g = 357.528 + 0.9856003 .* an;

g = rem(g,360);

ix = find(g < 0);

if length(ix) > 0,

	g(ix) = g(ix) + 360 .* ones(size(ix));

end



gr = g * d2r;



% SUN EARTH DISTANCE AND NORMAL INSOLATION

r = 1.00014 - 0.01671 .* cos(gr) - 0.00014 .* cos(2 * gr);

s = s0 ./ (r .* r);

%%fprintf('E-S distance, s0 = %10.3f,   %10.3f\n', r,s);



%% CALCULATE ECLIPTIC LONGITUDE

alambda = al + 1.915 .* sin(gr) + 0.020 .* sin(2*gr);

rlambda = alambda * d2r;



%% ECLIPTIC OBLIQUITY

epsilon = 23.439 - 0.0000004 .* an;



%% RIGHT ASCENSION

te = tan(epsilon * 0.5 .* d2r);

te = te .* te;

alpha = alambda - r2d .* te .* sin( 2 * rlambda ) +	0.5 .* r2d .* te .* te .* sin(4 * rlambda);



%% SOLAR DECLINATION FOR THIS DAY

sin_dec = sin(epsilon * d2r) .* sin(rlambda);

solar_dec = asin(sin_dec);

cos_dec = cos(solar_dec);



eq_time = ( al - alpha ) .* 4 ./ 60;

s_hour = tutc + eq_time;

hour_ang = pi * ( s_hour - 12 ) / 12 + rlon;

cos_h = cos(hour_ang);



c_solar_z = sin(rlat) .* sin_dec + cos(rlat) .* cos_dec .* cos_h;

s_z = acos(c_solar_z);





% CORRCT AZIMUTH -- locate slope of the zenith angle

zz = 0;

if length(s_z) > 1,

	zz = diff(s_z);                  % slope of s_z, increasing or decreasing

	zz = [zz; zz(length(zz))];       % to make the lengths the same

end



%%% ZENITH ANGLE

sz = s_z * r2d;

											% Correct for E-W



%%% AZIMUTH ANGLE

s_az = ( cos(pi * 0.5 - solar_dec) - cos(colatr) .* c_solar_z)./ ( sin(colatr) .* sin(s_z) );



saz = acos(s_az) * r2d;





% AZIMUTH CORRECTION -- where zenith angle is increasing

if length(sz) > 0,

	ix = find(zz > 0);

	if length(ix) > 0,

		saz(ix) = 360 .* ones(size(ix)) - saz(ix);

	end

end



%SINGLE POINT CALC

if tlength == 1

	saz = saz(1);

	sz = sz(1);

end





%%% SHORTWAVE IRRADIANCE

% the zenith angle = sz (degs),

% the integrated water vapor is set in the input parameters (TWP = 5).

%

%

if watvap(1) < 0,

   dir = s;

   sky = 0;

	sw = s .* cos(sz .* d2r);

else

	[dir,sky] = solflux(sz,watvap,p,k1,k2,l);

	sw = sky + dir .* cos(sz .* d2r);

end



ix = find(sw < 0);

if length(ix) > 0,

	sw(ix) = zeros(size(ix));

end

ix = find(sky < 0);

if length(ix) > 0,

	sky(ix) = zeros(size(ix));

end



return



