%% process2: process individual AUV mission
%%
%%
%% USAGE:
%% From within the AUV mission directory, issue process2 command.
%%
%% In the current data acquisition system, there are separate time stamps for [ctd], [lat lon], 
%% [roll pitch yaw], [pressure depth], [ice profiler].
%% 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 (0.25 seconds nominally)
%% dataIsus (nitrogen) ~ 0.4 Hz (2.5 seconds)
%% dataIps4 (ice profiler) ~ 0.5 Hz for pressure, 1 Hz for other data (2 seconds, 1 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.
%%
%% 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).
%%
%%
%% process2 calls other m-files: 
%% get_spatialNST.m, DirName.m
%% sw_dpth.m, sw_dens.m, sw_salt.m, sw_c3515.m, sw_cndr.m
%% Pres_cal.m, Temp_cal.m, Cond_cal.m, O2_cal_SBE43.m
%% calc_paros_pres.m, calc_temp.m, calc_cond.m, calc_O2_SBE43.m
%% plotsci2.m, NegTrunc.m
%% ncoutn2.m, ncout1.m, ncout2.m, ncoutbinned.m, ncoutips.m
%% 
%% ALTEX October 2001 cruise notes:
%% The drivers for the various science instruments were not well tested and finalized prior
%% to the cruise, so some of them changed over the course of the experiment.  In particular, 
%% the CTD drivers were not appropriately updated for at least the first half of the cruise.
%% For this reason, the column designators for flow1, flow2, humidity, dissolved oxygen, and
%% OBS had to be "hard-coded" into this script because of errors in the log files from the
%% experiment.
%%
%% Future use of this code will require modifications of the tubing volume calculations for the
%% pumped sensors as well as possible changes to the CT processing sections if a new sensor 
%% package is used for CT data collection.
%% 
%%
%% John Ryan 6/8/01 
%% Nicole Tervalon  Sep 2001, Aug 2002
%% 


%clear all; close all;

fprintf('\n Processing mission %13s on date %11s\n\n', DirName,date);

%%---------------------------------------------
%% Vehicle configuration and other user inputs
%%---------------------------------------------
% Horizontal length of vehicle midbody between tail section joining ring and science payload
% joining ring
midlength = 67.375;  % inches including battery section only
% midlength = 67.375 + 56.0;  % inches including battery and buoy section
% Horizontal distance from center of MVC sphere to outside of tail section joining ring
taillength = 15.5;  % inches
% Horizontal distance from center of conductivity and temperature sensors to outside of nose
% section joining ring
noselength = 36;  % inches
% Vertical distance between vehicle midline and main vehicle pressure sensor
pres1length = 8.5;  % inches
% Vertical distance bewteen vehicle midline and ice profiler pressure sensor
pres2length = 8.0;  % inches
% Horizontal distance from nose section joining ring to ice profiler pressure sensor
icelength = 4.5;  % inches
% Sampling interval for CT sensors
sampint = 0.25;  % seconds
% Alpha and beta constants for conductivity thermal mass calculation
alphac = 0.04;
betac = 1/8.0;
% 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 = 1;
% Averaging interval for binning data in seconds after calibration and temporal lag incorporation
%  entering a 0 will result in no averaging
avgint = input('---> Enter desired averaging time in seconds for the final processed data (default = 1): ');
if isempty(avgint)
%   avgint = 1.25;
   avgint = 1;
end
% Round to nearest integer if not an integer
avgint = round(avgint);
fprintf('---> The averaging interval that will be used to filter the data %3d seconds\n\n',avgint);
    
% 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;

%--------------------------------------------------------------
% Volume calculations for temporal lag between sensors (in mL)
%--------------------------------------------------------------
in2cm = 2.54*2.54*2.54;  % conversion from cubic inches to cubic centimeters
tcduct = 2.2;  % tc duct volume to first sensing element in cond cell
cvol = 2;  % active volume of conductivity cell (mL)
o2vol = 12.5;  % volume of dissolved oxygen sensor flow through (mL)
n2vol = 3;  % volume of nitrogen (ISUS) sensor flow through (mL)
tmp2c = tcduct + cvol/2;  % volume btw temp probe and center of cond cell
c2o2 = cvol/2 + (1.8+2.6) + 14 + o2vol/2; % vol btw center of cond and o2 cell
o22n2 = o2vol/2 + 27 + n2vol/2;  % vol btw o2 cell and ISUS sensor


%%------------------------
%% Load Spatial Data 
%% Note: yaw is heading
%%------------------------
% [glat,glon,gpstime,ndepth,nroll,npitch,nyaw,navtime,nposx,nposy] = get_spatial(1);
get_spatialNST
% move roll data to -180 to 180 range instead of 0 to 360 degrees
   indx=find(nroll >= 180);
   nroll(indx) = nroll(indx) - 360;
   clear indx
