%% process_auvmission: to process individual AUV mission.
%% Series of missions (e.g. for a transect) are appended with assemble_section.m
%% after application of process_mission to each AUV mission.
%%
%% Without LBL navigation, latitude and longitude information is limited to the 
%% last recorded lat/lon before vehicle dive, and lat/lon are simply interpolated
%% between successive mission starting points during application of 
%% assemble_section.m to define the X and Y spatial attributes aligned with 
%% depth and the instrument measurements. 
%%
%% USAGE:
%% From within the AUV mission directory, issue process_auvmission command.
%%
%% In the current DAS, there are separate time stamps for [ctd], [lat lon], 
%% [roll pitch yaw], [depth].
%% dataCTD ~ 4 Hz (0.25 second)
%% dataGPS (lat/lon) ~ 1 Hz (1 second)
%% dataN (depth, roll, pitch, yaw) ~ 5 Hz (0.2 second)
%% Also, the nominal 4 Hz sampling rate is imperfect (variable), so proper 
%% alignment accounting for instrument plumbing and response lags requires that 
%% a common time reference be established for binning and interpolating all 
%% data to a common reference frame.  1 second averaging is chosen here.
%% This version will yield variables: temperature, salinity, density, oxygen, 
%% optical backscatter, and pH (the instrument suite on the Odyssey IIB AUV 
%% used for day operations during MUSE).
%%
%% MUSE 1 notes for Temperature and Salinity
%% Some temperatures require correction at the temp_frequency level;
%% corrected temperature is then used to compute salinity.
%% Also, the instruments were changed out on the morning of 8/28 before
%%  the day AUV ops.  New coefficients are required for this.
%% process_auvmission calls other m-files: 
%% bin_mission.m, get_spatial.m, missiondata.m, salsp.m, sw_dens.m, calc_O2.m, 
%% calc_pH.m, plotall.m, bearing.m, giffig_wide
%%
%% John Ryan 6/8/01 


%%----------------------
%% Load CTD data
%%----------------------
clear all; close all; 
load dataCTD; ctime = ctdlog(:,time); 	
t = ctdlog(:,temperature); c=ctdlog(:,conductivity); phv=ctdlog(:,pH); obs=ctdlog(:,OBS); f=ctdlog(:,fluor); 
oc=ctdlog(:,o2current); ot=ctdlog(:,o2temp); 

%%-----------------------------------
%% Define Instrument Lags (seconds) 
%%-----------------------------------
fluor_lag = 2;	% plumbing lag only
O2_lag = [3 2];	% plumbing lag and response lag 

%%------------------------
%% Load Spatial Data 
%% Note: yaw is heading
%%------------------------
% [glat,glon,gpstime,ndepth,nroll,npitch,nyaw,navtime,nposx,nposy] = get_spatial(1);
get_spatial
%% Lat and lon fixed at last recorded position before mission submergence
startlat = glat(length(glat)); startlon = glon(length(glon));
missiondata; 	% end point navigation 

%%-----------------------------------------------------------
%% Define standard 1 second average time frame from CTD data 
%%-----------------------------------------------------------
avgint = 1;	% averaging interval, seconds
mtime = [ceil(min(ctime)):avgint:floor(max(ctime))]';

%%---------------------------
%% Bin all data to 1 second
%%---------------------------
%% All variables will be averaged to 1 second, except lat/lon which will be interpolated to ctd time.
%% For February 2000 Odyssey 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
disp(['Time bin (seconds) is ' int2str(avgint) '.  Binning data ...']);
[T,C,PHV,OBS,F,OC,OT,depth,pitch,roll,yaw,Xpos,Ypos] = bin_mission(t,c,phv,obs,f,oc,ot,ndepth,npitch,nroll,nyaw,nposx,nposy,mtime,ctime,navtime);
%% Align fluorescence (oxygen aligned when calculated below)
fshift = round(fluor_lag/avgint);     % plumbing lag and shift
F=[F([fshift+1:length(F)]); repmat(NaN,fshift,1)];

disp(['Calculating derived variables...']);

%%--------------------------
%% Salinity and density
%%--------------------------
pres = sw_pres(depth,repmat(startlat,size(depth)));	% pressure [db] from depth and lat
S = salsp(C,T,pres); D = sw_dens(S,T,pres);

