%% clean_mission: to assemble and clean AUV mission data.
%% choose to focus on either CTD or AUV control

clear all; addpath 'review/' '/u/ryjo/matlab/CSIRO/'
option = 'control'

%%------------------------
%% Define path for review
%%------------------------
mission = input('Enter mission: ','s'); eval(['load ' mission '/dataCTD']);
mtime = ctdlog(:,time); 

%%------------------------
%% Temperature 
%%------------------------
t = ctdlog(:,temperature);

%%------------------------
%% Spatial 
%%------------------------
%% Get spatial data and vehicle orientation.    Note:
%%  yaw is heading
%%  separate time stamps for [ctd], [lat lon], [roll pitch yaw], [depth]
[lat,lon,gpstime,depth,roll,pitch,yaw,navtime] = get_spatial(mission);
%% Lat and lon fixed at last recorded position before mission submergence
startlat = lat(length(lat)); startlon = lon(length(lon));
%% Interpolate depth to ctdtime
cdepth = interp1(navtime,depth,mtime); clear depth
croll = interp1(navtime,roll,mtime); clear roll;
cpitch = interp1(navtime,pitch,mtime); clear pitch;
cyaw = interp1(navtime,yaw,mtime); clear yaw;
%% Create lat vector for pressure calculation
clat = repmat(mean(lat),length(ctdlog),1);
%% Compute pressure (db) from depth
pres = sw_pres(cdepth,clat);		% [db]

%%------------------------
%% Time check 
%%------------------------
% minCTDtime = min(mtime); maxCTDtime = max(mtime);
% minGPStime = min(lon(:,1)); maxGPStime = max(lon(:,1));
% minNAVtime = min(roll(:,1)); maxNAVtime = max(roll(:,1));
% minNAVtime = min(lon(:,1)); maxNAVtime = max(lon(:,1));
% if (minCTDtime-minNAVtime)>0; disp(['CTD time begins after NAV time']);
% else disp(['CTD time begins before NAV time']); end
% if (maxCTDtime-maxNAVtime)>0; disp(['CTD time ends after NAV time']);
% else; disp(['CTD time ends before NAV time']); end
% disp([' ']);

%%------------------------
%% Salinity
%%------------------------
disp(['Computing Salinity']); [s] = salsp(ctdlog(:,conductivity),t,pres);

%%------------------------
%% Oxygen
%%------------------------
disp(['Computing Oxygen']);
%% oxygen current and oxygen membrane temperature 
oc = ctdlog(:,7); ot = ctdlog(:,6); 
%% Get coefficients
sensor = '368'; prepost='post';
[m,b,k,c,SOc,BOc] = O2_calibration(sensor,prepost);
%% Oxygen
%% Lag time estimate needed 
Olag = [2.6 0];	% Oxygen lags, first is plumbing; 
		% second is sensor time constant
disp(['Lags defined for oxygen are ' num2str(Olag) ' seconds']);
%% last input is either mll (ml/l) or umolkg (umol/kg)
if strcmp(option,'control')==1; ounits='mll'; else ounits='umolkg';end
o2 = calc_O2(oc,ot,t,s,pres,m,b,k,c,SOc,BOc,mtime,Olag,ounits);
idx=find(o2 <0 | o2 >10); o2(idx)=NaN;

%%------------------------
%% pH
%%------------------------
% Calibration for Sea-Bird SBE-18 Serial Number 180312, 
%  from NBS calibration 23, Feb 2000.
pHV = ctdlog(:,8); off = 2.4738; slp = 4.6399;
ph = calc_pH(pHV, slp, off, ctdlog(:,temperature));


figure(1); clf; paper;
dZdt=[NaN;diff(pres)./diff(mtime)]; up=find(dZdt<-.1); down=find(dZdt>.1);

%% Bin O2
o2up=o2(up); o2down=o2(down); Zup=pres(up); Zdown=pres(down);
Zbin = [15:10:max(pres)-5]; buf=5;
for N = 1:length(Zbin);
 uidx=find(Zup>=Zbin(N)-buf & Zup<Zbin(N)+buf);
 didx=find(Zdown>=Zbin(N)-buf & Zdown<Zbin(N)+buf);
 o2u(N)=mean(o2up(uidx)); o2d(N)=mean(o2down(didx));