% calculate approximate derivative of pitch variations
   fprintf(' Mean of vehicle pitch for this mission is %5.1f degrees\n\n',mean(npitch))
   difpitch = diff(npitch);
   diftime = diff(navtime);
   derivpitch = difpitch./diftime;
   derivpitchtime = navtime(1:(length(navtime)-1)) + diff(navtime)/2;
%% 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 = input('---> Please enter the latitude of the mission in decimal degrees North: ');
latrad = nomLAT*pi/180;
grav = 9.780318 * (1.0 + sin(latrad)*sin(latrad)*(5.2788E-3 + 2.36E-5 * sin(latrad)*sin(latrad)));
fprintf('---> At latitude %7.3f N, the gravitational constant is %6.4f m/s^2\n\n',nomLAT,grav); 
% nomLAT = startlat;

%%-----------------------
%% Load Pressure Data
%%-----------------------
if ( exist('dataPS') == 2 )
   fprintf(' Loading Pressure data\n')
   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);
   PresSerialNumber = 77988; %for ALTEX and previous data
%   PresSerialNumber = 87760; %for SPOKES and after data
   % pres_serialnum = PresSerialNumber;
   if calpress == 1
   % calibrate pressure sensor
     [pres,prestc] = calc_paros_pres(parosci(:,pres_period),parosci(:,temp_period),PresSerialNumber); % (in bars)
   else
      pres = parosci(:,pressure);
   end
   parostemp = parosci(:,temp);
   parosdpth = parosci(:,depth);
   if (length(pres) >= 50)
      atmpres = mean(pres(1:50));  % (bar)
   else
      atmpres = mean(pres(1:length(pres)));  % (bars)
   end
% Set pressure measurement to RELATIVE pressure, NOT ABSOLUTE (this assumes that the pressure sensor
%   is at atmospheric pressure when the vehicle is floating on the surface)
   fprintf('   MVC pressure sensor mean atmospheric pressure is %8.4f bars\n',atmpres);
% Request input from user of atmospheric pressure collected by instrument when it was deck tested
% prior to deployment
   atmospres = input('---> Enter the est atm pressure from AUV deck tests in bars [above value is default] : ');
   if (isempty(atmospres))
       atmospres = atmpres;
   end
   fprintf('---> Atmospheric pressure being used is : %8.4f bars\n\n',atmospres);
% Plot fft of pressure data to view aliasing
%   plotfft(pres,1024)
   absauvpres = pres;
% boxcar averge signal for approximately 2 seconds to smooth data (0 phase shift)
   if (length(pres) > 3*10)
      presfilt = filtfilt((ones(1,10)/length(ones(1,10))), 1, pres);
   else
      fprintf(' ** Not filtering pressure because file length too short\n');
      presfilt = pres;
   end
   pres = presfilt - atmospres;
   parosdpthcalc = sw_dpth(pres*10,nomLAT);
   navdepthp=interp1(navtime,ndepth,ptime);
   clear parosci;
else
   fprintf(' !!! No pressure data file found - exiting\n')
   return
end



%%----------------------
%% Load CTD data
%%----------------------
%% Look for PRIMARY CTD data
if ( exist('dataCTD') == 2 )
   fprintf(' Loading primary CTD data\n')
   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.1;
     CondSerialNumber = 2514.1;
   % use the following if missions run after 24 July 2002
 %    TempSerialNumber = 4122;
 %    CondSerialNumber = 2514;
  % Row numbers hardcoded because using old ctd drivers
     humidity = 8;
     flow1 = 9;
     OBS = 10;
     dissolvedO2 = 11;
     ctdDriver = ctdDriverctdlog;
   end   
   % FILTER DATA FOR BAD TIME STAMPING (for data collected prior to 2001, yearday 235)
      fprintf('     Filtering old data for bad timestamps if necessary\n')
 % The following lines commented out by NST 25 June 2002
%      diffcondtemp = (diff(ctdDriver(:,conductivity)) | diff(ctdDriver(:,temperature)));
%      ctdlogtemp = ctdDriver(find(diffcondtemp == 1),:);
%      ctdlog = ctdlogtemp(find(diff(ctdlogtemp(:,time)) ~= 0),:);
%      clear ctdDriver diffcondtemp ctdlogtemp;
   % OR use the following definition
   ctdlog = ctdDriver;
   diffcondtemp = (diff(ctdlog(:,cond_frequency)) | diff(ctdlog(:,temp_frequency)));
   ctdlogtemp1 = find(diffcondtemp ~= 1);
   clear diffcondtemp
%   timestamps = find(diff(ctdlog(:,time) ~= 0));
%   pause
   
   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);
% recalculate SBE25 timestamping based on actual data output, not MVC storage timestamp
   newc1time=[1:length(c1time)]'*0.25-0.25+c1time(1);
