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.   
 from the University of Maine Ocean Modelling Group http://rocky.umeoce.maine.edu/hjx/courses/SMS585/my2.5/air_sea/sunrise.m

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