end

subplot(311); 
p=plot(t(up),pres(up),'r.'); set(p,'markersize',4); hold on;
p=plot(t(down),pres(down),'b.'); set(p,'markersize',4); 
set(gca,'ydir','reverse','xlim',[min(t) max(t)],'ylim',[0 max(pres)]);
xlabel('Temp'); ylabel('pres');

subplot(312); 
p=plot(s(up),pres(up),'r.'); set(p,'markersize',4); hold on;
p=plot(s(down),pres(down),'b.'); set(p,'markersize',4); 
set(gca,'ydir','reverse','xlim',[33.2 34.2],'ylim',[0 max(pres)]);
xlabel('Salinity'); ylabel('pres');

subplot(313); 
% p=plot(o2(up),pres(up),'r.'); set(p,'markersize',4); hold on;
plot(o2u,Zbin,'r'); hold on;
% p=plot(o2(down),pres(down),'b.'); set(p,'markersize',4); 
plot(o2d,Zbin,'b');
set(gca,'ydir','reverse','xlim',[1 6],'ylim',[0 max(pres)]);
xlabel('O_2'); ylabel('pres');

mndif = mean(o2d-o2u);
figure(2); clf; giffig;
subplot(211); plot(o2d-o2u,Zbin);
title(['Oxygen difference, down-up; mean: ' num2str(mndif)]);
set(gca,'ydir','reverse');
subplot(212); plot(t,o2,'*')
idx = find(finite(o2));
c=corrcoef(t(idx),o2(idx)); c=c(2);
title(['Correlation at tlag : ' num2str(sum(Olag)) ' seconds = ' num2str(c)]);
xlabel('temperature'); ylabel('o2');


return

%%------------------------
%% Average 
%%------------------------
%% For February AUV missions, ctd horizontal resolution was ~ 3.5 measurements 
%%  per along-track meter, with a temporal resolution of .25 seconds.
%%  This translates to a vehicle speed of ~ 1 m/s
avgint = 2.5;	% averaging interval, seconds, matches O2 offset
[T,S,pH,O2,time,depth,epochsec,pitch,roll,yaw] = avgCTD(t,s,ph,o2,mtime,cdepth,cpitch,croll,cyaw,avgint);

%%------------------------
%% Plot 
%%------------------------
sax = [33.2 34.2]; pHax = [7.5 8.05]; O2ax = [0 5];
XL = 'mission time';

figure(1); clf; 
% paper;
homefig;

subplot(511); plot(time,depth);
xlabel(XL); ylabel('depth'); set(gca,'Ydir','reverse');

subplot(512); plot(time,T,'r'); xlabel(XL); ylabel('temperature');
set(gca,'Ycolor','r','box','off','YaxisLocation','left');
P = get(gca,'position'); axes('position',P);
plot(time,S,'b'); set(gca,'YaxisLocation','right','Ycolor','b','box','off');
set(gca,'Ylim',sax); ylabel('salinity'); set(gca,'color','none');

subplot(513); plot(time,pH); ylabel('pH'); hold on; set(gca,'Ylim',pHax); 

subplot(514); plot(time,O2); ylabel('Oxygen (ml/l)'); % set(gca,'Ylim',O2ax); 
disp([' ']); disp(['Correct for oxygen time offset !!!']);
disp(['How much difference does use of ambient temperature make?']); disp([' ']);

subplot(515); plot(time,pitch,'k',time,roll,'b',time,yaw,'r');
% legend('pitch','roll','yaw');

%%------------------------
%% Subset for ctd or not for control 
%%------------------------
if strcmp(option,'ctd')==1;
 disp(['Select time subset']); [x y]=ginput(2);
 idx = find(time>=x(1) & time<=x(2));
 T=T(idx); S=S(idx); pH=pH(idx); O2=O2(idx); depth=depth(idx); time=time(idx);
 epochsec=epochsec(idx); roll=roll(idx); pitch=pitch(idx); yaw=yaw(idx);
 disp(['Saving to file ctd.mat']);
 eval(['save ' mission '/ctd T S pH O2 epochsec depth startlat startlon'])
else
 disp(['Saving to file auv.mat']);
 eval(['save ' mission '/auv T S pH O2 epochsec depth startlat startlon roll pitch yaw'])

end
