function [In,Id, Tr, To, Tg, Tw, Ta] = solflux(zdeg,w, p, k1, k2, l)

%SOLFLUX - COMPUTES GLOBAL SOLAR FLUX, CLEAR SKY

%  [In, Id, Tr, To, Tg, Tw, Ta] = solflux(zdeg, w, p, k1, k2, l)

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

% Calculate direct solar irradiance using a parameterization technique

% from M. Iqbal in "Physical Climatology for Solar and Wind Energy",

% 1988, World Scientific, pp. 196-242.

%

%Input:

% zdeg = zenith angle in degrees

% w = integrated water vapor (g/cm^2) typically 5 for TWP

% p == surface pressure in hPa (optional, default = 1013)

% k1 == aerosol optical thickness at 380 nm (optional, default = 0.0)

% k2 == aerosol optical thickness at 500 nm (optional, default = 0.0)

% l == ozone-layer thickness in cm(NTP) (optional, default = 0.3)

%Output:

% In == Direct-Normal flux at surface

% Id == Diffuse flux at surface

% Tr == Transmittance by Rayleigh scattering

% To == Transmittance by ozone

% Tg == Transmittance by uniformly mixed gases

% Tw == Transmittance by water vapor

% Ta == Transmittance by aerosol

%Internal constants:

% Isc == Solar constant = 1360 W/m^2

% mu = cos(zenith);

% Ma == Relative optical air mass

%

%reynolds 981029

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



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

% CONSTANTS

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

Uone = ones(size(w));		% ones vector

r2d = 180 ./ pi;

%p = 1013.;

%l = 0.3;

%k1 = 0.0;

%k2 = 0.0;

Isc = 1360.;

w0aer = 0.0;

fc = 0.84;

ag = 0.2;





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

% DEFAULT PARAMETERS

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

if nargin > 6, error('solflux parameter error'); end

if nargin <= 5, l = 0.3; end

if nargin <= 4, k2 = 0; end

if nargin <= 3, k1 = 0.0; end

if nargin <= 2, p = 1013; end

if nargin <= 1, error('solflux too few parameters'); end







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

% CLIP ALL ZENITH ANGLES TO <= 89.99 DEG.

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

ix = find( zdeg > 89.99);

if length(ix) > 0,

	zdeg(ix) = 89.99 .* ones(size(ix));

	mu(ix) = cos( zdeg(ix)  ./ r2d );

end

zrad = zdeg ./ r2d;

mu = cos(zrad);



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

% COMPUTE AIR MASS AND OTHER PROPERTIES

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

Mr = 1.0 ./ (mu + 0.15 .* (93.885 - zdeg ) .^ (-1.253) );

Ma = Mr .* (p ./ 1013.25);

U1 = w .* Mr;

U3 = l .* Mr;

ka = .2758 .* k1 + .35 .* k2;



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

% First calculate the transmittances due to various species

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

Tr = ones(size(zdeg));



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

% RAYLEIGH SCATTERING

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

ix = find(zdeg <= 85);

if length(ix) > 0,

	Tr(ix) = exp(-0.0903 .* (Ma(ix)) .^ 0.84 .* (1.0 + Ma(ix) - (Ma(ix)) .^1.01));

end

ix = find(zdeg > 85);

if length(ix) > 0,

	Tr(ix)  = exp(-(0.0738 + (zdeg(ix)-85) .* .03) .* (Ma(ix)) .^ 0.84)';

end



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

% OZONE SCATTERING

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

To = 1 - (.1611 .* U3 .* (1+139.48 .* U3) .^ (-.3035) - 0.002715 .* U3 ./(1+.044 .* U3 +.0003 .* U3 .* U3));



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

% UNIFORM GASES SCATTERING

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

Tg = exp(-.0127 .* (Ma) .^.26);



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

% WATERVAPOR SCATTERING

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

Tw = 1 - 2.459 .* U1 ./ ((1+79.034 .* U1) .^ .6828 + 6.385 .* U1);



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

% AEROSOL SCATTERING

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

Ta = exp(-ka .^.873 .* (1 + ka - ka .^ .7088) .* Ma .^ .9108);



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

% direct normal flux

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

In = Isc .* Tr .* To .* Tg .* Tw .*Ta;



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

% diffuse component

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

Taa = 1 - (1 - w0aer) .* (1 - Ma + Ma .^ 1.06) .* (1-Ta);



Idr = 0.79 .* Isc .* mu .* To .* Tg .* Tw .* Taa .* (.5 .* (1-Tr)) ./(1-Ma+Ma .^ 1.02);

	

Tas = Ta ./ Taa;



Ida = 0.79 .* Isc .* mu .* To .* Tg .* Taa .* (fc .* (1-Tas)) ./ (1-Ma+Ma.^ 1.02);



aa = 0.0685+(1-fc) .* (1-Tas);



Idm = (In .* mu + Idr + Ida) .* (ag .* aa ./ (1-ag .* aa));



Id = Idr + Ida + Idm;



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

% we get negative diffuse components??

%  This is a quick fix

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

ix = find(Id < 0);

if length(ix) > 0

	Id(ix) = 0;

end





return;