% determine SBE25 boardset measured year and decimal day from time variable
   [c1year,c1month,c1day,c1hour,c1min,c1sec] = datevec((newc1time+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;
% calibrate flow data using current estimates (to L/min)
   rawflw1 = (exp((ctdlog(:,flow1)+4.7685)/1.4816)-3.92)/111.5;
   rawflw1(1) = 0;
   rawflw1(find(ctdlog(:,flow1) <= 0)) = 0;
%  rawflw1(find(isnan(ctdlog(:,flow1)))) = 0;
% boxcar averge flow for approximately 1 second to smooth data (0 phase shift)
   flw1 = filtfilt([0.2 0.2 0.2 0.2 0.2], 1, rawflw1);
% calculate time lags caused by flow rate and temporally offset data
   % don't use first 700 data points (approximately 150 seconds) because flow through pumped
   % system isn't stable yet and stop data points approximately 8 seconds before end of data
   % set to allow for time lag shifting
   scistart = max(min(find(flw1>=mean(flw1))),700) + 20;
   newt1time = newc1time(scistart:(length(newc1time)-(8/0.25)));
   % calculate measured year and decimal day from shifted and lagged time variable
   [c1Ayear,c1Amonth,c1Aday,c1Ahour,c1Amin,c1Asec]=datevec((newt1time+startTime)/86400 + ...
      datenum(1970,1,1,0,0,0));
   c1Aabstime = (datenum(c1Ayear,c1Amonth,c1Aday,c1Ahour,c1Amin,c1Asec)- ...
      datenum(startyear,startmonth,startday,0,0,0));
% Convert flow to mL/s and clip beginning of mission data to remove startup lag
   filtflw1 = flw1(scistart:(length(newc1time)-(8/0.25)))*1000/60;
% get sensor data
   t1freq = ctdlog(scistart:(length(newc1time)-(8/0.25)),temp_frequency);  % (Hz)
   c1freqraw = ctdlog(:,cond_frequency)./1000;  % (kHz)
% have to adjust obs start and stop times to align with other science sensors in order to
% generate consistent netCDF output file with depth, pressure, salinity, temperature, oxygen, 
% and obs data in it
   obsraw = ctdlog(:,OBS);
   obs = ctdlog(scistart:(length(newc1time)-(8/0.25)),OBS);  % (volts)
   o2vraw = ctdlog(:,dissolvedO2);
% calibrate temperature (deg C)
   fprintf('   Calibrating Temperature\n');
   t1 = calc_temp(t1freq,temp1_serialnum);
% determine time lags between sensors
   shiftA1a = tmp2c./filtflw1;
   shiftB1a = c2o2./filtflw1;
   shiftC1a = o22n2./filtflw1;

% integrate flow rate to get the time shift between each sensor given the volumes between them
% Temperature to Conductivity
   flowtmptime = [newt1time(1):0.01:(newt1time(length(newt1time))+30)];
   flowtmp = interp1(newt1time,filtflw1,flowtmptime);
   flowtmp(((length(newt1time)-1)*25+1):length(flowtmp))=flowtmp((length(newt1time)-1)*25+1);
   flowtmp(find(flowtmp < 5)) = 5;
   if (find(flowtmp < 5))
      fprintf(' WARNING!!!!!  Flow rate is low during some portion of this data set\n')
   end       
   cumvol = (cumsum(flowtmp)-flowtmp(1))*0.01;
   fprintf ('        Shifting raw conductivity timestamp \n')
   for i = 1:length(newt1time)
       shifttmp2c(i) = interp1((cumvol(((i-1)*25+1):((i-1)*25+1+100)) - ...
           cumvol((i-1)*25+1)), ...
           flowtmptime(((i-1)*25+1):((i-1)*25+1+100)),tmp2c);
   end
   shifttmp2c = shifttmp2c';

% Temperature to Oxygen
   fprintf ('        Shifting raw oxygen timestamp (used only if oxygen is being collected) \n')
   for i = 1:length(newt1time)
       shiftc2o2(i) = interp1((cumvol(((i-1)*25+1):((i-1)*25+1+500)) - ...
           cumvol((i-1)*25+1)), ...
           flowtmptime(((i-1)*25+1):((i-1)*25+1+500)),(tmp2c+c2o2));
   end
   shiftc2o2 = shiftc2o2';

% Temperature to Nitrate
   fprintf ('        Shifting raw nitrate timestamp (used only if nitrate is being collected)\n')
   for i = 1:length(newt1time)
       shifto22n2(i) = interp1((cumvol(((i-1)*25+1):((i-1)*25+1+1000)) - ...
           cumvol((i-1)*25+1)), ...
           flowtmptime(((i-1)*25+1):((i-1)*25+1+1000)),(tmp2c+c2o2+o22n2));
   end
   shifto22n2 = shifto22n2';
   clear flowtmptime flowtmp cumvol

       
% shift conductivity, oxygen, and nitrate uncalibrated data in time such that at the
% time that the temperature measurement is made, the lag accommodated data for the 
% other sensors is stored at that same timestamp
   c1freqtime = shifttmp2c-0.035;
   shiftA1 = shifttmp2c - newt1time;
   o2vtime = shiftc2o2;
   shiftB1 = shiftc2o2 - newt1time - shiftA1;
   n2vtime = shifto22n2;
   shiftC1 = shifto22n2 - newt1time - shiftA1 - shiftB1;
   clear shifttmp2c shiftc2o2 shifto22n2

   c1freq = interp1(newc1time, c1freqraw, c1freqtime);
   o2v = interp1(newc1time, o2vraw, o2vtime);
   
% calculate/interpolate pressure at CT sensors for each data point and account for pitch of vehicle
% (pressure in bars, pitch in degrees)
% ptime and pres are from parosci.log.  
    pres1  = interp1( ptime,   pres,   newt1time );
    pitch1 = interp1( navtime, npitch, newt1time );
% calibrate conductivity (Siemens/meter) [using uncorrected pressure because conductivity
% dependence on pressure is very small (on the order of 10^-7)]
   fprintf('   Calibrating Conductivity\n');
   c1uncorrected = calc_cond(c1freq,t1,pres1,cond1_serialnum);
% correct conductivity for cell thermal mass (see Seabird documentation for explanation of equations)
   ctm1(1) = 0;
   for i = 1:(length(c1uncorrected)-1)
      ctm1(i+1) = (-1.0*(1-(2*(2*alphac/(sampint*betac+2))/alphac))*ctm1(i)) + ...
          (2*(alphac/(sampint*betac+2))*(0.1*(1+0.006*(t1(i)-20)))*(t1(i+1)-t1(i)));
   end
   c1 = c1uncorrected + ctm1';
   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   %   assumption is that the CT sensor intake is aligned with the center of the MVC sphere in 
   %   AUV body space
   poffset = abs(pres1length * cos(pitch1*pi/180));
   % Vertical distance from center of MVC sphere to CT sensors (inches) 
   MVCoffset = abs((taillength + midlength + noselength) * sin(pitch1*pi/180));
   % Total vertical distance between pressure sensor and CT sensors (meters)
   % Note that pitch is negative when nose is down and we want the pressure to be MORE when the
   % vehicle is pitch down (note that pitch1 is negative when the vehicle is nose down)
   totdeltad1 = - sign(pitch1).*(MVCoffset - sign(pitch1).*poffset)/39.37; 
   % Calculate approximate saltwater density (in order to calculate delta pressure)
   %   sw_c3515 in 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 * grav .* 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 nominal error in salinity (in psu) if there is an error in conductivity of +0.003
% using the maximum and minimum salinity values collected during the mission
   aa = find(s1 == max(s1));
   cratioaa=(c1(aa)+0.002)*10/sw_c3515;
   saa = sw_salt(cratioaa,t1(aa),pres1(aa)*10);
   densaa = sw_dens(saa,t1(aa),pres1(aa)*10);
   p1aa = pres1(aa) + (densaa * 9.81 .* totdeltad1(aa))/100000;
   s1aa = sw_salt(cratioaa,t1(aa),p1aa*10);
% calculate errors in salinity based on errors in conductivity values of order 0.0003 (from specs) 
   bb = find(s1 == min(s1(600:length(s1))));
   cratiobb=(c1(bb)+0.002)*10/sw_c3515;
   sbb = sw_salt(cratiobb,t1(bb),pres1(bb)*10);
   densbb = sw_dens(sbb,t1(bb),pres1(bb)*10);
   p1bb = pres1(bb) + (densbb * 9.81 .* totdeltad1(bb))/100000;
   s1bb = sw_salt(cratiobb,t1(bb),p1bb*10);
   fprintf('       Change in salinity if max measured conductivity is modified by +0.002 S/m : %8.4f psu\n',(s1(aa)-s1aa));
   fprintf('                          if min measured conductivity is modified by +0.002 S/m : %8.4f psu\n',s1(bb)-s1bb);
% calculate depth (meters) and set gauge pressure to 0 if it is negative
   p1(find(p1<0)) = 0;
   dpth1 = sw_dpth(p1*10,nomLAT);
% calibrate oxygen
   fprintf('   Computing Oxygen\n');
   % 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,newt1time,ounits);
% clear data array  
   clear ctdlog;
% print calibrated data to netcdf file
   ncout1;
% check for bad ctd dacq
   if( 0 )
     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
   end
else
   fprintf(' Primary CTD log file not found\n');
end

fprintf('\n');
%% Look for SECONDARY CTD data

if ( exist('dataCTD2') == 2 )
   fprintf (' Loading secondary CTD data\n')
   dataCTD2;
   if( exist('seabirdctdlog')==1 )
     % Serial Numbers: Hardcoded due to the old code on the vehicle. rsm
     % 01/10/9
     TempSerialNumber = 4123.1;
     CondSerialNumber = 2478.1;
   % use the following if missions run after 24 July 2002
%     TempSerialNumber = 4123;
%     CondSerialNumber = 2478;
     flow2 = 9;
     ctdDriver2 = seabirdctdlog;
   else
   % use the following if missions run after 24 July 2002
     TempSerialNumber = 4123;
     CondSerialNumber = 2478;
     flow2 = 9;
    ctdDriver2 = ctdDriverctdlog;
   end
   % FILTER DATA FOR BAD TIME STAMPING (for data collected prior to 2001, yearday 235)
      fprintf ('     Filtering old data for bad timestamps if necessary\n')
 % The following lines commented out by NST 25 June 2002
%      diffcondtemp = (diff(ctdDriver2(:,conductivity)) | diff(ctdDriver2(:,temperature)));
%      ctdlogtemp = ctdDriver2(find(diffcondtemp == 1),:);
%      ctdlog = ctdlogtemp(find(diff(ctdlogtemp(:,time)) ~= 0),:);
%      clear ctdDriver2 diffcondtemp ctdlogtemp;
   % OR use the following definition
   ctdlog = ctdDriver2;
   diffcondtemp = (diff(ctdlog(:,cond_frequency)) | diff(ctdlog(:,temp_frequency)));
   ctdlogtemp2 = find(diffcondtemp ~= 1);
   clear diffcondtemp
%   timestamps = find(diff(ctdlog(:,time) ~= 0)):
%   pause
   
   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);
