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;
outfile = sprintf('%s_lag',in_file);
fid=fopen(outfile,'wt');

%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;

%define the periods that ship data should be compared to inst data
shift=[datenum('00:00:10'):datenum('00:00:05'):datenum('00:10: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;

%define normal increments
norm_inc=datenum('00:00:10');
[y,m,d,h,mi,s]=datevec(norm_inc);
s=round(s);
norm_inc=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:50: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 x minute windows
window=(now_time-win_add);%x minute data window

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

%loop through the entire merged file
while end_time>=now_time
    %normal pco2 time stamps
    time_vec=[window:norm_inc:now_time];
    
    %grab an instrument time window as defined above by window and now_time
    for i=1:length(time_vec)
    inst_time_window(i)=near(data(:,1),time_vec(i),1);
    end;
    
    %for each inst window grab a ship time window for every period defined by shift above
    for j=1:length(shift)
        
        for i=1:length(time_vec)
        ship_offset_window(i)=near(data(:,1),(time_vec(i)-shift(j)),1);
        end;
        %find the difference in temperatures
        temp_offset=data(inst_time_window,2)-data(ship_offset_window,4);
        %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);
    clear inst_time_window;

    
        
end;

fprintf(fid,'%Time Time_Lag SDN_Time SDN_Time_Lag\n');
for i=1:length(time_store)
fprintf(fid,'%s %s %f %f\n',datestr(time_store(i)),datestr(lag_store(i),13),time_store(i),lag_store(i));
end;

fclose(fid);
open(outfile);
