function [daytime, nightime,rs]=daynight(lat,lon,date)

% Flags profiles by daytime/nightime depending on location and date of the
% year
% This function computes sunrise and sunset times using suncycle.m from the
% AIR_SEA TOOLBOX (version 2.0: 8/9/99) adapted by Rich Pawlowicz

% Input parameters:
% lat= scalar or vector with latitudes
% lon= scalar or vector with longitudes
% date =[ YYYY MM DD HH MM SS] 

% Output parameters:
% daytime =1 if profile time id during daylight hours or 0 otherwise
% nightime =1 if profile time is during night hours or 0 otherwise
% rs= is the output from suncycle (in GMT): sunrise time in the first column and
% sunset time in te second column. Special values: [0 24] is polar day, 
% [24 0] is polar night

% Jan 23, code modified to use +/- 1h from sunrise/sunset time to define
% day/night

np=length(date(:,1));
rs=nan(np,2);

% Compute sunset and sunrise time in GMT
for i=1:np;
    [rs(i,:)]=suncycle(lat(i),lon(i),date(i,1:3)); %% gives sunrise and sunset in GMT
    
end

% Convert time of day of profiles to decimal hours
dechr=date(:,4)+date(:,5)./60;

% Flag profiles by daytime/nightime
daytime=zeros(np,1); % using zeros to be able to convert to logical
nightime=zeros(np,1); % I tried with NaN and I get the same number of night vs day profiles

% % Original code using sunrise/sunset time to delimit day
% for i=1:length(dechr)
%     if rs(i,1)<rs(i,2) % sunrise time < sunset time
%         if dechr(i)>rs(i,1) && dechr(i)<rs(i,2) % profile time is > sunrise and < sunset
%             daytime(i)=1;
%         elseif dechr(i)<rs(i,1) || dechr(i)>rs(i,2)
%             nightime(i)=1;
%         end
%     elseif rs(i,1)>rs(i,2) %sunrise time > sunset time
%         if dechr(i)>rs(i,2) && dechr(i)<rs(i,1) % profile time is > sunset and <sunrise
%             nightime(i)=1;
%         elseif dechr(i)>rs(i,1) || dechr(i)<rs(i,2)
%             daytime(i)=1;
%         end
%     % Special values: [0 24] Polar Day, [24 0] is Polar Night]
%     elseif rs(i,1)==0 && rs(i,2)==24
%         daytime(i)=1; % Polar day
%     elseif rs(i,1)==24 && rs(i,2)==0
%         nightime(i)=1; % Polar night
%     end
% end

% New version of code using sunrise/sunset -/+ 1hr to delimit day:
for i=1:length(dechr)
    
    % Redefine sunrise and sunset times:
    if rs(i,1)>=1
        sunrise(i)=rs(i,1)-1;
    elseif rs(i,1)<1
        sunrise(i)=floor(rs(i,1)); % to avoid negative values
    end
        
    if rs(i,2)<=23
        sunset(i)=rs(i,2)+1;
    elseif rs(i,2)>23
        sunset(i)=ceil(rs(i,2));
    end
    
    
    % Daytime vs nightime
    if sunrise(i)<sunset(i) % sunrise time < sunset time
        if dechr(i)>sunrise(i) && dechr(i)<sunset(i) % profile time is > sunrise and < sunset
            daytime(i)=1;
        elseif dechr(i)<sunrise(i) || dechr(i)>sunset(i)
            nightime(i)=1;
        end
        
    elseif sunrise(i)>sunset(i) %sunrise time > sunset time (because it is GMT!)
        if dechr(i)>sunset(i) && dechr(i)<sunrise(i) % profile time is > sunset and <sunrise
            nightime(i)=1;
        elseif dechr(i)>sunrise(i) || dechr(i)<sunset(i)
            daytime(i)=1;
        end
        
    % Special values: [0 24] Polar Day, [24 0] is Polar Night]
    elseif rs(i,1)==0 && rs(i,2)==24
        daytime(i)=1; % Polar day
    elseif rs(i,1)==24 && rs(i,2)==0
        nightime(i)=1; % Polar night
    end
end

nightime=logical(nightime);
daytime=logical(daytime);



