sunrise

 SUNRISE: computes sunrise and sunset times for specified day and location. 
 [rhr,rmin,shr,smin] = SUNRISE(mon,da,yr,lat,long) computes the time 
 of sunrise rhr:rmin and sunset shr:smin to the nearest minute in GMT 
 for a calendar day(s) and a specified (scalar) position.   
 
 INPUT:  mon - month (e.g., Jan is 1)
         da - day (e.g., Jan 10th is 10)
         yr - year (e.g., 1995)
         lat - latitude [deg]
         long - longitude (west is positive)  [deg] 

 OUTPUT: rhr,rmin  - sunrise in GMT hours and minutes
          shr,smin  - sunset in GMT hours and minutes.

0001 function [rhr,rmin,shr,smin]=sunrise(mon,da,yr,lat,long)
0002 % SUNRISE: computes sunrise and sunset times for specified day and location.
0003 % [rhr,rmin,shr,smin] = SUNRISE(mon,da,yr,lat,long) computes the time
0004 % of sunrise rhr:rmin and sunset shr:smin to the nearest minute in GMT
0005 % for a calendar day(s) and a specified (scalar) position.
0006 %
0007 % INPUT:  mon - month (e.g., Jan is 1)
0008 %         da - day (e.g., Jan 10th is 10)
0009 %         yr - year (e.g., 1995)
0010 %         lat - latitude [deg]
0011 %         long - longitude (west is positive)  [deg]
0012 %
0013 % OUTPUT: rhr,rmin  - sunrise in GMT hours and minutes
0014 %          shr,smin  - sunset in GMT hours and minutes.
0015 
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 % 8/28/98: version 1.1 (contributed by RP)
0018 % 8/5/99: version 2.0
0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 
0021 % convert calender time to julian yd
0022 j=julianmd(yr,mon,da,0);
0023 j0=julianmd(yr,1,1,0);
0024 yd=j(:)-j0(:);
0025 
0026 % compute solar altitude for entire day
0027 dt=1./2880;
0028 
0029 % we don't want abs(long)>180...
0030 if long<-180, long=long+360; end;
0031 if long>180,  long=long-360; end;
0032 
0033 time=dt.*[0:2879]'+long/360; % have a whole day, beginning at midnight (near enough)
0034 yday=yd(ones(1,2880),:)+time(:,ones(length(yd),1));
0035 
0036 if length(yr)>1,
0037   yr=yr(:,ones(1,2880))';
0038 end;
0039 
0040 [z,sorad]=soradna1(yday(:),yr(:),long,lat);
0041 
0042 z=reshape(z,2880,length(yd));
0043 sorad=reshape(sorad,2880,length(yd));
0044 
0045 [ir,jr]=find(sorad(1:2879,:)==0 & sorad(2:2880,:)>0);
0046 [is,js]=find(sorad(2:2880,:)==0 & sorad(1:2879,:)>0);
0047 
0048 srise=zeros(length(yd),1);
0049 %SO: bugfix
0050 sset=double(any(sorad>0));
0051 sset=any(sorad>0);
0052 
0053 srise(jr)=yday(ir+(jr-1)*2880);
0054 sset(js) =yday(is+(js-1)*2880);
0055 
0056 rhr=fix(rem(srise,1)*24);
0057 rmin=rem(rem(srise,1)*1440,60);
0058 shr=fix(rem(sset,1)*24);
0059 smin=rem(rem(sset,1)*1440,60);
0060 
0061 
0062 
0063