%This m file helps determine the time offset between the inlet to the ship and the measurement
%at the licor system. It can be used to create a lookup table to determine change in temperature
%that occurs in the ships plumbing and effects pCO2. The temperature offset of a period of 
%instrument data is examined and compared to different windows of ship data. The point at
%which the standard deviation of the instrument data subtracted from the ships data is minimized
%is returned as the time offset for that period. Modified from code written by Pete Strutton.

data=[];
for i = [10:11]     %These are Julian days - Jan 1 = 1 (Duh!). Change this for the days you want to process.
    
    in_file = sprintf('pco2_2002%03d_merge.txt', i);     %open the merged pCO2 input file
    fprintf('starting julian day %d\n', i);
    pco2_num = load(in_file);
    data=[data; pco2_num(:,[2 10 15 24])];%get inst_time,inst_temp,ship_time, ship_temp
    clear pco2_num ;
end;

%round off the in time so that its easier to use
[y,m,d,h,mi,s]=datevec(data(:,1));
s=round(s);
data(:,1)=datenum(y,m,d,h,mi,s);
clear y;clear m;clear d;clear h;clear mi;clear s;


shift=[datenum('00:00:10'):datenum('00:00:10'):datenum('00:12:00')];
[y,m,d,h,mi,s]=datevec(shift);
s=round(s);
shift=datenum(y,m,d,h,mi,s);
clear y;clear m;clear d;clear h;clear mi;clear s;

%for those non perfect measurement times 9.2138 seconds for example
[y,m,d,h,mi,s]=datevec(datenum('00:00:04'));
s=round(s);
buf=datenum(y,m,d,h,mi,s);
clear y;clear m;clear d;clear h;clear mi;clear s;

%how often to look for signatures
[y,m,d,h,mi,s]=datevec(datenum('01:00:00'));
s=round(s);
now_add=datenum(y,m,d,h,mi,s);
clear y;clear m;clear d;clear h;clear mi;clear s;

%the amount of time over which to look for signatures
[y,m,d,h,mi,s]=datevec(datenum('00:45:00'));
s=round(s);
win_add=datenum(y,m,d,h,mi,s);
clear y;clear m;clear d;clear h;clear mi;clear s;

%set some useful vars
end_time=max(data(:,1));%last time period
now_time=data(1,1)+now_add;%start of time and 20 minute windows
window=(now_time-win_add);%10 minute data window

%declare variables
dev_temp=[];
time_shift=[];
lag_store=[];
time_store=[];

%loop through the entire merged file
while end_time>=now_time
    
    %grab an instrument time window as defined above by window and now_time
    inst_time_window=find(data(:,1)>=(window-buf)&data(:,1)<=((now_time)+buf));
    
    %for each inst window grab a ship time window every period defined by shift above
    for j=1:length(shift)
        ship_offset_window=find(data(:,1)>=(window-(shift(j)+buf))&data(:,1)<=(now_time-(shift(j)-buf)));
        
        %in case of missing times such as 10-Jan-2002 18:46:00 & 18:46:20 where there are no values
        if length(ship_offset_window)~=length(inst_time_window)
            clear ship_offset_window;
            continue;
        end;
            
        %find the difference in temperatures
        temp_offset=data(inst_time_window,2)-data(ship_offset_window,4);%find the difference between the temps
        %find the standard dev in those temperatures & store it
        dev_temp=[dev_temp;std(temp_offset)]
        %store that time_shift
        time_shift=[time_shift;shift(j)];
        %next shift
        clear temp_offset;clear ship_offset_window;%again
    end;
    
    %find in what shift the standard deviation is smallest
    dev_ind=find(dev_temp==min(dev_temp));
    %store shift with min deviation
    lag_store=[lag_store;time_shift(dev_ind)];
    %what time of day was that
    time_store=[time_store;now_time];
    %next instrument window
    dev_temp=[];time_shift=[];
    now_time=now_time+now_add;
    window=(now_time-win_add);

    
        
end;
