function [out,rlags,RU,RL] = XcorelNaN(yy1_in,yy_in, nlags_in, fig)
%Calculates the lagged cross correlation of two variables, yy1, and yy,
%over nlags:
% out(tau) = integral{(yy1(t) * yy(t+tau)} dt

% note + lag will be yy1 leading, negative will be yy leading
% fig = figure number to plot, fig=0 means no plots

% make sure the vectors are equal size
sz1=size(yy1_in,1);
sz2=size(yy_in,1);
npts=min(sz1,sz2);
yy1=yy1_in(1:npts);
yy = yy_in(1:npts);

%remove the mean values from yy1 and yy, assuming there there be NaN's 
yy1=yy1-nanmean(yy1);
yy =yy -nanmean(yy);

% make sure the number of lags is <= number of points-4 (to ensure at least 2 points
% overlapped in the correlation integral.
if nlags_in==0
    nlags_in=npts-4;
end
nlags=min(nlags_in, npts-4);

% just use this command to set the variable LAGS.
[C,LAGS] = xcorr(yy,yy1,nlags,'coeff');
rlags = LAGS';

% set up the output arrays
out=zeros(2*nlags+1,1);p=zeros(2*nlags+1,1);RL=zeros(2*nlags+1,1);RU=zeros(2*nlags+1,1);

% first the negative lags, where we use corrcoef 'rows' with either 'complete' or
% 'pairwise' (for just 2 columns it does not matter) to eliminate any NaN's
% in the correlation. This will yield Pearson's linear correlation coefficient. With
% corrcoef we can also calculate the probability and upper and lower
% confidence bounds for alpha = 0.05 -> 95% confidence
for i=nlags+1:-1:2
    [cccoef ppp cclow ccup ]= corrcoef(yy1(i:npts),yy(1:npts-i+1),'alpha',0.05,'rows','pairwise');
    out(nlags+2-i)=cccoef(2,1); p(nlags+2-i)=ppp(2,1); RU(nlags+2-i)=ccup(2,1); RL(nlags+2-i)=cclow(2,1); 
end
%then the positive
for i=1:nlags+1
    [cccoef ppp cclow ccup ]= corrcoef(yy(i:npts),yy1(1:npts-i+1),'alpha',0.05,'rows','pairwise');
    out(i+nlags)=cccoef(2,1); p(i+nlags)=ppp(2,1); RU(i+nlags)=ccup(2,1); RL(i+nlags)=cclow(2,1);
end

% plot the output if requested
if fig ~= 0
    figure(fig)
    % plot the confidence band as a filled colour between the upper and lower
    % confidence interval of the estimated coefficients, with the actual values
    % on top
    fill([LAGS,fliplr(LAGS)],[RU',fliplr(RL')],[.702 0.78 1.0],'EdgeColor',[0.702 0.78 1.0]); hold;
    plot(-nlags:nlags,out,'k')
    plot([-(nlags+5) nlags+5], [0  0]);
    plot([0  0], [ -1 1 ]);
    xlim([-(nlags+5) (nlags+5)]);
    xlabel('Lag (days)');
    ylabel('Correlation Coefficient, r, and 95% confidence interval');
    hold
end

end
