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]0001 function [z,sorad]=soradna1(yd,yr,long,lat) 0002 % SORADNA1: computes no-sky solar radiation and solar altitude. 0003 % [z,sorad]=SORADNA1(yd,yr,long,lat) computes instantaneous values of 0004 % solar radiation and solar altitude from yearday, year, and position 0005 % data. It is put together from expressions taken from Appendix E in the 0006 % 1978 edition of Almanac for Computers, Nautical Almanac Office, U.S. 0007 % Naval Observatory. They are reduced accuracy expressions valid for the 0008 % years 1800-2100. Solar declination computed from these expressions is 0009 % accurate to at least 1'. The solar constant (1368.0 W/m^2) represents a 0010 % mean of satellite measurements made over the last sunspot cycle (1979-1995) 0011 % taken from Coffey et al (1995), Earth System Monitor, 6, 6-10. Assumes 0012 % yd is either a column or row vector, the other input variables are scalars, 0013 % OR yd is a scalar, the other inputs matrices. 0014 % 0015 % INPUT: yd - decimal yearday (e.g., 0000Z Jan 1 is 0.0) 0016 % yr - year (e.g., 1995) 0017 % long - longitude (west is positive!) [deg] 0018 % lat - latitude [deg] 0019 % 0020 % OUTPUT: z - solar altitude [deg] 0021 % sorad- no atmosphere solar radiation [W/m^2] 0022 0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0024 % 3/8/97: version 1.0 0025 % 8/28/98: version 1.1 (vectorized by RP) 0026 % 8/5/99: version 2.0 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 0029 % get constants 0030 as_consts; 0031 0032 % convert yd to column vector if necessary 0033 [n,m]=size(yd); 0034 if m > n 0035 yd=yd'; 0036 end 0037 0038 % convert yearday to calender time 0039 gtime=greg2(yd,yr); 0040 0041 SC=gtime(:,6); 0042 MN=fix(gtime(:,5)); 0043 HR=fix(gtime(:,4)); 0044 D=fix(gtime(:,3)); 0045 M=fix(gtime(:,2)); 0046 Y=fix(gtime(:,1)); 0047 0048 % convert to new variables 0049 LONG=long; 0050 LAT=lat; 0051 0052 % two options - either long/lat are vectors, time is a scalar 0053 0054 if length(LONG)==1 & length(LAT)>1, 0055 LONG=LONG(ones(size(LAT))); 0056 elseif length(LONG)>1 & length(LAT)==1, 0057 LAT=LAT(ones(size(LAT))); 0058 end; 0059 0060 if length(SC)==1, 0061 osiz=ones(size(LONG)); 0062 SC=SC(osiz); 0063 MN=MN(osiz); 0064 HR=HR(osiz); 0065 D=D(osiz); 0066 M=M(osiz); 0067 Y=Y(osiz); 0068 elseif length(LONG)==1, 0069 LONG=LONG(ones(size(SC))); 0070 LAT=LAT(ones(size(SC))); 0071 end; 0072 0073 % constants 0074 DTR=3.14159265/180; 0075 RTD=1./DTR; 0076 0077 % compute Universal Time in hours 0078 UT = HR+(MN+SC./60.)./60; 0079 0080 % compute Julian ephemeris date in days (Day 1 is 1 Jan 4713 B.C.=-4712 Jan 1) 0081 JED=367.*Y-fix(7.*(Y+fix((M+9)./12))./4)+fix(275.*M./9)+D+1721013 + UT./24; 0082 0083 % compute interval in Julian centuries since 1900 0084 T=(JED-2415020.0)./36525; 0085 0086 % compute mean anomaly of the sun 0087 G=358.475833+35999.049750.*T-.000150.*T.^2; 0088 NG=fix(G./360); 0089 G=(G-NG.*360).*DTR; 0090 0091 % compute mean longitude of sun 0092 L=279.696678+36000.768920.*T+.000303.*T.^2; 0093 NL=fix(L./360); 0094 L=(L-NL.*360).*DTR; 0095 0096 % compute mean anomaly of Jupiter 0097 JUP=225.444651+2880.0.*T+154.906654.*T; 0098 NJUP=fix(JUP/360); 0099 JUP=(JUP-NJUP.*360).*DTR; 0100 0101 % compute longitude of the ascending node of the moon's orbit 0102 NM=259.183275-1800.*T-134.142008.*T+.002078.*T.^2; 0103 NNM=fix(NM./360); 0104 NM=(NM-NNM.*360+360).*DTR; 0105 0106 % compute mean anomaly of Venus 0107 V=212.603219+58320.*T+197.803875.*T+.001286.*T.^2; 0108 NV=fix(V/360); 0109 V=(V-NV.*360.).*DTR; 0110 0111 % compute sun theta 0112 THETA=.397930.*sin(L)+.009999.*sin(G-L)+.003334.*sin(G+L)... 0113 -.000208.*T.*sin(L)+.000042.*sin(2.*G+L)-.000040.*cos(L)... 0114 -.000039.*sin(NM-L)-.000030.*T.*sin(G-L)-.000014.*sin(2.*G-L)... 0115 -.000010.*cos(G-L-JUP)-.000010.*T.*sin(G+L); 0116 0117 % compute sun rho 0118 RHO=1.000421-.033503.*cos(G)-.000140.*cos(2*G)... 0119 +.000084.*T.*cos(G)-.000033.*sin(G-JUP)+.000027.*sin(2.*G-2.*V); 0120 0121 % compute declination 0122 DECL=asin(THETA./sqrt(RHO)); 0123 0124 % compute equation of time (in seconds of time) (L in degrees) 0125 L = 276.697+0.98564734.*(JED-2415020.0); 0126 L = (L - 360.*fix(L./360.)).*DTR; 0127 EQT = -97.8.*sin(L)-431.3.*cos(L)+596.6.*sin(2.*L)-1.9.*cos(2.*L)... 0128 +4.0.*sin(3.*L)+19.3.*cos(3.*L)-12.7.*sin(4.*L); 0129 EQT = EQT./60; 0130 L = L.*RTD; 0131 0132 % compute local hour angle 0133 GHA = 15.*(UT-12.) + 15.*EQT./60; 0134 LHA = GHA - LONG; 0135 0136 % compute radius vector 0137 RV=sqrt(RHO); 0138 0139 % compute solar altitude 0140 SZ=sin(DTR.*LAT).*sin(DECL)+cos(DTR.*LAT).*cos(DECL).*cos(DTR.*LHA); 0141 z=RTD.*asin(SZ); 0142 0143 % compute solar radiation outside atmosphere 0144 [n,m]=size(z); 0145 sorad=zeros(n,m); 0146 ii=z>0; 0147 sorad(ii)=(Solar_const./RV(ii).^2).*sin(DTR.*z(ii)); 0148