% Process SBE 56 Thermistor Data
% this file takes .cnv files, with 
% colum 1 = time in julian days
% column 2 = temperature
run presentation_plot
clc, close all, clear all;

%% load data files
start_year=2012;
start_time=datenum(2012,9,19,13,0,0);
end_time=datenum(2012,9,25,17,0,0);
files=dir('SBE056*');

for i=1:size(files,1); % get order for files in increasing mab
filename=files(i).name;
temp = sscanf(filename,'SBE056%f _%f-%f-%f.cnv',[1 Inf]);
sn_v(i)=temp(1);
end
[mab,ii]=sort(sn_v,2);


%% load and process data

for i=1:size(files,1);
filename=files(ii(i)).name;
fid=fopen(filename);
DATA = textscan(fid,'%f %f %f','delimiter',',','headerlines',88);

temp_raw{i}=DATA{1,2};
mtime_raw{i}=jday2matlab(DATA{1,1},start_year); % convert julian dats to mtime

time_id=find(mtime_raw{i}>=start_time & mtime_raw{i}<=end_time);
temp_proc{i}=temp_raw{i}(time_id);
mtime_proc{i}=mtime_raw{i}(time_id);

% low pass filter temperature data
t_window = 10*60/3600/24; % averaging 10 min window
dt = mtime_proc{1}(2)-mtime_proc{1}(1); % data time step
window = 2.*round(t_window/dt/2);

tempT = 0.*temp_proc{i};
    if window<=1 % if filter is less than timestep, use orig data
        tempT=temp_proc{i};
    else % filter data otherwise
        for jj = window/2+1:size(tempT,1)-window/2-1;
            tempT(jj,:) = mean(temp_proc{i}(jj-window/2:jj+window/2,:),1);
        end
    end
temp_filt{i} = tempT;
clear tempT
i
end
%%
% interpolate temp data to a common time
dt=10*60*(1/86400); % interpolate to 10 min window
timevec=start_time:dt:end_time;
Temp_mat=nan(size(files,1),length(timevec));

for i=1:size(mab,2)
    Temp_mat(i,:)=interp1(mtime_proc{i}, temp_proc{i},timevec,'linear','extrap');
    Temp_filt(i,:)=interp1(mtime_proc{i}, temp_filt{i},timevec,'linear','extrap');
end

clear temp fid filename files mab_temp DATA i ii ff mab_v time_id temp_common;

save Processed_SBE56_data

%% plot data
clc, clear all, close all;
load Processed_SBE56_data

p=figure(1);
for i=1:size(mab,2)
plot(mtime_raw{i},temp_raw{i}); hold all;
leg{i}=sprintf('SN%g',mab(i));
datetickzoom('x');
end
legend((leg),'location','nw')
xlabel('Date'), ylabel('Temperature ($^{\circ}C$)');
title('Figure 1.b raw data');
hgsave(p,'figure1.fig');

p=figure(2);
for i=1:size(mab,2)
plot(mtime_proc{i},temp_proc{i}); hold all;
leg{i}=sprintf('SN%g',mab(i));
datetickzoom('x');
end
legend((leg),'location','nw')
xlabel('Date'), ylabel('Temperature ($^{\circ}C$)');
title('Figure 1.c truncated data');
hgsave(p,'figure2.fig');

p=figure(3);
for i=1:size(mab,2)
plot(mtime_proc{i},temp_filt{i}); hold all;
leg{i}=sprintf('SN%g',mab(i));
datetickzoom('x');
end
legend((leg),'location','nw')
xlabel('Date'), ylabel('Temperature ($^{\circ}C$)');
title('Figure 1.d filtered data');
hgsave(p,'figure2.fig');