% recalculate SBE25 timestamping based on actual data output, not MVC storage timestamp
   newc2time=[1:length(c2time)]'*0.25-0.25+c2time(1);
% determine SBE25 boardset measured year and decimal day from time variable
   [c2year,c2month,c2day,c2hour,c2min,c2sec] = datevec((newc2time+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;
% calibrate flow data using current estimates (to L/min)
   rawflw2 = (exp((ctdlog(:,flow2)+4.7685)/1.4816)-3.92)/111.5;
   rawflw2(1) = 0;
   rawflw2(find(ctdlog(:,flow2) <= 0)) = 0;
% boxcar average flow for approximately 1 second to smooth data (0 phase shift)
   flw2 = filtfilt([0.2 0.2 0.2 0.2 0.2], 1, rawflw2);
% calculate time lags caused by flow rate and temporally offset data
   % don't use first 600 data points (approximately 150 seconds) because flow through pumped
   % system isn't stable yet and stop data points approximately 8 seconds before end of data
   % set to allow for time lag shifting
   % Only use flow data after it reaches 80% of the mean or after the initial startup (160 data points)
%   scistart = max(min(find(diff(flw1)>=0.3)),160) + 20;
   scistart2 = max(min(find(flw2>=mean(flw2))),700) + 20;
   newt2time = newc2time(scistart2:(length(newc2time)-(8/0.25)));
   % calculate measured year and decimal day from shifted and lagged time variable
   [c2Ayear,c2Amonth,c2Aday,c2Ahour,c2Amin,c2Asec]=datevec((newt2time+startTime)/86400 + ...
      datenum(1970,1,1,0,0,0));
   c2Aabstime = (datenum(c2Ayear,c2Amonth,c2Aday,c2Ahour,c2Amin,c2Asec)- ...
      datenum(startyear,startmonth,startday,0,0,0));
% Convert flow to mL/s and clip beginning of mission data to remove startup lag
   filtflw2 = flw2(scistart2:(length(newc2time)-(8/0.25)))*1000/60;
% rename variables
   t2freq = ctdlog(scistart2:(length(newc2time)-(8/0.25)),temp_frequency);  % (Hz)
   c2freqraw = ctdlog(:,cond_frequency)./1000;  % (kHz)
% calibrate temperature (deg C)
   fprintf('   Calibrating Temperature\n');
   t2 = calc_temp(t2freq,temp2_serialnum);
% determine time lags between sensors in a general sense
   shiftA2a = tmp2c./filtflw2;

% integrate flow rate to get the time shift between each sensor given the volumes between them
% Temperature to Conductivity
   fprintf ('        Shifting raw conductivity data \n')
   flowtmptime = [newt2time(1):0.01:(newt2time(length(newt2time))+30)];
   flowtmp = interp1(newt2time,filtflw2,flowtmptime);
   flowtmp(find(flowtmp < 5)) = 5;
   if (find(flowtmp < 5))
      fprintf(' WARNING!!!!!  Flow rate is 0 or very low during some portion of this data set\n')
   end       
   flowtmp(((length(newt2time)-1)*25+1):length(flowtmp))=flowtmp((length(newt2time)-1)*25+1);
   cumvol = (cumsum(flowtmp)-flowtmp(1))*0.01;
   for i = 1:length(newt2time)
       shifttmp2c(i) = interp1((cumvol(((i-1)*25+1):((i-1)*25+1+100)) - ...
           cumvol((i-1)*25+1)), ...
           flowtmptime(((i-1)*25+1):((i-1)*25+1+100)),tmp2c);
   end
   clear flowtmptime flowtmp cumvol
   shifttmp2c = shifttmp2c';
 
   
% shift conductivity uncalibrated data in time such that at the
% time that the temperature measurement is made, the lag accommodated data for the 
% other sensors is stored at that same timestamp
   c2freqtime = shifttmp2c-0.035;
   shiftA2 = shifttmp2c - newt2time;
   clear shifttmp2c
%   c2freqtime = newt2time + shiftA2;
   c2freq = interp1(newc2time, c2freqraw, c2freqtime);
   
% calculate/interpolate pressure at CT sensors for each data point and account for pitch of vehicle
% (pressure in bars, pitch in degrees)
% ptime and pres are from parosci.log.  
   pres2  = interp1( ptime,   pres,   newt2time );
   pitch2 = interp1( navtime, npitch, newt2time );
   
% calibrate conductivity (Siemens/meter) [using uncorrected pressure
% because conductivity dependence on pressure is very small (on the order
% of 10^-7)] 
   fprintf('   Calibrating Conductivity\n');
   c2uncorrected = calc_cond(c2freq,t2,pres2,cond2_serialnum);
% correct conductivity for cell thermal mass (see Seabird documentation for explanation of equations)
   ctm2(1) = 0;
   for i = 1:(length(c2uncorrected)-1)
      ctm2(i+1) = (-1.0*(1-(2*(2*alphac/(sampint*betac+2))/alphac))*ctm2(i)) + ...
          (2*(alphac/(sampint*betac+2))*(0.1*(1+0.006*(t2(i)-20)))*(t2(i+1)-t2(i)));
   end
   c2 = c2uncorrected + ctm2';

   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   %   assumption is that the CT sensor intake is aligned with the center of the MVC sphere in 
   %   AUV body space
   poffset = abs(pres1length * cos(pitch2*pi/180));
   % Vertical distance from center of MVC sphere to CT sensors (inches) 
   MVCoffset = abs((taillength + midlength + noselength) * sin(pitch2*pi/180));
   % Total vertical distance between pressure sensor and CT sensors (meters)
   % Note that pitch is negative when nose is down and we want the pressure to be MORE when the
   % vehicle is pitch down (note that pitch1 is negative when the vehicle is nose down)
   totdeltad2 = - sign(pitch2).*(MVCoffset - sign(pitch2).*poffset)/39.37; 
   % Calculate approximate saltwater density (in order to calculate delta pressure)
   %    sw_c3515 in 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 * grav .* 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);
% calculate nominal error in salinity (in psu) if there is an error in conductivity of +0.0003
% using the maximum and minimum salinity values collected during the mission
   aa = find(s2 == max(s2));
   cratioaa=(c2(aa)+0.002)*10/sw_c3515;
   saa = sw_salt(cratioaa,t2(aa),pres2(aa)*10);
   densaa = sw_dens(saa,t2(aa),pres2(aa)*10);
   p2aa = pres2(aa) + (densaa * 9.81 .* totdeltad2(aa))/100000;
   s2aa = sw_salt(cratioaa,t2(aa),p2aa*10);
   % don't use first 600 data points (approximately 150 seconds) because flow through pumped
   % system isn't stable yet
   bb = find(s2 == min(s2(600:length(s2))));
   cratiobb=(c2(bb)+0.002)*10/sw_c3515;
   sbb = sw_salt(cratiobb,t2(bb),pres2(bb)*10);
   densbb = sw_dens(sbb,t2(bb),pres2(bb)*10);
   p2bb = pres2(bb) + (densbb * 9.81 .* totdeltad2(bb))/100000;
   s2bb = sw_salt(cratiobb,t2(bb),p2bb*10);
   fprintf('       Change in salinity if max measured conductivity is modified by +0.002 S/m:  %8.4f psu\n',(s2(aa)-s2aa));
   fprintf('                          if min measured conductivity is modified by +0.002 S/m:  %8.4f psu\n',s2(bb)-s2bb);
% clear variables
   clear aa bb cratioaa cratiobb saa sbb densaa densbb p1aa p2aa p1bb p2bb;

% clear data array  
   clear ctdlog;
% print calibrated data to netcdf file
   ncout2;
% check for bad ctd dacq
   if( 0 )
     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
   end
else
   fprintf(' Secondary CTD log file not found\n')
end

if ((exist('dataCTD') == 0) & (exist('dataCTD2') == 0))
    fprintf(' !!! No CTD data found - exiting\n')
    return
end


fprintf('\n');
%%-----------------------------------
%% Load ISUS data if it exists
%%-----------------------------------
if ( exist('dataIsus') == 2 )
   fprintf(' Loading ISUS data\n');
   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);

   % Some of the data files have corrupted timestamps at the beginning of the file
   % for about 13 or 14 timesteps.  To account for this bad data, I have started
   % the data sets at a point at least 5 seconds after the beginning of the file.
   fprintf('   The first good ISUS data point used is %4d\n',max(find((isuslog(:,time)-isuslog(1,time))<5)));
   n2 = isuslog(max(find((isuslog(:,time)-isuslog(1,time))<5)):length(isuslog(:,1)),isusNitrate);
   n2qual = isuslog(max(find((isuslog(:,time)-isuslog(1,time))<5)):length(isuslog(:,1)),isusQuality);
   n2time = isuslog(max(find((isuslog(:,time)-isuslog(1,time))<5)):length(isuslog(:,1)),time);   
   n2new = interp1(n2time,n2,n2vtime);
   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)
   %   assumption is that the CT sensor intake is aligned with the center of the MVC sphere in 
   %   AUV body space
   poffset = abs(pres1length * cos(pitchn2*pi/180));
   % Vertical distance from center of MVC sphere to CT sensor intake (inches) 
   MVCoffset = abs((taillength + midlength + noselength) * sin(pitchn2*pi/180));
   % Total vertical distance between pressure sensor and CT sensor intake (meters)
   % Note that pitch is negative when nose is down and we want the pressure to be MORE when the
   % vehicle is pitch down (note that pitchn2 is negative when the vehicle is nose down)
   totdeltadn2 = - sign(pitchn2).*(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
   fprintf(' ISUS data file not found\n')
   n2time = 0;
end
fprintf('\n');


%%-------------------------------------------
%% Load Ice Profiler (IPS) data if it exists
%%-------------------------------------------
if ( exist('dataIPS') == 2 )
   fprintf(' Loading Ice Profiler (IPS) data\n')
   dataIPS
elseif ( exist('dataIps4') == 2 )
   fprintf(' Loading Ice Profiler (IPS) data\n')
   dataIps4
else
   fprintf(' Ice Profiler (IPS) data file not found\n')
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 IPS pressure sensor offset (due to pressure sensor offset error)
   ipspres = ipspres - 0.09;
  % subtract out average atmospheric pressure from IPS pressure
   fprintf('   Atmospheric pressure from Ice Profiler (bars):  %8.3f bars \n',mean(ipspres(1:20))/10);
%   disp(mean(ipspres(1:20))/10);
   ipsatmospres = input('---> Enter the est atm pressure from IPS deck tests (bars) [above value is default] : ');
   if (isempty(ipsatmospres))
       ipsatmospres = mean(ipspres(1:20))/10;
   end
   fprintf('---> The approximate IPS atmospheric pressure is %7.4f bars\n\n',ipsatmospres);
   ipspres = ipspres - ipsatmospres*10;
  %clear ips4log;
  % calculate or estimate speed of sound in water
   if (calcsoundvel == 1)
       svel1 = sw_svel(s1,t1,p1*10);
       svel1av = mean(svel1);
       fprintf('     Average sound velocity from CT 1 : %7.1f \n', svel1av);
       svel1dif = max(svel1) - min(svel1);
       svel2 = sw_svel(s2,t2,p2*10);
       svel2av = mean(svel2);
       fprintf('     Average sound velocity from CT 2 : %7.1f \n', svel2av);
       svel2dif = max(svel2) - min(svel2);
       soundvel = svel1av;
   else
      if (nomLAT > 70)
         soundvel = 1455;  %(m/s)
      else
         soundvel = 1450;  % (m/s)    
      end
   end
   ipsrange = (ipstravtime * soundvel)/2;
   % interpolate navigation pitch and roll to ips data timestamps
   ipspitch = interp1( navtime, npitch, ipstime );
   ipsphi   = interp1( navtime, nroll, ipstime );
  % calculate depth of Ice Profiler
   ipsdepth = sw_dpth(ipspres,nomLAT);
  % calculate depth of vehicle at MVC pressure sensor location from Ice Profiler pressure
  % sensor
   % Vertical distance from center of MVC sphere to pressure sensor (inches)
   poffset = abs(pres1length * cos(ipspitch*pi/180));
   % Vertical distance from center of MVC sphere to center of IPS (inches) 
   MVCoffset = abs((taillength + midlength + icelength) * sin(ipspitch*pi/180));
   % Ice profiler pressure sensor vertical offset from midline (inches)
   p2offset = abs(pres2length * cos(ipspitch*pi/180));
   % Total vertical distance between pressure sensor and IPS pressure sensor (meters)
   % (note that ipspitch is neg when nose down)
   % (want ipsdepthauv to be less when vehicle is nose down)
   totdeltadips = sign(ipspitch).*(MVCoffset - (sign(ipspitch).*poffset) + (sign(ipspitch).*p2offset))/39.37; 
   ipsdepthauv = ipsdepth + totdeltadips;
  % calucate total draft
   ipsdraftraw = ipsdepth - ipsrange.*cos(ipspitch*pi/180).*cos(ipsphi*pi/180);
  % determine Ice Profiler timeset year and decimal day from ipstime 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;
  % interpolate depth of vehicle from nav and paros sensor at IPS data points
    auvdepthn=interp1(navtime,ndepth,ipstime);
    auvdepthp=interp1(ptime,parosdpthcalc,ipstime);
  % interpolate estimate of pitch derivative at IPS data points
    derivpitchips = interp1(derivpitchtime,derivpitch,ipstime);
  % convert and store data in netcdf format
   ncoutips;

end


fprintf('\n');
%%-----------------------------------------------------------
%% Define averaging time frame from primary CTD data 
%%-----------------------------------------------------------
% See user defined variables in beginning of file for definition of avgint
%      mtime = [ceil(min(newt1time(1),newt2time(1))):avgint: ...
%              floor(max(newt1time(length(newt1time)),newt2time(length(newt2time))))]';
if ((exist('newt1time') == 1) & (exist('newc2time') == 1))
   if (length(newt1time) > length(newc2time))
      mtime = [ceil(min(newt1time)):avgint:floor(max(newt1time))]';
  else
      mtime = [ceil(min(newc2time)):avgint:floor(max(newc2time))]';
  end      
elseif (exist('newt1time') == 1)
   mtime = [ceil(min(newt1time)):avgint:floor(max(newt1time))]';
elseif (exist('newc2time') == 1)
   mtime = [ceil(min(newc2time)):avgint:floor(max(newc2time))]';
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;



%%-----------------------------------
%% Filter all data
%%-----------------------------------

%% All variables will be averaged to a user defined number of seconds, except lat/lon 
% and nitrogen 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.5 m/s

% divide the averaging time by the timestep of the CTD data acquisition
datafilt = ones(1,avgint/sampint)/(avgint/sampint);
  fprintf(' Time bin is %3u seconds.  Averaging/filtering data ...\n',avgint);

% subsample averaged data
if (exist('newt1time') == 1)
  T1 = interp1(newt1time,filtfilt(datafilt,1,t1),mtime);
  P1 = interp1(newt1time,filtfilt(datafilt,1,p1),mtime);
  C1 = interp1(newt1time,filtfilt(datafilt,1,c1),mtime);
  S1 = interp1(newt1time,filtfilt(datafilt,1,s1),mtime);
  D1 = interp1(newt1time,filtfilt(datafilt,1,d1),mtime);
  OBSbin = interp1(newt1time,filtfilt(datafilt,1,obs),mtime);
  O2 = interp1(newt1time,filtfilt(datafilt,1,o2),mtime);
  FLOW1 = interp1(newt1time,filtfilt(datafilt,1,filtflw1),mtime);
  depth1 = interp1(newt1time,filtfilt(datafilt,1,dpth1),mtime);
end
if (exist('newt2time') == 1)
  T2 = interp1(newt2time,filtfilt(datafilt,1,t2),mtime);
  P2 = interp1(newt2time,filtfilt(datafilt,1,p2),mtime);
  C2 = interp1(newt2time,filtfilt(datafilt,1,c2),mtime);
  S2 = interp1(newt2time,filtfilt(datafilt,1,s2),mtime);
  D2 = interp1(newt2time,filtfilt(datafilt,1,d2),mtime);
  FLOW2 = interp1(newt2time,filtfilt(datafilt,1,filtflw2),mtime);
  depth2 = interp1(newt2time,filtfilt(datafilt,1,dpth2),mtime);
end
if( exist('n2') == 1 )
  N2 = interp1(newt1time,filtfilt(datafilt,1,n2new),mtime);
end
  
% divide the averaging time by the timestep of the navigation data acquisition
  datafiltn = ones(1,avgint/0.2)/(avgint/0.2);
  
% subsample averaged data
  depthnav = interp1(navtime,filtfilt(datafiltn,1,ndepth),mtime);
  pitch = interp1(navtime,filtfilt(datafiltn,1,npitch),mtime);
  roll = interp1(navtime,filtfilt(datafiltn,1,nroll),mtime);
  yaw = interp1(navtime,filtfilt(datafiltn,1,nyaw),mtime);
  Xpos = interp1(navtime,filtfilt(datafiltn,1,nposx),mtime);
  Ypos = interp1(navtime,filtfilt(datafiltn,1,nposy),mtime);

saveIps = 0;
if ( exist('ips4log') == 1 )
if ( isempty(ips4log) ~= 1 )
   if avgint > 1
    % divide the averaging time by the timestep of the IPS data acquisition
      datafiltips = ones(1,avgint/1)/(avgint/1);
  
    % subsample averaged data
      IPSPRES = interp1(ipstime,filtfilt(datafiltips,1,ipspres),mtime);
      IPSDEPTH = interp1(ipstime,filtfilt(datafiltips,1,ipsdepth),mtime);
      IPSDRAFT = interp1(ipstime,filtfilt(datafiltips,1,ipsdraft),mtime);
      IPSDRAFTRAW = interp1(ipstime,filtfilt(datafiltips,1,ipsdraftraw),mtime);
   else
      IPSPRES = interp1(ipstime,ipspres,mtime);
      IPSDEPTH = interp1(ipstime,ipsdepth,mtime);
      IPSDRAFT = interp1(ipstime,ipsdraft,mtime);
      IPSDRAFTRAW = interp1(ipstime,ipsdraftraw,mtime);
   end
 
  saveIps = 1;
end
end
  
  ncoutbinned;
% Plot
% set(gcf,'color',[1 1 1],'units','pixels','position',[10 40 1200 800]);
% plotallNST
  plotsci2
  
  save(DirName,'-mat')
  fprintf('-----DONE PROCESSING----');

  return
