%% process_auvmissionNST: process individual AUV mission
%%
%%
%% USAGE:
%% From within the AUV mission directory, issue process_auvmissionNST command.
%%
%% In the current DAS, there are separate time stamps for [ctd], [lat lon], 
%% [roll pitch yaw], [pressure depth].
%% dataCTD ~ 4 Hz (0.25 second)
%% dataGPS (lat/lon) ~ 1 Hz (1 second)
%% dataN (roll, pitch, yaw) ~ 5 Hz (0.2 second)
%% dataPS (pressure, depth) ~ 3 to 5 Hz
%% 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.
%%
%% It should be noted that temporal and spatial lags associated with the pumped data streams
%% have NOT been accommodated in this version of the code.  BEWARE
%%
%% This version will yield variables: temperature (deg C IPTS-68), salinity (psu PSS-78),
%% density (kg/m**3), oxygen(mL/L), optical backscatter (raw volts), and nitrate (uM)
%% (the instrument suite on the ALTEX AUV used for day operations during the 
%% Arctic October 2001 test cruise).
%%
%% ALTEX October 2001 cruise notes
%% 
%%
%% process_auvmission calls other m-files: 
%% bin_missionNST.m, get_spatialNST.m, missiondata.m, sw_dpth.m, sw_dens.m, 
%% sw_salt.m, sw_c3515.m, sw_cndr.m, O2_cal_SBE43.m, calc_O2_SBE43.m
%% Pres_cal.m, Temp_cal.m, Cond_cal.m
%% plotallNST.m, giffig_wide, NegTrunc.m
%% ncoutn2.m, ncout1.m, ncout2.m, ncoutbinned.m, ncoutips.m
%% 
%%
%% John Ryan 6/8/01 
%% Nicole Tervalon  Sep 2001


clear all; close all;

%%---------------------------------------------
%% Vehicle configuration and other user inputs
%%---------------------------------------------
% Length of vehicle midbody between tail section joining ring and science payload
% joining ring
midlength = 53;  % inches
% Distance from center of MVC sphere to outside of tail section joining ring
taillength = 18;  % inches
% Distance from center of conductivity and temperature sensors to outside of nose
% section joining ring
noselength = 36;  % inches
% Calibrate pressure sensor or use calibrated AUV output?  1 = calibrate it, 0 = use AUV output
calpress = 1;
% Calculate sound velocity of water column from CTD data?  1 = calculate it, 0 = estimate it
calcsoundvel = 0;
% Start year for time reference (generally the year in which the data was collected)
startyear = 2001;
% Start month for time reference (generally the month in which the data was collected)
startmonth = 8;
% Start day for time reference (generally a day prior to the day the data was collected)
% For ALTEX mission leave this day set to 1
startday = 1;

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

%%---------------------------------
%% Other user defined coefficients
%%---------------------------------
% Nominal latitude of mission
nomLAT = 80.0;  % (deg North)
% nomLAT = startlat;


%%-----------------------
%% Load Pressure Data
%%-----------------------
if ( exist('dataPS') == 2 )
   disp(' Loading Pressure data ')
   dataPS;
  % truncate pressure data
   parosci(:,time) = parosci(:,time) - startTime;
   parosci = NegTrunc(parosci);
  % filter parosci data for bad data hits using the following methods
  %  set bad hits to previous pressure value
   parosci = parosci(find(diff(parosci(:,pres_period)) > -0.01),:);
   ptime = parosci(:,time);
   pres_serialnum = 77988;
   % pres_serialnum = PresSerialNumber;
   if calpress == 1
    % calibrate pressure sensor
      [pu0, py1, py2, py3, pc1, pc2, pc3, pd1, pd2, pt1, pt2, pt3, pt4, pt5] = ...
           Pres_cal(pres_serialnum);
       pu = parosci(:,temp_period)-pu0;
       pc = pc1 + pc2*pu + pc3*(pu.^2);
       pd = pd1 + pd2*pu;
       pt0 = pt1 + pt2*pu + pt3*(pu.^2) + pt4*(pu.^3) + pt5*(pu.^4);
       pres = pc.*(1-(pt0.^2)./(parosci(:,pres_period).^2)).* ...
           (1-pd.*(1-(pt0.^2)./(parosci(:,pres_period).^2)));  % (in psi)
       pres = pres * 0.06894757;  % (in bars)
       clear pu pc pd pt0;
   else
      pres = parosci(:,pressure);
   end
   atmpres = median(pres(1:20));  % (bar)
