soradna1.m

function [z,sorad]=soradna1(yd,yr,long,lat)
% SORADNA1: computes no-sky solar radiation and solar altitude.
% [z,sorad]=SORADNA1(yd,yr,long,lat) computes instantaneous values of 
% solar radiation and solar altitude from yearday, year, and position 
% data. It is put together from expressions taken from Appendix E in the
% 1978 edition of Almanac for Computers, Nautical Almanac Office, U.S.
% Naval Observatory. They are reduced accuracy expressions valid for the
% years 1800-2100. Solar declination computed from these expressions is
% accurate to at least 1'. The solar constant (1368.0 W/m^2) represents a 
% mean of satellite measurements made over the last sunspot cycle (1979-1995) 
% taken from Coffey et al (1995), Earth System Monitor, 6, 6-10.  Assumes 
% yd is either a column or row vector, the other input variables are scalars,
% OR yd is a scalar, the other inputs matrices.
%
%  INPUT:   yd   - decimal yearday (e.g., 0000Z Jan 1 is 0.0)
%           yr   - year (e.g., 1995)
%           long - longitude (west is positive!) [deg] 
%           lat  - latitude  [deg]
%
%  OUTPUT:  z    - solar altitude [deg]
%           sorad- no atmosphere solar radiation  [W/m^2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/28/98: version 1.1 (vectorized by RP)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get constants
as_consts; 

% convert yd to column vector if necessary
[n,m]=size(yd);
if m > n
  yd=yd';
end

% convert yearday to calender time
gtime=greg2(yd,yr);

SC=gtime(:,6);
MN=fix(gtime(:,5));
HR=fix(gtime(:,4));
D=fix(gtime(:,3));
M=fix(gtime(:,2));
Y=fix(gtime(:,1));

% convert to new variables
LONG=long;
LAT=lat;

% two options - either long/lat are vectors, time is a scalar

if length(LONG)==1 & length(LAT)>1,
 LONG=LONG(ones(size(LAT)));
elseif length(LONG)>1 & length(LAT)==1,
 LAT=LAT(ones(size(LAT)));
end;
  
if length(SC)==1,
 osiz=ones(size(LONG));
 SC=SC(osiz);
 MN=MN(osiz);
 HR=HR(osiz);
 D=D(osiz);
 M=M(osiz);
 Y=Y(osiz);
elseif length(LONG)==1,
 LONG=LONG(ones(size(SC)));
 LAT=LAT(ones(size(SC)));
end;

% constants
   DTR=3.14159265/180;
   RTD=1./DTR;

% compute Universal Time in hours
   UT = HR+(MN+SC./60.)./60;

% compute Julian ephemeris date in days (Day 1 is 1 Jan 4713 B.C.=-4712 Jan 1)
  JED=367.*Y-fix(7.*(Y+fix((M+9)./12))./4)+fix(275.*M./9)+D+1721013 + UT./24;

% compute interval in Julian centuries since 1900
    T=(JED-2415020.0)./36525;

% compute mean anomaly of the sun
    G=358.475833+35999.049750.*T-.000150.*T.^2;
   NG=fix(G./360);
    G=(G-NG.*360).*DTR;

% compute mean longitude of sun
    L=279.696678+36000.768920.*T+.000303.*T.^2;
   NL=fix(L./360);
    L=(L-NL.*360).*DTR;

% compute mean anomaly of Jupiter
  JUP=225.444651+2880.0.*T+154.906654.*T;
 NJUP=fix(JUP/360);
  JUP=(JUP-NJUP.*360).*DTR;

% compute longitude of the ascending node of the moon's orbit
   NM=259.183275-1800.*T-134.142008.*T+.002078.*T.^2;
  NNM=fix(NM./360);
   NM=(NM-NNM.*360+360).*DTR;

% compute mean anomaly of Venus
    V=212.603219+58320.*T+197.803875.*T+.001286.*T.^2;
   NV=fix(V/360);
    V=(V-NV.*360.).*DTR;

% compute sun theta
 THETA=.397930.*sin(L)+.009999.*sin(G-L)+.003334.*sin(G+L)...
     -.000208.*T.*sin(L)+.000042.*sin(2.*G+L)-.000040.*cos(L)...
     -.000039.*sin(NM-L)-.000030.*T.*sin(G-L)-.000014.*sin(2.*G-L)...
     -.000010.*cos(G-L-JUP)-.000010.*T.*sin(G+L);

% compute sun rho
  RHO=1.000421-.033503.*cos(G)-.000140.*cos(2*G)...
     +.000084.*T.*cos(G)-.000033.*sin(G-JUP)+.000027.*sin(2.*G-2.*V);

% compute declination
   DECL=asin(THETA./sqrt(RHO));

% compute equation of time (in seconds of time) (L in degrees)
    L = 276.697+0.98564734.*(JED-2415020.0);
    L = (L - 360.*fix(L./360.)).*DTR;
  EQT = -97.8.*sin(L)-431.3.*cos(L)+596.6.*sin(2.*L)-1.9.*cos(2.*L)...
         +4.0.*sin(3.*L)+19.3.*cos(3.*L)-12.7.*sin(4.*L);
  EQT = EQT./60;
    L = L.*RTD;

% compute local hour angle
  GHA = 15.*(UT-12.) + 15.*EQT./60;
  LHA = GHA - LONG;

% compute radius vector
   RV=sqrt(RHO);

% compute solar altitude
   SZ=sin(DTR.*LAT).*sin(DECL)+cos(DTR.*LAT).*cos(DECL).*cos(DTR.*LHA);
    z=RTD.*asin(SZ);

% compute solar radiation outside atmosphere
[n,m]=size(z);
sorad=zeros(n,m);
ii=z>0;
sorad(ii)=(Solar_const./RV(ii).^2).*sin(DTR.*z(ii));