function [rhr,rmin,shr,smin]=sunrise(mon,da,yr,lat,long)
% 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/28/98: version 1.1 (contributed by RP)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert calender time to julian yd
j=julianmd(yr,mon,da,0);
j0=julianmd(yr,1,1,0);
yd=j(:)-j0(:);

% compute solar altitude for entire day
dt=1./2880;

% we don't want abs(long)>180...
if long<-180, long=long+360; end;
if long>180,  long=long-360; end;

time=dt.*[0:2879]'+long/360; % have a whole day, beginning at midnight (near enough)
yday=yd(ones(1,2880),:)+time(:,ones(length(yd),1));

if length(yr)>1,
  yr=yr(:,ones(1,2880))';
end;

[z,sorad]=soradna1(yday(:),yr(:),long,lat);

z=reshape(z,2880,length(yd));
sorad=reshape(sorad,2880,length(yd));

[ir,jr]=find(sorad(1:2879,:)==0 & sorad(2:2880,:)>0);
[is,js]=find(sorad(2:2880,:)==0 & sorad(1:2879,:)>0);

srise=zeros(length(yd),1);
sset=any(sorad>0);

srise(jr)=yday(ir+(jr-1)*2880);
sset(js) =yday(is+(js-1)*2880);

rhr=fix(rem(srise,1)*24);
rmin=rem(rem(srise,1)*1440,60);
shr=fix(rem(sset,1)*24);
smin=rem(rem(sset,1)*1440,60);