% If initial location of vehicle is not at the water surface, set the atmospheric pressure here
%   atmpres = 1.01325;  % (bars)
% Set pressure measurement to RELATIVE pressure, NOT ABSOLUTE
   pres = pres - atmpres;
   clear parosci;
else
   disp(' No pressure data file found - exiting ')
   return
end



disp (' ');
%%----------------------
%% Load CTD data
%%----------------------
%% Look for PRIMARY CTD data
if ( exist('dataCTD') == 2 )
   disp (' Loading primary CTD data ')
   dataCTD;
   % FILTER DATA FOR BAD TIME STAMPING (for data collected prior to 2001, yearday 235)
      disp ('   FILTERING OLD DATA FOR BAD TIMESTAMPS ')
      diffcondtemp = (diff(ctdDriver(:,conductivity)) | diff(ctdDriver(:,temperature)));
      ctdlogtemp = ctdDriver(find(diffcondtemp == 1),:);
      ctdlog = ctdlogtemp(find(diff(ctdlogtemp(:,time)) ~= 0),:);
   % OR use the following definition
   % ctdlog = ctdDriver;
   clear ctdDriver diffcondtemp ctdlogtemp;
   temp1_serialnum = TempSerialNumber;
   cond1_serialnum = CondSerialNumber;
% truncate CTD data if it started recording before the navigation data did
   ctdlog(:,time) = ctdlog(:,time) - startTime;
   ctdlog = NegTrunc(ctdlog);
   c1time = ctdlog(:,time);
% determine SBE25 boardset measured year and decimal day from time variable
   [c1year,c1month,c1day,c1hour,c1min,c1sec] = datevec((ctdlog(:,time)+startTime)/86400 + ...
    datenum(1970,1,1,0,0,0));
   c1abstime = (datenum(c1year,c1month,c1day,c1hour,c1min,c1sec)-datenum(startyear,startmonth,startday,0,0,0));
   c1yearday = datenum(c1year,c1month,c1day,c1hour,c1min,c1sec) - datenum(c1year,1,1,0,0,0) + 1;
   clear c1month c1day c1hour c1min c1sec;
