%% 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 = 10;
% 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 = 10;
%
% Fixes to so that this will work with the old CTD code that
% is still on the vehicle.
%
% Serial Numbers:  Hardcoded for now.  rsm 01/10/9.
TempSerialNumber = 4122;
CondSerialNumber = 2514;
%
humidity = 8;
flow1 = 9;
flow2 = 9;

%%------------------------
%% 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 = 82.5;  % (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;
   if( exist('ctdDriverctdlog')==1 )
     %
     % Fixes to so that this will work with the old CTD code that
     % is still on the vehicle.
     %
     % Serial Numbers:  Hardcoded for now.  rsm 01/10/9.
     TempSerialNumber = 4122;
     CondSerialNumber = 2514;
     %
     humidity = 8;
     flow1 = 9;
     ctdDriver = ctdDriverctdlog;
   end   
   % 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';

%
% ptime and pres are from parosci.log.  
%
    pres1  = interp1( ptime,   pres,   c1time );
    pitch1 = interp1( navtime, npitch, c1time );

% 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);
   if( 0 )
   figure; clf; 
     plot(c1timestep) 
     xlabel('datapoint')
     ylabel('Primary Boardset Timestep (s)')
     title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
     grid
   end
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;
   if( exist('seabirdctdlog')==1 )
     % Serial Numbers: Hardcoded due to the old code on the vehicle. rsm
     % 01/10/9
     TempSerialNumber = 4123;
     CondSerialNumber = 2478;
     flow2 = 9;
     ctdDriver2 = seabirdctdlog;
   end
   % 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';

% calculate/interpolate pressure at CT sensors for each data point and
% account for pitch of vehicle

   pres2  = interp1( ptime,   pres,   c2time );
   pitch2 = interp1( navtime, npitch, c2time );
   
% 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);
   if( 0 )
   figure; clf; 
     plot(c2timestep) 
     xlabel('datapoint')
     ylabel('Secondary Boardset Timestep (s)')
     title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
     grid
   end
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

   pitchn2 = interp1( navtime, npitch, n2time);
   n2depth = interp1( navtime, ndepth, n2time);
   
   % 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;
    
   if( 0 )
   figure;
   plot(itimestep,'b')
   xlabel('datapoint')
   ylabel('ISUS Timestep (s)')
   title(['Time between datapoints, Data file: ',missionID,'  Plotted on ',dateS])
   grid
   end
% 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 & length(ips4log) > 20)
  % 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)
   ipspresraw = ipspres;
  % 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 = 1455;  %(m/s)
      else
         soundvel = 1450;  % (m/s)    
      end
%   end
   ipsrange = (ipstravtime * soundvel)/2;
   ipspitch = interp1( navtime, dataNav(:,mTheta),  ipstime );
   ipsphi   = interp1( navtime, dataNav(:,mPhi),    ipstime );
  % calculate depth of Ice Profiler
   ipsdepth = sw_dpth(ipspres,nomLAT);
   ipsdraftraw = ipsdepth - ipsrange.*cos(ipspitch).*cos(ipsphi);
  % 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)); 
    P1=T1; 
    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;
    N2     = T1;
    N2Qual = T1;
    N2Depth= T1;
  end

  if ( exist('c2time') == 1)
    T2=repmat(NaN,size(mtime)); 
    P2=T2; 
    C2=T2; 
    S2=T2; 
    D2=T2; 
    FLOW2=T2; 
    depth2=T2;
  end

  buf=median(diff(mtime))/2;

  for ii=1:length(mtime)
    %%------
    %% CTD1
    %%------
    if( exist('c1time') == 1 )
      idx = find(c1time >= mtime(ii)-buf & c1time < mtime(ii)+buf);
      if (isempty(idx) == 0)
	P1(ii,1)=median(p1(idx));
	T1(ii,1)=median(t1(idx));
	C1(ii,1)=median(c1(idx));
	S1(ii,1)=median(s1(idx));
	D1(ii,1)=median(d1(idx));
	OBSbin(ii,1)=median(obs(idx));
	O2(ii,1)=median(o2(idx));
	FLOW1(ii,1)=median(flw1(idx));
	depth1(ii,1)=median(dpth1(idx));
      end
    end %if( exist('c1time') == 1 )
    %%------
    %% CTD2
    %%------
    if( exist('c2time') == 1 )
      idx = find(c2time >= mtime(ii)-buf & c2time < mtime(ii)+buf);
      if (isempty(idx) == 0)
	P2(ii,1)=median(p2(idx));
	T2(ii,1)=median(t2(idx));
	C2(ii,1)=median(c2(idx));
	S2(ii,1)=median(s2(idx));
	D2(ii,1)=median(d2(idx));
	FLOW2(ii,1)=median(flw2(idx));
	depth2(ii,1)=median(dpth2(idx));
      end
    end %if( exist('c2time') == 1 )

    %%----------
    %% Spatial 
    %%----------
    idx = find(navtime >= mtime(ii)-buf & navtime < mtime(ii)+buf);
    if (isempty(idx) == 0)
      depthnav(ii,1)=median(ndepth(idx));
      pitch(ii,1)=median(npitch(idx));
      roll(ii,1)=median(nroll(idx));
      yaw(ii,1)=median(nyaw(idx));
      Xpos(ii,1)=median(nposx(idx));
      Ypos(ii,1)=median(nposy(idx));
    else
      depthnav(ii,1)=0;
      pitch(ii,1)=0;
      roll(ii,1)=0;
      yaw(ii,1)=0;
      % Xpos(ii,1)=Xpos(1,1);
      % Ypos(ii,1)=Ypos(1,1);
      Xpos(ii,1)=0;
      Ypos(ii,1)=0;
    end

    %%------
    %% ISUS
    %%------
    if( exist('n2') == 1 )
      idx = find(n2time >= mtime(ii)-buf & n2time < mtime(ii)+buf);
      if (isempty(idx) == 0)
	N2(ii,1)      = median(n2(idx));
	N2Qual(ii,1)  = median(n2qual(idx));
	N2Depth(ii,1) = median(n2depth(idx));
      end
    end %if( exist('n2') == 1 )


    %%------
    %% Ips
    %%------
    if( exist('ips4log') == 1 )
      idx = find(ipstime >= mtime(ii)-buf & ipstime < mtime(ii)+buf);
      if (isempty(idx) == 0)
	IPSPRES(ii,1)     = median(ipspres(idx));
	IPSDEPTH(ii,1)    = median(ipsdepth(idx));
	IPSDRAFT(ii,1)    = median(ipsdraft(idx));
	IPSDRAFTRAW(ii,1) = median(ipsdraftraw(idx));
      end
    end %if( exist('ips4log') == 1 )

  end %for  
% print calibrated and binned data to netcdf file
  ncoutbinned;
% Plot
% set(gcf,'color',[1 1 1],'units','pixels','position',[10 40 1200 800]);
% plotallNST
  plotsci

  return
