function latency= Latency(GPS,Echo)
%%% determine latency between echosounder and GPS

% Load echosounder and GPS values as acquired by hypack into workspace
% first by running readHypackraw.m
% e.g. [GPS,Echo]=readHypackraw('00X_XXXX.RAW')

%% Show figure echosounder values
figure; 
 plot(Echo.GPSMatlabtime,-Echo.depth,'g')
 datetick
 xlabel('time')
 ylabel('depth below tranducer (m)')

%% Select a sample by picking an area of interest (AOI)
disp('Pick an area within the survey to determine the latency')
AOI=ginput(2); % mouse click area of interest
AOIbegin=min(AOI(:,1));
AOIend=max(AOI(:,1));
 
%%% Resample GPS on Echosounder (10Hz) time
 GPS_spline = interp1(GPS.Matlabtime, GPS.elevation, Echo.GPSMatlabtime, 'spline');  

%%% find sample indices
 sample1= find(Echo.GPSMatlabtime > AOIbegin & Echo.GPSMatlabtime < AOIend);
 
%%% Show full Sample 
 figure; 
 plot(Echo.GPSMatlabtime(sample1),GPS_spline(sample1)-43.11,'r')
 hold on;
 plot(Echo.GPSMatlabtime(sample1),-Echo.depth(sample1),'g')
datetick
grid on;
xlabel('time')
ylabel('Echosounder (g) and interpolated GPS (r) values')
title(['Sample as clicked'])
 
%%% Show detail of sample
 figure; 
 plot(Echo.GPSMatlabtime(sample1(1:300)),GPS_spline(sample1(1:300))-43.11,'r')
 hold on;
 plot(Echo.GPSMatlabtime(sample1(1:300)),-Echo.depth(sample1(1:300)),'g')
datetick
grid on;
xlabel('time')
ylabel('Echosounder (g) and interpolated GPS (r) values')
title(['Small detail of sample'])

%% Input demeaned GPS and echosounder values in to crosscorrelation function 
a=GPS_spline(sample1)-mean(GPS_spline(sample1));
b=Echo.depth(sample1)-mean(Echo.depth(sample1));

maxlag=30;      % Maximum time lag in sec
dt=0.1;         % sampling rate
subwindow=50;    % number of subseries

[c, lags] = jxcorr(a,b,maxlag,dt,subwindow);
[temp,ind]=max(c);      %% find max correlation peak
latency=lags(ind)
        
%%% plot crosscorrelation
figure
plot(lags, c)
axlimits=axis;
hold on
plot([latency latency], [axlimits(3:4)],'r') % line at maximum correlation
grid on;
xlabel('latency in seconds')
ylabel('Correlation value')
title(['maximum correlation for a latency of ' num2str(latency) ' sec'])

if (Echo.GPStime_s(sample1(end)) - Echo.GPStime_s(sample1(end)))/subwindow <= 10
    text(-15,0.021,'Warning, clicked sample very small')
end

%% 2. Analyse pc clock deviation (Transfer this to another file that analyses the results and the latency error)
figure; 
plot (GPS.Matlabtime,GPS.GPStime_s_midnight-GPS.PCtime_s_midnight - (GPS.GPStime_s_midnight(1)-GPS.PCtime_s_midnight(1)) )
datetick('x')
ylabel('pc clock deviation [s] since the start of hypack')
xlabel('GPS time');

%% End of file..functions below

%%
function [c, lags] = jxcorr(a,b,maxlag,dt,subwindow);

%[c lags]=jxcorr(a,b,maxlag,dt,subwindow)
%
%a,b       - input signals
%dt        - sampling time (in seconds)
%maxlag    - window of interest (in seconds)
%subwindow - # of subseries to average
%
%This program is uses xcorr.m  %% Signal Processing toolbox!
%Jamie MacMahan 3/14/02

if  length(a)<=1,
   error('1st input argument must be a vector or matrix.')
end
if nargin < 5, subwindow = 1; end
if nargin < 4, dt = 1; end
if nargin < 3, maxlag = length(a); end
if nargin < 2, b=a; end

subwindow=subwindow+1;
index=round(linspace(1,length(a),subwindow));
maxlag=round((maxlag-1)./dt);
for i=1:length(index)-1;
         AA=detrend(a(index(i):index(i+1)));
         BB=detrend(b(index(i):index(i+1)));
         [C,LAGS]=xcorr(AA,BB,maxlag,'biased');
     if i==1
         C1=zeros(size(C));
     end
     C1=C./sqrt(var(AA).*var(BB))+C1;
	 %C1=C+C1;
end
c=C1./(subwindow-1);
lags=LAGS.*dt;