%%------------------------------
%% Oxygen (and alignment)
%%------------------------------
disp(['Computing Oxygen']);
%% Get coefficients
o2serial = 544; 	% oxygen sensor serial number
[m,b,k,c,SOc,BOc] = O2_calibration(o2serial);
%% Oxygen
%% Lag time estimate needed 
%% [Plumbing (+ chamber mixing) lag , response lag];  (seconds)
%% Sampling frequency taken into account in calc_O2.m
%% last input is either mll (ml/l) or umolkg (umol/kg)
% if strcmp(option,'control')==1; ounits='mll'; else ounits='umolkg';end
ounits='mll';
O2 = calc_O2(OC,OT,T,S,pres,m,b,k,c,SOc,BOc,mtime,O2_lag,ounits);

%%------------------------
%% pH
%%------------------------
% Calibration for Sea-Bird SBE-18 Serial Number 180312, 
% from NBS calibration 23, Feb 2000.
% off = 2.4738; slp = 4.6399; 
%% Note that pH sensor # 180312 was crushed early in experiment
%% Used # 180280
%% Calibration of 9/22/00
% off = 2.76; slp=4.53;
% off = 2.55; slp=1.24;
%% Factory calibration from etp3;
off = 2.4873; slp=4.5665;
PH = calc_pH(PHV, slp, off, T);

%% Plot
figure(1); clf; 
set(gcf,'color',[1 1 1],'units','pixels','position',[10 40 1200 800]);
plotall

%%------------------------
%% Cleanup Options
%%------------------------
showmenu=1;
mfsize=5;  % 1-D median filter size, for noisy variables (typically salinity)

while showmenu
A=menu('Choices','Crop time limit.','Excise a portion of a variable.','Filter a variable.','Start processing over.','Plot 3D path.','Save.','Keyboard.');

if A==1			% Crop based on time limits
 disp(['Select time subset']); figure(1); [x y]=ginput(2);
 idx = find(mtime>=x(1) & mtime<=x(2));
 T=T(idx); S=S(idx); PH=PH(idx); O2=O2(idx); depth=depth(idx); mtime=mtime(idx);
 F=F(idx); OBS=OBS(idx); D=D(idx);
 roll=roll(idx); pitch=pitch(idx); yaw=yaw(idx); Xpos=Xpos(idx); Ypos=Ypos(idx);
 plotall;
 
elseif A==2		% Excise portion of a variable	 
   varname=input('Which variable (T, S, PH, O2, F, OBS)?','s');
   eval(['wvar=' varname ';']);
   figure(2); clf; set(gcf,'color',[1 1 1],'units','pixels','position',[10 40 1200 800]);
   plot(wvar);
   B=1; 
   while B
      C=menu('Excise portion menu','Choose portion of record','done');
      if C==1
        figure(2); [x y]=ginput(2); x=round(sort(x));
        if x(1)<1; x(1)=1; end; if x(2)>length(wvar); x(2)=length(wvar); end
      	wvar(x(1):x(2))=NaN;
      	clf; plot(wvar);
      else
         figure(2); clf;
         B=0; eval([varname '=wvar;']); plotall
      end
   end
   
elseif A==3		% Filter a variable
   varname=input('Which variable (T, S, PH, O2, F, OBS)?','s');
   eval(['wvarO=' varname ';']); %% retain original
   wvarf=wvarO; 
   figure(2); clf; set(gcf,'color',[1 1 1],'units','pixels','position',[10 40 1200 800]);
   plot(wvarO);
	B=1; 
   while B
      C=menu('Filter Options','median filter','restore original','save change');
      if C==1	% median filter
         wvarf=medfilt1(wvarf,mfsize);      	
         figure(2); clf; plot(wvarO,'k'); hold on; plot(wvarf,'m');
      elseif C==2;
         wvarf=wvarO;
         figure(2); clf; plot(wvarf,'k');
      else
         figure(2); clf;
         B=0; eval([varname '=wvarf;']); plotall
      end
   end
		           
elseif A==4; 	% Start over
		process_auvmission
      
elseif A==5; 	% Plot path
   figure(3); clf; comet3(Xpos,Ypos,-1*depth); title('AUV path');
   
elseif A==6		% Save
 disp(['Saving to file mission_sd.mat']);
 %% Clear density where T/S has been deleted
 idx=find(isnan(T) | isnan(S)); D(idx)=NaN;
 save mission_sd T S D PH O2 F OBS mtime depth startlat startlon roll pitch yaw Xpos Ypos missionID ep 
 disp(['Plotting processed mission data...']);
 showmenu=0;
 plotall % check_mission
 figure(2); clf; plot(Xpos,Ypos); title('Dead reckoning Position'); ylabel('Y'); xlabel('X');
else				% Keyboard Control
 keyboard
end

end

