% 1) Import Giovanni data using our new import function
%  Note that the imported values are in a variable called data.
import_Giovanni_file('Giovanni_chl_timeseries.dat');
chl = data;
import_Giovanni_file('Giovanni_SST_timeseries.dat');
SST = data;
% Create time vector of months (6 years total)
t = 1:12*6;
% Plot for a first visual inspection
figure
subplot(2,1,1)
plot(t,chl,'go-')
ylabel('chl [mg/m**3]')
subplot(2,1,2)
plot(t,SST,'ro-')
ylabel('SST [^oC]')
xlabel(' time [months since Jan 2003]')

% 2) Note that there is a missing data point in the 
% temperature series (it's a NaN in the vector).
% Remove NaN and fill the missing data point using
% interpolation 
bad = find(isnan(SST));
SST(bad) = [];
t_SST = t;
t_SST(bad) = [];
% Now interpolate to fill the empty spot
SST_new = interp1(t_SST,SST,t);
% Add this new SST data to our plot 
% (the second panel is still active)
hold on
plot(t,SST_new,'m.:')

% 3) Calculate cross-correlation
C = xcov(SST_new,chl,'coeff'); % option 'coeff' yields correlation
% Plot it for lags of -24 to +24 lags ( -2 year to +2 years )
N = length(t);
lags = [-24:24];
figure
plot(lags,C(N-lags),'-o')
hold on
% plot confidence limits
plot(lags,1.96./sqrt(N-lags-3),'g--')
plot(lags,-1.96./sqrt(N-lags-3),'g--')
xlabel('lag [months]')
ylabel('crosscorrelation')