% rename variables (use the commented out line if looking at data collected prior to 4 August 2001
   t1freq = ctdlog(:,temp_frequency);  % (Hz)
   c1freq = ctdlog(:,cond_frequency)./1000;  % (kHz)
   obs = ctdlog(:,OBS); o2v = ctdlog(:,dissolvedO2);
% calibrate flow data using current estimates (to L/min)
   flw1 = (exp((ctdlog(:,flow1)+4.7685)/1.4816)-3.92)/111.5;
   flw1(find(ctdlog(:,flow1) == 0)) = 0;
   % obs=ctdlog(:,dissolvedO2); o2v=ctdlog(:,OBS); flw1 = ctdlog(:,flow1);
% calibrate temperature (deg C)
   disp(' Calibrating Temperature ');
   [t1a,t1b,t1c,t1d,t1f0] = Temp_cal(temp1_serialnum);
   t1 = 1./(t1a + t1b*(log(t1f0./t1freq)) + t1c*((log(t1f0./t1freq)).^2) + ...
           t1d*((log(t1f0./t1freq)).^3)) - 273.15;
   clear t1freq;
% calculate/interpolate pressure at CT sensors for each data point and account for pitch of vehicle
% (pressure in bars, pitch in degrees)
   buffer = median(diff(c1time))*10;
   for N=1:length(c1time)
      idxii = find(ptime >= c1time(N)-buffer & ptime < c1time(N)+buffer);
      pres1(N)=interp1(ptime(idxii),pres(idxii),c1time(N));
      idxi = find(navtime >= c1time(N)-buffer & navtime < c1time(N)+buffer);
      pitch1(N)=interp1(navtime(idxi),npitch(idxi),c1time(N));
   end
   pres1 = pres1'; pitch1 = pitch1';
% calibrate conductivity (Siemens/meter) [using uncorrected pressure because conductivity
% dependence on pressure is very small (on the order of 10^-7)]
   disp(' Calibrating Conductivity ');
   [c1a,c1b,c1c,c1d,c1m,c1Cor] = Cond_cal(cond1_serialnum);
   c1 = (c1a*(c1freq.^c1m)+c1b*(c1freq.^2)+c1c+c1d*t1)./(10*(1+c1Cor*10*pres1));
   clear c1freq;
   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   poffset = 10.5 * cos(pitch1*pi/180);
   % Vertical distance from center of MVC sphere to CT sensors (inches) (neg when nose down)
   MVCoffset = (taillength + midlength + noselength) * sin(pitch1*pi/180);
   % Total vertical distance between pressure sensor and CT sensors (meters) (pos when nose down)
   totdeltad1 = (- MVCoffset - sign(pitch1).*poffset)/39.37; 
   % Calculate approximate saltwater density (in order to calculate delta pressure)
   c3515 = sw_c3515;  % (mS/cm)
   cratio = c1*10/sw_c3515;
   s = sw_salt(cratio,t1,pres1*10); % (psu)
   dens = sw_dens(s,t1,pres1*10); % (kg/m^3)
   % Pressure at CT sensors based on parosci pressure and vertical offset
   p1 = pres1 + (dens * 9.81 .* totdeltad1)/100000; % (bars)
% calculate actual salinity and density
   s1 = sw_salt(cratio,t1,p1*10);  % (psu)
   d1 = sw_dens(s1,t1,p1*10);  % (kg/m^3)
% calculate depth (meters) and set gauge pressure to 0 if it is negative
%   p1gauge = p1-atmpres;
%   p1gauge(find(p1gauge<0)) = 0;
%   dpth1 = sw_dpth(p1gauge,nomLAT);
   p1(find(p1<0)) = 0;
   dpth1 = sw_dpth(p1*10,nomLAT);
% calibrate oxygen
   disp(' Computing Oxygen ');
   O2_lag = [3 2];	% plumbing lag and response lag 
   % Get coefficients
   o2serial = 103; 	% oxygen sensor serial number
   [O2_offset, Tcor, Pcor, SOc, BOc] = O2_cal_SBE43(o2serial);
   % if strcmp(option,'control')==1; ounits='mll'; else ounits='umolkg';end
   ounits='mll';
   o2 = calc_O2_SBE43(o2v,t1,s1,p1*10,O2_offset,Tcor,Pcor,SOc,BOc,c1time,O2_lag,ounits);
% clear data array  
   clear ctdlog;
% print calibrated data to netcdf file
   ncout1;
% check for bad ctd dacq
   c1timestep = diff(c1time);
   figure; clf; 
     plot(c1timestep) 
     xlabel('datapoint')
     ylabel('Primary Boardset Timestep (s)')
     title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
     grid
else
   disp (' Primary CTD log file not found ')
end

disp (' ');
%% Look for SECONDARY CTD data
if ( exist('dataCTD2') == 2 )
   disp (' Loading secondary CTD data ')
   dataCTD2;
   % FILTER DATA FOR BAD TIME STAMPING (for data collected prior to 2001, yearday 235)
      disp ('   FILTERING OLD DATA FOR BAD TIMESTAMPS ')
      diffcondtemp = (diff(ctdDriver2(:,conductivity)) | diff(ctdDriver2(:,temperature)));
      ctdlogtemp = ctdDriver2(find(diffcondtemp == 1),:);
      ctdlog = ctdlogtemp(find(diff(ctdlogtemp(:,time)) ~= 0),:);
   % OR use the following definition
   % ctdlog = ctdDriver2;
   clear ctdDriver2 diffcondtemp ctdlogtemp;
   temp2_serialnum = TempSerialNumber;
   cond2_serialnum = CondSerialNumber;
% truncate CTD data if it started recording before the navigation data did
   ctdlog(:,time) = ctdlog(:,time) - startTime;
   ctdlog = NegTrunc(ctdlog);
   c2time = ctdlog(:,time);
% determine SBE25 boardset measured year and decimal day from time variable
   [c2year,c2month,c2day,c2hour,c2min,c2sec] = datevec((ctdlog(:,time)+startTime)/86400 + ...
    datenum(1970,1,1,0,0,0));
   c2abstime = (datenum(c2year,c2month,c2day,c2hour,c2min,c2sec)-datenum(startyear,startmonth,startday,0,0,0));
   c2yearday = datenum(c2year,c2month,c2day,c2hour,c2min,c2sec) - datenum(c2year,1,1,0,0,0) + 1;
   clear c2month c2day c2hour c2min c2sec;
% rename variables
   t2freq = ctdlog(:,temp_frequency);  % (Hz)
   c2freq = ctdlog(:,cond_frequency)./1000;  % (kHz)
% calibrate flow data using current estimates (to L/min)
   flw2 = (exp((ctdlog(:,flow2)+4.7685)/1.4816)-3.92)/111.5;
   flw2(find(ctdlog(:,flow2) == 0)) = 0;
% calibrate temperature (deg C)
   disp(' Calibrating Temperature ');
   [t2a,t2b,t2c,t2d,t2f0] = Temp_cal(temp2_serialnum);
   t2 = 1./(t2a + t2b*(log(t2f0./t2freq)) + t2c*((log(t2f0./t2freq)).^2) + ...
           t2d*((log(t2f0./t2freq)).^3)) - 273.15;
   clear t2freq;
% calculate/interpolate pressure at CT sensors for each data point and account for pitch of vehicle
% (pressure in bars, pitch in degrees)
   buffer = median(diff(c2time))*10;
   for N=1:length(c2time)
      idxii = find(ptime >= c2time(N)-buffer & ptime < c2time(N)+buffer);
      pres2(N)=interp1(ptime(idxii),pres(idxii),c2time(N));
      idxi = find(navtime >= c2time(N)-buffer & navtime < c2time(N)+buffer);
      pitch2(N)=interp1(navtime(idxi),npitch(idxi),c2time(N));
   end
   pres2 = pres2'; pitch2 = pitch2';
% calibrate conductivity (Siemens/meter) [using uncorrected pressure because conductivity
% dependence on pressure is very small (on the order of 10^-7)]
   disp(' Calibrating Conductivity ');
   [c2a,c2b,c2c,c2d,c2m,c2Cor] = Cond_cal(cond2_serialnum);
   c2 = (c2a*(c2freq.^c2m)+c2b*(c2freq.^2)+c2c+c2d*t2)./(10*(1+c2Cor*10*pres2));
   clear c2freq;
   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   poffset = 10.5 * cos(pitch2*pi/180);
   % Vertical distance from center of MVC sphere to CT sensors (inches) (neg when nose down)
   MVCoffset = (taillength + midlength + noselength) * sin(pitch2*pi/180);
   % Total vertical distance between pressure sensor and CT sensors (meters) (pos when nose down)
   totdeltad2 = (- MVCoffset - sign(pitch2).*poffset)/39.37; 
   % Calculate approximate saltwater density (in order to calculate delta pressure)
   c3515 = sw_c3515;  % (mS/cm)
   cratio = c2*10/sw_c3515;
   s = sw_salt(cratio,t2,pres2*10); % (psu)
   dens = sw_dens(s,t2,pres2*10); % (kg/m^3)
   % Pressure at CT sensors based on parosci pressure and vertical offset
   p2 = pres2 + (dens * 9.81 .* totdeltad2)/100000; % (bars)
% calculate actual salinity and density
   s2 = sw_salt(cratio,t2,p2*10);  % (psu)
   d2 = sw_dens(s2,t2,p2*10);  % (kg/m^3)
% calculate depth (meters) and set gauge pressure to 0 if it is negative
   p2(find(p2<0)) = 0;
   dpth2 = sw_dpth(p2*10,nomLAT);
% clear data array  
   clear ctdlog;
% print calibrated data to netcdf file
   ncout2;
% check for bad ctd dacq
   c2timestep = diff(c2time);
   figure; clf; 
     plot(c2timestep) 
     xlabel('datapoint')
     ylabel('Secondary Boardset Timestep (s)')
     title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
     grid
else
   disp (' Secondary CTD log file not found ')
end

if ((exist('dataCTD') == 0) & (exist('dataCTD2') == 0))
    disp (' No CTD data found - exiting ')
    return
end

% Update the matrix name depending on which driver is being used
% if ( exist('ctdDriverctdlog') == 1 )
%   ctdlog = ctdDriverctdlog;
% else ( exist('seabirdctdlog') == 1 )
%   ctdlog = seabirdctdlog;
% end


disp (' ');
%%-----------------------------------
%% Load ISUS data if it exists
%%-----------------------------------
if ( exist('dataIsus') == 2 )
   disp (' Loading ISUS data, matrix = isuslog ')
   dataIsus
   %% filter ISUS data (look for n2qual <= some value)
   nqualcut = 0.005;
   isuslog = isuslog(find(isuslog(:,isusQuality) <= nqualcut),:);
   % truncate ISUS data if it started recording before the navigation data did
   isuslog(:,1) = isuslog(:,1) - startTime;
   isuslog = NegTrunc(isuslog);

   n2 = isuslog(:,isusNitrate);
   n2qual = isuslog(:,isusQuality);
   n2time = isuslog(:,time);
   clear isuslog;
   
   % check for bad ISUS dacq
   itimestep = diff(n2time);

   % calculate/interpolate pressure at n2 sensor for each data point and account for pitch of vehicle
   buf = median(diff(n2time))/2;
   if ( length(n2time) >= 2 )
      for N=1:length(n2time)
         idxi = find(navtime >= n2time(N)-buf & navtime < n2time(N)+buf);
        % N2(N)=median(n2(idxi));
         pitchn2(N)=interp1(navtime(idxi),npitch(idxi),n2time(N));
         n2depth(N)=interp1(navtime(idxi),ndepth(idxi),n2time(N));
      end
   end
   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   poffset = 10.5 * cos(pitchn2*pi/180);
   % Vertical distance from center of MVC sphere to N2 sensor (inches) (neg when nose down)
   MVCoffset = (taillength + midlength + noselength) * sin(pitchn2*pi/180);
   % Total vertical distance between pressure sensor and N2 sensor (meters) (pos when nose down)
   totdeltadn2 = (- MVCoffset - sign(pitchn2).*poffset)/39.37; 

   n2depth = n2depth + totdeltadn2;
    
   figure;
   plot(itimestep,'b')
   xlabel('datapoint')
   ylabel('ISUS Timestep (s)')
   title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
   grid
% determine n2 timeset year and decimal day from n2time variable
   [n2year,n2month,n2day,n2hour,n2min,n2sec] = datevec((n2time+startTime)/86400 + ...
    datenum(1970,1,1,0,0,0));
   n2abstime = (datenum(n2year,n2month,n2day,n2hour,n2min,n2sec)- ...
       datenum(startyear,startmonth,startday,0,0,0));
   n2yearday = datenum(n2year,n2month,n2day,n2hour,n2min,n2sec) - ...
       datenum(n2year,1,1,0,0,0) + 1;
   clear n2month n2day n2hour n2min n2sec;
   ncoutn2;
else
   disp (' ISUS data file not found ')
   n2time = 0;
end


%%-------------------------------------------
%% Load Ice Profiler (IPS) data if it exists
%%-------------------------------------------
if ( exist('dataIPS') == 2 )
   disp (' Loading Ice Profiler (IPS) data, matrix = ips4log ')
   dataIPS
elseif ( exist('dataIps4') == 2 )
   disp (' Loading Ice Profiler (IPS) data, matrix = ips4log ')
   dataIps4
else
   disp (' Ice Profiler (IPS) data file not found ')
end
if ( exist('ips4log') == 1 )
  % truncate Ice Profiler (IPS) data if it started recording before the navigation data did
   ips4log(:,1) = ips4log(:,1) - startTime;
   ips4log = NegTrunc(ips4log);
  % define variables of interest
   ipstime = ips4log(:,time);
   ipstravtime = ips4log(:,ipsTravelTime);
   ipsdraft = ips4log(:,ipsDraft);
   ipspres = ips4log(:,ipsPressure); % (db)
  % subtract out average atmospheric pressure from IPS pressure
   ipspres = ipspres - mean(ipspres(1:20));
   clear ips4log;
  % calculate or estimate speed of sound in water
%   if (calcsoundvel == 1)
%   else
      if (nomLAT > 70)
         soundvel = 1410;  %(m/s)
      else
         soundvel = 1450;  % (m/s)    
      end
%   end
   ipsrange = (ipstravtime * soundvel)/2;
  % calculate depth of Ice Profiler
   ipsdepth = sw_dpth(ipspres,nomLAT);
   ipsdraftraw = ipsdepth - ipsrange;
  % determine Ice Profiler timeset year and decimal day from n2time variable
   [ipsyear,ipsmonth,ipsday,ipshour,ipsmin,ipssec] = datevec((ipstime+startTime)/86400 + ...
    datenum(1970,1,1,0,0,0));
   ipsabstime = (datenum(ipsyear,ipsmonth,ipsday,ipshour,ipsmin,ipssec)- ...
       datenum(startyear,startmonth,startday,0,0,0));
   ipsyearday = datenum(ipsyear,ipsmonth,ipsday,ipshour,ipsmin,ipssec) - ...
       datenum(ipsyear,1,1,0,0,0)+1;
   clear ipsmonth ipsday ipshour ipsmin ipssec;
  % convert and store data in netcdf format
   ncoutips;

end


%%-----------------------------------------------------------
%% Define averaging time frame from primary CTD data 
%%-----------------------------------------------------------
% Averaging time is the maximum of the following
%     1
%     integer representing the maximum SBE25 timestep from the primary boardset data file
%     integer representing the maximum SBE25 timestep from the secondary boardset data file
% avgint = max( 1, ceil(max(c1timestep)), ceil(max(c2timestep)) );
avgint = 1;
if ((exist('c1time') == 1) & (exist('c2time') == 1))
   if (length(c1time) > length(c2time))
      mtime = [ceil(min(c1time)):avgint:floor(max(c1time))]';
   else
      mtime = [ceil(min(c2time)):avgint:floor(max(c2time))]';
   end      
elseif (exist('c1time') == 1)
   mtime = [ceil(min(c1time)):avgint:floor(max(c1time))]';
elseif (exist('c2time') == 1)
   mtime = [ceil(min(c2time)):avgint:floor(max(c2time))]';
else
    disp(' CAN''T FIND TIME VARIABLE FOR BINNING - exiting ');
    return
end
% determine binned timeset year and decimal day from mtime variable
   [binyear,binmonth,binday,binhour,binmin,binsec] = datevec((mtime+startTime)/86400 + ...
    datenum(1970,1,1,0,0,0));
   binabstime = (datenum(binyear,binmonth,binday,binhour,binmin,binsec)- ...
       datenum(startyear,startmonth,startday,0,0,0));
   binyearday = datenum(binyear,binmonth,binday,binhour,binmin,binsec) - ...
       datenum(binyear,1,1,0,0,0) + 1;



%%-----------------------------------
%% Bin all data
%%-----------------------------------
%% All variables will be averaged to 1 second, except lat/lon which will be interpolated to ctd time.
%% For October 2001 ALTEX AUV missions, ctd temporal resolution was ~.25 seconds.
%%  The vehicle speed relative to the water mass was ~ 1 m/s

disp([' Time bin (seconds) is ' int2str(avgint) '.  Binning data ...']);
if ( exist('c1time') == 1)
   T1=repmat(NaN,size(mtime)); C1=T1; S1=T1; D1=T1; OBSbin=T1; O2=T1; FLOW1=T1; 
   depthnav=T1; depth1=T1; pitch=T1; roll=T1; yaw=T1; Xpos=T1; Ypos=T1;
   if ( exist('c2time') == 1)
      T2=repmat(NaN,size(mtime)); C2=T2; S2=T2; D2=T2; FLOW2=T2; depth2=T2;
[T1,T2,C1,C2,S1,S2,D1,D2,OBSbin,O2,FLOW1,FLOW2,depthnav,depth1,depth2,pitch,roll,yaw,Xpos,Ypos] = ...
    bin_mission2NST(t1,t2,c1,c2,s1,s2,d1,d2,obs,o2,flw1,flw2,dpth1,dpth2, ...
    ndepth,npitch,nroll,nyaw,nposx,nposy,mtime,c1time,c2time,navtime);
   else
      [T1,C1,S1,D1,OBSbin,O2,FLOW1,depthnav,depth1,pitch,roll,yaw,Xpos,Ypos] = ...
        bin_mission1aNST(t1,c1,s1,d1,obs,o2,flw1,dpth1,ndepth,npitch,nroll,nyaw, ...
        nposx,nposy,mtime,c1time,navtime);
   end
elseif ( exist('c2time') == 1)
   T2=repmat(NaN,size(mtime)); C2=T2; S2=T2; D2=T2; FLOW2=T2; depthnav=T2; 
   depth2=T2; pitch=T2; roll=T2; yaw=T2; Xpos=T2; Ypos=T2;
    
   [T2,C2,S2,D2,FLOW2,depthnav,depth2,pitch,roll,yaw,Xpos,Ypos] = ...
     bin_mission1bNST(t2,c2,s2,d2,flw2,dpth2,ndepth,npitch,nroll,nyaw, ...
     nposx,nposy,mtime,c2time,navtime);
end

% print calibrated and binned data to netcdf file
ncoutbinned;

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

return
