function parseSKQunderway_1min_1sec(CruiseName,DOY,parse,OW,minAve,secAve)
%{
MATLAB script to read R/V Sikuliaq logged underway data. 

*** NOTE: BE PATIENT!  IT TAKES 30 MIN TO PROCESS ONE DAY'S DATA
    TO 1MIN & 1HZ AVERAGES

 Input parameters:
    CruiseName:  UAF cruise name, e.g. 'SKQ201512S'
    DOY: Day-of-year, ('Julian Date'), e.g. 32 for Feb 1
    parse: 1 = parse raw data files, anything else = skip parse and 
               assume parsed files are in current directory
    OW: 1 = overwrite existing files, anything else = Don't overwrite
    minAve: 1 = Compute 1-min averages, anything else = Don't compute
    secAve: 1 = Compute 1Hz averages, anything else = Don't compute

 Example Syntax: parseSKQunderway_1min_1sec('SKQ201512S',274,1,1,1,0);

 Parsed data is written out to one file for each sensor AND a single 
 file that combines all selected sensors' data into one-minute averages.
 In the average file, missing data are filled with NaN's.

 If you want to process the current day's data, copy the /lds/raw/ files
 onto your local computer so you don't try to read from a file that is
 actively being written to.
 
 Edit the Base_Path and/or the StreamsToProcess list to select new sensors
 or data locations.

 Caveats:
   Not all SKQ sensors are currently being parsed out. Modify code below
        to add additional sensors.
   This is coded for a MATLAB on a MAC. For a PC, you will 
        need to reverse the signs of the slashes in the 
        pathnames and replace "/Volumes/" with "Z:\" or whatever 
        letter that you have mapped the /archive drive to.
   This script has not been vetted exhaustively. Please treat
        the data with suspicion until you are convinced that all is
        well. I imagine there still exist some unknown bugs and limitations.

                                            Seth Danielson 03/15/2015

 update 3.19.2015: added new sensors and made code more general for
 different NMEA codes.  Note special syntax of NMEA strings that are 
 handled with multiple passes of one input file (to be made more 
 efficient in future versions). SLD

*** If you want to process the current day's data, copy the /lds/raw/ files
    onto your local computer so you don''t try to read from a file that is
    actively being written to.

    Edit the StreamsToProcess list to select new sensors.
    Edit the Base_Path to select new data locations.


 update 10.15.2015:
    - modified to process one day of data at a time
    - averaging is tolerant of NaNs in the raw input
    - added 1Hz averaging option for wind/motion variables (slow!)
    - function is OS-independent, except for Base_Path definition below
    - time stamp for averaged data indicates start of averaging interval
    - all velocities converted to m/s
    - fluoro_turner variable commented out of this version
    - added following variables to StreamsToProcess: thermo_sbe38_incubator,
        wind_gill_fwdmast, wind_mast_port_true, wind_mast_stbd_true
    - tolerates missing files, inserting NaNs into the 1min/1sec averages
    BWB
%}

tic
%% initialize paths and variables

%Base_Path = ['Z:\',CruiseName,'\lds\raw\'];               % for PC remote mount
Base_Path = ['/Volumes/',CruiseName,'/lds/raw/'];          % for MAC remote mount
%Base_Path = ['~/projects/SKQ/',CruiseName,'/lds/raw/'];   % for MAC local

% Data will be included in the combined datafile based on the order
% they appear within StreamsToProcess
%	unused stream: 'fluoro_turner-c6',... 

 StreamsToProcess = strvcat('ins_seapath_position_$GPRMC',...
                            'ins_seapath_position_$GPHDT',...
                            'ins_seapath_position_$PSXN23',...
                            'mb_em302_centerbeam', ...
                            'met_ptu307' , ...
                            'tsg_sbe45_fwd',...
                            'wind_gill_fwdmast',...
                            'wind_gill_fwdmast_true',...
                            'rad_psp-pir',...
                            'rad_qsr2150a',... 
                            'thermo_pyrometer-ct15',...
                            'flow_krohne_fwd',...
                            'gyro_1_$HEHDT',...
                            'gyro_1_$TIROT',...
                            'speedlog',...
                            'wind_mast_port_true',...
                            'wind_mast_stbd_true');

year = str2num(CruiseName(4:7));
[month,day] = yd2md(year, DOY);
disp(['Processing for Date ',num2str(year),'/',num2str(month),'/',num2str(day)]);
timeStr = [sprintf('%04i',year),sprintf('%02i',month),sprintf('%02i',day)];

if OW == 1
    OW = 'w';
else
    OW = 'a';
end

%% Initialize what we know about these files:

CombinedFORMAT = '%4i %2i %2i %2i %2i %2i ';
NumFields = 6;
VarNames = strvcat('Year','Month','Day','Hour','Minute','Second');
[mP,nP] = size(StreamsToProcess);
Headings = zeros(mP,1);

for F = 1:mP
    DataStream = strtrim(StreamsToProcess(F,:));
    if strcmp(DataStream,'met_ptu307') % as required for each type
        CombinedFORMAT = [CombinedFORMAT '%8.2f %6.1f %6.1f '];
        NumFields = NumFields + 3;
        VarNames = strvcat(VarNames,'met_ptu307 - Atmospheric pressure[hPa]',...
            'met_ptu307 - Air Temperature[C]',...
            'met_ptu307 - Relative Humidity[%]');
        Nvars_met_ptu307 = 9;
        
    elseif ~isempty(findstr(DataStream,'flow_krohne_fwd'))
        CombinedFORMAT = [CombinedFORMAT '%6.3f %6.2f %6.2f %6.5f '];
        NumFields = NumFields + 4;
        VarNames = strvcat(VarNames,'flow_krohne_fwd - flow speed [m/s]',...
            'flow_krohne_fwd - volume flow [l/min]',...
            'flow_krohne_fwd - coil temperature [C]',...
            'flow_krohne_fwd - conductivity [S/m]');
        Nvars_flow_krohne_fwd = 10;
%     elseif ~isempty(findstr(DataStream,'fluoro_turner-c6'))
%           CombinedFORMAT = [CombinedFORMAT '%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f '];
%           NumFields = NumFields + 7;
%           VarNames = strvcat(VarNames,'fluoro_turner-c6 - Phycoeryth(RFU)',...
%                                       'fluoro_turner-c6 - CDOM(RFU)',...
%                                       'fluoro_turner-c6 - Chlorophyll_a(RFU)',...
%                                       'fluoro_turner-c6 - Crude Oil(RFU)',...
%                                       'fluoro_turner-c6 - Turbidity(RFU)',...
%                                       'fluoro_turner-c6 - Depth(m) ',...
%                                       'fluoro_turner-c6 - Temp(C)');
%    Nvars_fluoro_turner-c6 = 13;
    elseif ~isempty(findstr(DataStream,'mb_em302_centerbeam'))
        CombinedFORMAT = [CombinedFORMAT '%10.6f %11.6f %8.2f '];
        NumFields = NumFields + 3;
        VarNames = strvcat(VarNames,'mb_em302_centerbeam - Latitude [Decimal Degrees]',...
            'mb_em302_centerbeam - Longitude [Decimal Degrees]',...
            'mb_em302_centerbeam - Depth [Meters]');
        Nvars_mb_em302_centerbeam = 9;
    elseif strcmp(DataStream,'rad_psp-pir')
        CombinedFORMAT = [CombinedFORMAT '%8.2f %8.2f %8.2f %8.2f %8.2f '];
        NumFields = NumFields + 5;
        VarNames = strvcat(VarNames,'rad_psp-pir - PIR (thermopile, mV)',...
            'rad_psp-pir - LW (computed longwave downwelling irradiance,Wm**2)',...
            'rad_psp-pir - TCASE (PIR case temperature, C)',...
            'rad_psp-pir - TDOME (PIR dome temperature, C)',...
            'rad_psp-pir - SW (computed shortwave downwelling irradiance,Wm**2)');
        Nvars_rad_psp_pir = 11;
    elseif ~isempty(findstr(DataStream,'rad_qsr2150a'))
        CombinedFORMAT = [CombinedFORMAT '%8.2f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'rad_qsr2150a - PAR (uE/m2 sec)');
        Nvars_rad_qsr2150a = 7;
    elseif strcmp(DataStream,'thermo_pyrometer-ct15')
        CombinedFORMAT = [CombinedFORMAT '%8.2f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'thermo_pyrometer-ct15 - Sea Surface skin temperature (deg C)');
        Nvars_thermo_pyrometer_ct15 = 7;
    elseif ~isempty(findstr(DataStream,'tsg_sbe45_fwd'))
        CombinedFORMAT = [CombinedFORMAT '%8.4f %9.5f %8.4f %8.3f '];
        NumFields = NumFields + 4;
        VarNames = strvcat(VarNames,'tsg_sbe45_fwd - Temperature [C]',...
            'tsg_sbe45_fwd - Conductivity [S/m]',...
            'tsg_sbe45_fwd - Salinity [psu]',...
            'tsg_sbe45_fwd - Speed of Sound [m/s]');
        Nvars_tsg_sbe45_fwd = 10;
    elseif ~isempty(findstr(DataStream,'thermo_sbe38_incubator'))
        CombinedFORMAT = [CombinedFORMAT '%8.4f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'thermo_sbe38_incubator - Temperature [C]');
        Nvars_thermo_sbe38_incubator = 7;
    elseif strcmp(DataStream,'wind_gill_fwdmast')
        CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
        NumFields = NumFields + 2;
        VarNames = strvcat(VarNames,'wind_gill_fwdmast - Relative wind direction, 0 to 359 degrees',...
            'wind_gill_fwdmast - Relative wind speed, M[m/s]');
        Nvars_wind_gill_fwdmast = 8;
    elseif strcmp(DataStream,'wind_gill_fwdmast_true')
        CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
        NumFields = NumFields + 2;
        VarNames = strvcat(VarNames,'wind_gill_fwdmast_true - Wind direction, 0 to 359 degrees,T[True]',...
            'wind_gill_fwdmast_true - Wind speed, M[m/s]');
        Nvars_wind_gill_fwdmast_true = 8;
    elseif strcmp(DataStream,'ins_seapath_position_$GPRMC')
        CombinedFORMAT = [CombinedFORMAT '%10.6f %11.6f %5.1f %6.2f '];
        NumFields = NumFields + 4;
        VarNames = strvcat(VarNames,'ins_seapath_position - Latitude [Decimal Degrees]',...
            'ins_seapath_position - Longitude [Decimal Degrees]',...
            'ins_seapath_position - Speed over ground [m/s]',...
            'ins_seapath_position - Course over Ground [Decimal Degrees True]');
        Nvars_ins_seapath_position_GPRMC = 10;
    elseif strcmp(DataStream,'ins_seapath_position_$GPHDT')
        CombinedFORMAT = [CombinedFORMAT '%6.2f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'ins_seapath_position - Heading [Degrees]');
        Nvars_ins_seapath_position_GPHDT = 7;
        Headings(F) = 7;
    elseif strcmp(DataStream,'ins_seapath_position_$PSXN23')
        CombinedFORMAT = [CombinedFORMAT '%6.2f %6.2f %6.2f %6.2f '];
        NumFields = NumFields + 4;
        VarNames = strvcat(VarNames,'ins_seapath_position - Roll [Degrees]');
        VarNames = strvcat(VarNames,'ins_seapath_position - Pitch [Degrees]');
        VarNames = strvcat(VarNames,'ins_seapath_position - Heading [Degrees]');
        VarNames = strvcat(VarNames,'ins_seapath_position - Heave  ');
        Nvars_ins_seapath_position_PSXN23 = 10;
        Headings(F) = 9;
    elseif strcmp(DataStream,'gyro_1_$HEHDT')
        CombinedFORMAT = [CombinedFORMAT '%6.2f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'gyro_1 - Heading [Degrees]');
        Nvars_gyro_1_HEHDT = 7;
        Headings(F) = 7;
    elseif strcmp(DataStream,'gyro_1_$TIROT')
        CombinedFORMAT = [CombinedFORMAT '%7.1f '];
        NumFields = NumFields + 1;
        VarNames = strvcat(VarNames,'gyro_1 - Rate of Turn');
        Nvars_gyro_1_TIROT = 7;
    elseif strcmp(DataStream,'speedlog')
        CombinedFORMAT = [CombinedFORMAT '%7.1f %7.1f %7.1f %7.1f '];
        NumFields = NumFields + 4;
        VarNames = strvcat(VarNames,'speedlog - Longitudinal water speed [m/s]');
        VarNames = strvcat(VarNames,'speedlog - Transverse water speed [m/s]');
        VarNames = strvcat(VarNames,'speedlog - Longitudinal ground speed [m/s]');
        VarNames = strvcat(VarNames,'speedlog - Transverse ground speed [m/s]');
        Nvars_speedlog = 10;
    elseif strcmp(DataStream,'wind_mast_stbd_true')
        CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
        NumFields = NumFields + 2;
        VarNames = strvcat(VarNames,'wind_mast_stbd_true - Wind direction, 0 to 359 degrees,T[True]',...
            'wind_mast_stbd_true - Wind speed, M[m/s]');
        Nvars_wind_mast_stbd_true = 8;
    elseif strcmp(DataStream,'wind_mast_port_true')
        CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
        NumFields = NumFields + 2;
        VarNames = strvcat(VarNames,'wind_mast_port_true - Wind direction, 0 to 359 degrees,T[True]',...
            'wind_mast_port_true - Wind speed, M[m/s]');
        Nvars_wind_mast_port_true = 8;
    end
end

%% Loop through each file type
if parse
    
    for F = 1:mP,
        
        DataStream = strtrim(StreamsToProcess(F,:));
        if ~isempty(strfind(DataStream,'$')),
            DataStreamF = DataStream(1:strfind(DataStream,'$')-2);
        else
            DataStreamF = DataStream;
        end
        
        files = dir(fullfile(Base_Path,DataStreamF,[DataStreamF '*']));
        [mF,nF] = size(files);
        
        for in = 1:mF  % locate file for the specified day
            if findstr(files(in).name,timeStr);
                break;
            end
            if in == mF % end of list reached, file not found
                in = -1;
            end
        end
        
        if in == -1
            disp(['File ',DataStreamF,'.',timeStr,' not found.']);
            continue
        end
        
        fidOUT = fopen([DataStream,'_',sprintf('%03i',DOY),'.dat'],OW);
        fidIN = fopen(fullfile(Base_Path,DataStreamF,files(in).name),'r');
        disp(['Working on ' files(in).name]);
        lineCnt = 0;
        
        while 1
            LOD = fgetl(fidIN);  % get a Line Of Data from file
            lineCnt = lineCnt+1;
            if ~isstr(LOD),
                fclose(fidIN);
                break,
            end
            
            if isempty(LOD)
                beep; pause(0.2); beep;
                disp([' '])
                disp(['Warning: Sensor input file ',files(in).name,' unreadable']),
                disp(['   at line number ',str2num(lineCnt)])
                disp(['Fix file permissions before continuing ']),
                disp(['   or remove this sensor from the input list.']),
                disp(['   Don''t read from a file being actively written to.']),
                disp([' '])
                break,
            end
            
            wTAB = strfind(LOD,'	');
            SkipWriting = 0;
            
            if isempty(strfind(LOD(1),'#')) && ...
                    isempty(strfind(LOD,'metadata test')) && ... % check for header info
                    ~isempty(wTAB),
                DateTime = LOD(wTAB(1)+1:wTAB(2)-1);
                DateTime(strfind(DateTime,':')) = ' ';
                DateTime(strfind(DateTime,'Z')) = ' ';
                DateTime(strfind(DateTime,'T')) = ' ';
                DateTime(strfind(DateTime,'-')) = ' ';
                DT = str2num(DateTime);
                
                if strcmp(DataStream,'met_ptu307') % parse as required for each type
                    wP = strfind(LOD,'P=');
                    wT = strfind(LOD,'T=');
                    wR = strfind(LOD,'RH=');
                    
                    SLP = str2num(LOD(wP+2:wP+10));
                    T = str2num(LOD(wT+2:wT+7));
                    RH = str2num(LOD(wR+3:wR+8));
                    
                    OUTDATA = [DT SLP T RH];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.2f %6.1f %6.1f \r\n';
                    if length(OUTDATA) ~= 9,
                        SkipWriting = 1;
                    end
                elseif strcmp(DataStream,'flow_krohne_fwd')
                    FLOW = str2num(LOD(wTAB(2):end));
                    OUTDATA = [DT FLOW];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %6.3f %6.2f %6.2f %6.5f \r\n';
                    if length(OUTDATA) ~= 10,
                        SkipWriting = 2;
                    end
                elseif strcmp(DataStream,'fluoro_turner-c6')
                    wSPACE = strfind(LOD,' ');
                    FLUOR = str2num(LOD(wSPACE(4):end));
                    OUTDATA = [DT FLUOR];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \r\n';
                    if length(OUTDATA) ~= 13,
                        SkipWriting = 3;
                    end
                elseif strcmp(DataStream,'mb_em302_centerbeam')
                    wCOMMA = strfind(LOD,',');
                    LATLONDEP = str2num(LOD(wCOMMA(5):wCOMMA(8)));
                    OUTDATA = [DT LATLONDEP];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %10.6f %11.6f %8.2f \r\n';
                    if length(OUTDATA) ~= 9,
                        SkipWriting = 4;
                    end
                elseif strcmp(DataStream,'rad_psp-pir')
                    wCOMMA = strfind(LOD,',');
                    RAD = str2num([LOD(wCOMMA(4)+1:wCOMMA(5)-1),' ',LOD(wCOMMA(5)+1:wCOMMA(6)-1),...
                        ' ',LOD(wCOMMA(6)+1:wCOMMA(7)-1),' ',LOD(wCOMMA(7)+1:wCOMMA(8)-1),...
                        ' ',LOD(wCOMMA(8)+1:wCOMMA(9)-1) ]);
                    OUTDATA = [DT RAD];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.2f %8.2f %8.2f %8.2f %8.2f \r\n';
                    if length(OUTDATA) ~= 11,
                        SkipWriting = 6;
                    end
                elseif strcmp(DataStream,'rad_qsr2150a')
                    wCOMMA = strfind(LOD,',');
                    RAD = str2num(LOD(wTAB(2)+1:wCOMMA(1)-1));
                    OUTDATA = [DT RAD];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.2f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 7;
                    end
                elseif strcmp(DataStream,'thermo_pyrometer-ct15')
                    wC = strfind(LOD,'C');
                    wSPACE = strfind(LOD,' ');
                    RAD = str2num(LOD(wC-8:wC-1));
                    OUTDATA = [DT RAD];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.2f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 8;
                    end
                elseif strcmp(DataStream,'tsg_sbe45_fwd')
                    TSG = str2num(LOD(wTAB(2):end));
                    OUTDATA = [DT TSG];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.4f %9.5f %8.4f %8.3f \r\n';
                    if length(OUTDATA) ~= 10,
                        SkipWriting = 9;
                    end
                elseif strcmp(DataStream,'thermo_sbe38_incubator')
                    SST = str2num(LOD(wTAB(2):end));
                    OUTDATA = [DT SST];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %8.4f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 10;
                    end
                elseif strcmp(DataStream,'wind_gill_fwdmast')
                    wCOMMA = strfind(LOD,',');
                    WIND = str2num([ LOD(wCOMMA(1)+1:wCOMMA(2)-1), ' '...
                        LOD(wCOMMA(2)+1:wCOMMA(3)-1)]);
                    if isempty(WIND)  % no data, just commas
                        WIND = [NaN NaN];
                    end
                    OUTDATA = [DT WIND]; % DT DIR SPEED
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %5.1f %6.1f \r\n';
                    if length(OUTDATA) ~= 8,
                        SkipWriting = 11;
                    end
                elseif strcmp(DataStream,'wind_gill_fwdmast_true')
                    wCOMMA = strfind(LOD,',');
                    WIND = str2num([ LOD(wCOMMA(1)+1:wCOMMA(2)-1), ' '...
                        LOD(wCOMMA(7)+1:wCOMMA(8)-1)]);
                    if isempty(WIND)  % no data, just commas
                        WIND = [NaN NaN];
                    end
                    OUTDATA = [DT WIND]; % DT DIR SPEED
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %5.1f %6.1f \r\n';
                    if length(OUTDATA) ~= 8,
                        SkipWriting = 12;
                    end
                elseif ~isempty(findstr(StreamsToProcess(F,:),'$GPRMC')) && ~isempty(findstr(LOD,'$GPRMC'))
                    wCOMMA = strfind(LOD,',');
                    LA = str2num(LOD(wCOMMA(3)+1:wCOMMA(4)-1));
                    LO = str2num(LOD(wCOMMA(5)+1:wCOMMA(6)-1));Y(1) = floor(LA/100);
                    Y(2) = LA-Y(1)*100;
                    wCOMMA = strfind(LOD,',');
                    LA = str2num(LOD(wCOMMA(3)+1:wCOMMA(4)-1));
                    LO = str2num(LOD(wCOMMA(5)+1:wCOMMA(6)-1));
                    X(1) = floor(LO/100);
                    X(2) = LO-X(1)*100;
                    X = X(1) + X(2)/60;
                    Y = Y(1) + Y(2)/60;
                    NS = LOD(wCOMMA(4)+1:wCOMMA(5)-1);
                    EW = LOD(wCOMMA(6)+1:wCOMMA(7)-1);
                    if NS == 'S',
                        Y = -Y;
                    end
                    if EW == 'W',
                        X = -X;
                    end
                    SOG = str2num(LOD(wCOMMA(7)+1:wCOMMA(8)-1))*0.514; % in m/s
                    COG = str2num(LOD(wCOMMA(8)+1:wCOMMA(9)-1));
                    OUTDATA = [DT Y X SOG COG];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %10.6f %11.6f %5.1f %6.2f \r\n';
                    if length(OUTDATA) ~= 10,
                        SkipWriting = 13;
                    end
                elseif ~isempty(findstr(StreamsToProcess(F,:),'$GPHDT')) && ~isempty(findstr(LOD,'$GPHDT'))
                    wCOMMA = strfind(LOD,',');
                    HDT = str2num(LOD(wCOMMA(1)+1:wCOMMA(2)-1));
                    OUTDATA = [DT HDT];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %6.2f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 14;
                    end
                elseif ~isempty(findstr(StreamsToProcess(F,:),'$PSXN23')) && ~isempty(findstr(LOD,'$PSXN,23'))
                    wCOMMA = strfind(LOD,',');
                    wSTAR = strfind(LOD,'*');
                    RPHH = str2num(LOD(wCOMMA(2)+1:wSTAR-1));
                    OUTDATA = [DT RPHH];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %6.2f %6.2f %6.2f %6.2f \r\n';
                    if length(OUTDATA) ~= 10,
                        SkipWriting = 15;
                    end
                elseif ~isempty(findstr(StreamsToProcess(F,:),'$HEHDT')) && ~isempty(findstr(LOD,'$HEHDT'))
                    wCOMMA = strfind(LOD,',');
                    HDT = str2num(LOD(wCOMMA(1)+1:wCOMMA(2)-1));
                    OUTDATA = [DT HDT];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %6.2f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 16;
                    end
                elseif ~isempty(findstr(StreamsToProcess(F,:),'$TIROT')) && ~isempty(findstr(LOD,'$TIROT'))
                    wCOMMA = strfind(LOD,',');
                    HDT = str2num(LOD(wCOMMA(1)+1:wCOMMA(2)-1));
                    OUTDATA = [DT HDT];
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %7.1f \r\n';
                    if length(OUTDATA) ~= 7,
                        SkipWriting = 17;
                    end
                elseif ~isempty(findstr(DataStream,'speedlog'))  && ~isempty(findstr(LOD,'$VDVBW'))
                    try
                        wCOMMA = strfind(LOD,',');
                        SLOG1 = str2num(LOD(wCOMMA(1)+1:wCOMMA(3)-1))*0.514; % in m/s
                        SLOG2 = str2num(LOD(wCOMMA(4)+1:wCOMMA(6)-1))*0.514; % in m/s
                        if length(SLOG1) ~= 2 || LOD(wCOMMA(3)+1:wCOMMA(4)-1) ~= 'A',
                            SLOG1 = [NaN NaN];
                        end
                        if length(SLOG2) ~= 2 || LOD(wCOMMA(6)+1:wCOMMA(7)-1) ~= 'A',
                            SLOG2 = [NaN NaN];
                        end
                        OUTDATA = [DT SLOG1 SLOG2];
                        OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %4.1f %4.1f %4.1f %4.1f \r\n';
                        if length(OUTDATA) ~= 10,
                            SkipWriting = 18;
                        end
                    catch
                        SkipWriting = 18;
                    end
                elseif strcmp(DataStream,'wind_mast_stbd_true')
                    wCOMMA = strfind(LOD,',');
                    WIND = str2num([ LOD(wCOMMA(1)+1:wCOMMA(2)-1), ' '...
                        LOD(wCOMMA(7)+1:wCOMMA(8)-1)]);
                    if isempty(WIND)  % no data, just commas
                        WIND = [NaN NaN];
                    end
                    OUTDATA = [DT WIND]; % DT DIR SPEED
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %5.1f %6.1f \r\n';
                    if length(OUTDATA) ~= 8,
                        SkipWriting = 19;
                    end
                elseif strcmp(DataStream,'wind_mast_port_true')
                    wCOMMA = strfind(LOD,',');
                    WIND = str2num([ LOD(wCOMMA(1)+1:wCOMMA(2)-1), ' '...
                        LOD(wCOMMA(7)+1:wCOMMA(8)-1)]);
                    if isempty(WIND)  % no data, just commas
                        WIND = [NaN NaN];
                    end
                    OUTDATA = [DT WIND]; % DT DIR SPEED
                    OUTFORMAT = '%4i %2i %2i %2i %2i %6.3f %5.1f %6.1f \r\n';
                    if length(OUTDATA) ~= 8,
                        SkipWriting = 20;
                    end
                else
                    SkipWriting = 99;
                end
                if ~SkipWriting,
                    fprintf(fidOUT,OUTFORMAT,OUTDATA);
                elseif SkipWriting < 99,
                    disp(['Found non-conforming data: ',DataStream,':line ',num2str(lineCnt)]);
                end
            end
        end
    end
    
    fclose all;
    
end

%% Combine data from all files and compute 1-min averages

if minAve

    [mF,nF] = size(StreamsToProcess);
    AllTimes = [];
    
    for F = 1:mF
        disp(F)
        DataStream = strtrim(StreamsToProcess(F,:));
        if ~exist([DataStream,'_',sprintf('%03i',DOY),'.dat'],'file');
            disp(['File ',DataStream,'.dat not found, Setting values to NaN']);
            Dstr = DataStream; Dstr(strfind(DataStream,'-')) = '_';
            Dstr = DataStream; Dstr(strfind(DataStream,'$')) = '';
            eval(['Ncols = Nvars_',Dstr,';'])             % number of columns for the missing file
            eval(['D',num2str(F),' = NaN(86400,Ncols);']) % dummay array of 1 Hz NaNs
            % default 1 Hz timestamp
            eval(['Dtime',num2str(F),' = datenum(year,month,day):1/86400:datenum(year,month,day)+1-(1/86400);']);
            eval(['AllTimes = Dtime',num2str(F),';'])
        else
            eval(['load ',DataStream,'_',sprintf('%03i',DOY),'.dat']);
            DataStream(strfind(DataStream,'-')) = '_';
            DataStream(strfind(DataStream,'$')) = '_';
            eval(['D',num2str(F),' = ',DataStream,'_',sprintf('%03i',DOY),';'])
            eval(['[m,n] = size(D',num2str(F),');']);
            eval(['Dtime',num2str(F),' = datenum(D',num2str(F),'(:,1:6));'])
            eval(['AllTimes = Dtime',num2str(F),';'])
            eval(['clear ',DataStream,'_',sprintf('%03i',DOY)])
            
            % don't average directions - convert magnitude/direction to vector format
            if ~isempty(strfind(DataStream,'wind_gill_fwdmast'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Relative wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                % relative to ship's bow = 0 deg
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            elseif ~isempty(strfind(DataStream,'wind_gill_fwdmast_true'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            elseif ~isempty(strfind(DataStream,'GPRMC'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,9);              % SOG
                DIR = DD(:,10);             % COG
                U = MAG.*sin(DIR*pi/180);   % direction is where the boat is GOING TOWARD
                V = MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,9) = U;']);
                eval(['D' num2str(F) '(:,10) = V;']);
            elseif ~isempty(strfind(DataStream,'wind_mast_stbd_true'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            elseif ~isempty(strfind(DataStream,'wind_mast_port_true'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            end
        end
    end
    
    startTime = round(min(floor(AllTimes*24*60)))/24/60;
    endTime = round(max(ceil(AllTimes*24*60))-1)/24/60;
    DL = length(startTime:1/24/60:endTime);
    fidOUT = fopen([CruiseName,'_underway_',sprintf('%03i',DOY),'_1min.txt'],OW);
    
    [m,n] = size(VarNames);
    fprintf(fidOUT,'%s\r\n',['% R/V Sikuliaq Underway Data: 1 min averages from ' CruiseName]);
    fprintf(fidOUT,'%s\r\n',['% Time stamp at start of averaging interval']);
    for in = 1:m
        fprintf(fidOUT,'%s\r\n',['% Column #' num2str(in) ': ' VarNames(in,:)]);
    end
    
    disp(['Averaging and Exporting data to File: ',CruiseName,'_underway_',sprintf('%03i',DOY),'_1min.txt']);
    for in = startTime:1/24/60:endTime
        AveDATA = str2num(datestr(in,'yyyy mm dd HH MM SS'));
        for F = 1:mF
            DataStream = strtrim(StreamsToProcess(F,:));
            eval(['w = find(Dtime' num2str(F) ' >= in & ' ...
                'Dtime' num2str(F) ' <  in + 1/24/60);'])
            
            if length(w) > 1
                eval(['Ave =  nanmean1(D' num2str(F) '(w,:));']);
                AveDATA = [AveDATA Ave(7:end)];
            elseif length(w) == 1,
                eval(['Ave =  D' num2str(F) '(w,:);']);
                AveDATA = [AveDATA Ave(7:end)];
            else
                eval(['[m,n] = size(D' num2str(F) ');']);
                AveDATA = [AveDATA NaN*ones(1,n-6)];
            end
            if ~isempty(strfind(DataStream,'wind_gill_fwdmast'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            elseif ~isempty(strfind(DataStream,'wind_gill_fwdmast_true'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            elseif ~isempty(strfind(DataStream,'GPRMC'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi; % direction is where the boat is GOING TOWARD
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [MAG DIR];
            elseif ~isempty(strfind(DataStream,'wind_mast_stbd_true'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            elseif ~isempty(strfind(DataStream,'wind_mast_port_true'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            end
            
            % special averaging for directions to account for 0/360 phase wraparound.
            % This ignores speed, so it is just a time-weighted average, NOT a vetor-averaged direction.
            
            if Headings(F) > 0,
                
                eval(['[m,n] = size(D' num2str(F) ');']);
                Nback = n-Headings(F);
                if length(w) > 1
                    eval(['A =  D' num2str(F) '(w,Headings(F));']);
                    n0 = find(A <= 90 | A > 270);
                    n180 = find(A > 90 & A <= 270);
                    if length(n0) > length(n180)
                        wQ1 = find(A < 180 & A >= 0);
                        A(wQ1) = A(wQ1) + 360;
                    end
                    Ave = mean(A);
                    if Ave > 360,
                        Ave = Ave - 360;
                    end
                    AveDATA(end-Nback) = Ave; % place average heading in appropriate spot
                elseif length(w) == 1,
                    eval(['Ave =  D' num2str(F) '(w,Headings(F));']);
                    AveDATA(end-Nback) = Ave; % place average heading in appropriate spot
                end
            end
        end
        if length(AveDATA) == NumFields,
            fprintf(fidOUT,[CombinedFORMAT '\r\n'],AveDATA);
        end
    end
    fclose all;
end

%% Do 1-sec averages

if secAve
    
    % define subset of variables for 1Hz averaging
     StreamsToProcess = strvcat('ins_seapath_position_$GPRMC',...
                                'ins_seapath_position_$GPHDT',...
                                'ins_seapath_position_$PSXN23',...
                                'wind_gill_fwdmast',...
                                'wind_gill_fwdmast_true',...
                                'gyro_1_$HEHDT',...
                                'gyro_1_$TIROT',...
                                'speedlog');
    
    % Initialize what we know about these files, again...
    CombinedFORMAT = '%4i %2i %2i %2i %2i %2i ';
    NumFields = 6;
    VarNames = strvcat('Year','Month','Day','Hour','Minute','Second');
    [mP,nP] = size(StreamsToProcess);
    Headings = zeros(mP,1);

    for F = 1:mP
        DataStream = strtrim(StreamsToProcess(F,:));
        if strcmp(DataStream,'wind_gill_fwdmast')
            CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
            NumFields = NumFields + 2;
            VarNames = strvcat(VarNames,'wind_gill_fwdmast - Relative wind direction, 0 to 359 degrees',...
                'wind_gill_fwdmast - Relative wind speed, M[m/s]');
            Nvars_wind_gill_fwdmast = 8;
        elseif strcmp(DataStream,'wind_gill_fwdmast_true')
            CombinedFORMAT = [CombinedFORMAT '%6.1f %5.1f '];
            NumFields = NumFields + 2;
            VarNames = strvcat(VarNames,'wind_gill_fwdmast_true - Wind direction, 0 to 359 degrees,T[True]',...
                'wind_gill_fwdmast_true - Wind speed, M[m/s]');
            Nvars_wind_gill_fwdmast_true = 8;
        elseif strcmp(DataStream,'ins_seapath_position_$GPRMC')
            CombinedFORMAT = [CombinedFORMAT '%10.6f %11.6f %5.1f %6.2f '];
            NumFields = NumFields + 4;
            VarNames = strvcat(VarNames,'ins_seapath_position - Latitude [Decimal Degrees]',...
                'ins_seapath_position - Longitude [Decimal Degrees]',...
                'ins_seapath_position - Speed over ground [m/s]',...
                'ins_seapath_position - Course over Ground [Decimal Degrees True]');
            Nvars_ins_seapath_position_GPRMC = 10;
        elseif strcmp(DataStream,'ins_seapath_position_$GPHDT')
            CombinedFORMAT = [CombinedFORMAT '%6.2f '];
            NumFields = NumFields + 1;
            VarNames = strvcat(VarNames,'ins_seapath_position - Heading [Degrees]');
            Nvars_ins_seapath_position_GPHDT = 7;
            Headings(F) = 7;
        elseif strcmp(DataStream,'ins_seapath_position_$PSXN23')
            CombinedFORMAT = [CombinedFORMAT '%6.2f %6.2f %6.2f %6.2f '];
            NumFields = NumFields + 4;
            VarNames = strvcat(VarNames,'ins_seapath_position - Roll [Degrees]');
            VarNames = strvcat(VarNames,'ins_seapath_position - Pitch [Degrees]');
            VarNames = strvcat(VarNames,'ins_seapath_position - Heading [Degrees]');
            VarNames = strvcat(VarNames,'ins_seapath_position - Heave  ');
            Nvars_ins_seapath_position_PSXN23 = 10;
            Headings(F) = 9;
        elseif strcmp(DataStream,'gyro_1_$HEHDT')
            CombinedFORMAT = [CombinedFORMAT '%6.2f '];
            NumFields = NumFields + 1;
            VarNames = strvcat(VarNames,'gyro_1 - Heading [Degrees]');
            Nvars_gyro_1_HEHDT = 7;
            Headings(F) = 7;
        elseif strcmp(DataStream,'gyro_1_$TIROT')
            CombinedFORMAT = [CombinedFORMAT '%7.1f '];
            NumFields = NumFields + 1;
            VarNames = strvcat(VarNames,'gyro_1 - Rate of Turn');
            Nvars_gyro_1_TIROT = 7;
        elseif strcmp(DataStream,'speedlog')
            CombinedFORMAT = [CombinedFORMAT '%7.1f %7.1f %7.1f %7.1f '];
            NumFields = NumFields + 4;
            VarNames = strvcat(VarNames,'speedlog - Longitudinal water speed [m/s]');
            VarNames = strvcat(VarNames,'speedlog - Transverse water speed [m/s]');
            VarNames = strvcat(VarNames,'speedlog - Longitudinal ground speed [m/s]');
            VarNames = strvcat(VarNames,'speedlog - Transverse ground speed [m/s]');
            Nvars_speedlog = 10;
        end
    end
    
    [mF,nF] = size(StreamsToProcess);
    AllTimes = [];
    
    for F = 1:mF
        disp(F)
        DataStream = strtrim(StreamsToProcess(F,:));
        if ~exist([DataStream,'_',sprintf('%03i',DOY),'.dat'],'file');
            disp(['File ',DataStream,'.dat not found, Setting values to NaN']);
            Dstr = DataStream; Dstr(strfind(DataStream,'-')) = '_';
            Dstr = DataStream; Dstr(strfind(DataStream,'$')) = '';
            eval(['Ncols = Nvars_',Dstr,';'])               % number of columns for the missing file
            eval(['D',num2str(F),' = NaN(86400*2,Ncols);']) % dummay array of 10 Hz NaNs
            % default 2 Hz timestamp
            eval(['Dtime',num2str(F),' = datenum(year,month,day):1/864000:datenum(year,month,day)+1-(1/(86400*2));']);
            eval(['AllTimes = Dtime',num2str(F),';'])
        else
            eval(['load ',DataStream,'_',sprintf('%03i',DOY),'.dat']);
            DataStream(strfind(DataStream,'-')) = '_';
            DataStream(strfind(DataStream,'$')) = '_';
            eval(['D',num2str(F),' = ',DataStream,'_',sprintf('%03i',DOY),';'])
            eval(['[m,n] = size(D',num2str(F),');']);
            eval(['Dtime',num2str(F),' = datenum(D',num2str(F),'(:,1:6));'])
            eval(['AllTimes = Dtime',num2str(F),';'])
            eval(['clear ',DataStream,'_',sprintf('%03i',DOY)])
            
            % don't average directions - convert magnitude/direction to vector format
            if ~isempty(strfind(DataStream,'wind_gill_fwdmast'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Relative wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                % relative to ship's bow = 0 deg
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            elseif ~isempty(strfind(DataStream,'wind_gill_fwdmast_true'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,8);              % Wind Speed
                DIR = DD(:,7);              % direction is where the wind is COMING FROM
                U = -MAG.*sin(DIR*pi/180);
                V = -MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,7) = U;']);
                eval(['D' num2str(F) '(:,8) = V;']);
            elseif ~isempty(strfind(DataStream,'GPRMC'))
                eval(['DD = D' num2str(F) ';']);
                MAG = DD(:,9);              % SOG
                DIR = DD(:,10);             % COG
                U = MAG.*sin(DIR*pi/180);   % direction is where the boat is GOING TOWARD
                V = MAG.*cos(DIR*pi/180);
                eval(['D' num2str(F) '(:,9) = U;']);
                eval(['D' num2str(F) '(:,10) = V;']);
            end
        end
    end
    
    startTime = round(min(floor(AllTimes*24*60*60)))/24/60/60;
    endTime = round(max(ceil(AllTimes*24*60*60))-1)/24/60/60;
    DL = length(startTime:1/24/60/60:endTime);
    fidOUT = fopen([CruiseName,'_underway_',sprintf('%03i',DOY),'_1sec.txt'],OW);
    
    [m,n] = size(VarNames);
    fprintf(fidOUT,'%s\r\n',['% R/V Sikuliaq Underway Data: 1 sec averages from ' CruiseName]);
    fprintf(fidOUT,'%s\r\n',['% Time stamp at start of averaging interval']);
    for in = 1:m
        fprintf(fidOUT,'%s\r\n',['% Column #' num2str(in) ': ' VarNames(in,:)]);
    end
    
    disp(['Averaging and Exporting data to File: ',CruiseName,'_underway_',sprintf('%03i',DOY),'_1sec.txt']);
    for in = startTime:1/24/60/60:endTime
        AveDATA = str2num(datestr(in,'yyyy mm dd HH MM SS'));
        for F = 1:mF
            DataStream = strtrim(StreamsToProcess(F,:));
            eval(['w = find(Dtime' num2str(F) ' >= in & ' ...
                'Dtime' num2str(F) ' <  in + 1/24/60/60);'])
            
            if length(w) > 1
                eval(['Ave =  nanmean1(D' num2str(F) '(w,:));']);
                AveDATA = [AveDATA Ave(7:end)];
            elseif length(w) == 1,
                eval(['Ave =  D' num2str(F) '(w,:);']);
                AveDATA = [AveDATA Ave(7:end)];
            else
                eval(['[m,n] = size(D' num2str(F) ');']);
                AveDATA = [AveDATA NaN*ones(1,n-6)];
            end
            if ~isempty(strfind(DataStream,'wind_gill_fwdmast'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            elseif ~isempty(strfind(DataStream,'wind_gill_fwdmast_true'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi-180; % direction is where the wind is COMING FROM
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [DIR MAG];
            elseif ~isempty(strfind(DataStream,'GPRMC'))
                U = AveDATA(end-1);
                V = AveDATA(end);
                MAG = sqrt(U.^2 + V.^2);
                DIR = 90-atan2(V,U)*180/pi; % direction is where the boat is GOING TOWARD
                wx = find(DIR < 0);
                DIR(wx) = DIR(wx)+360;
                AveDATA(end-1:end) = [MAG DIR];
            end
            
            % special averaging for directions to account for 0/360 phase wraparound.
            % This ignores speed, so it is just a time-weighted average, NOT a vetor-averaged direction.
            
            if Headings(F) > 0,
                
                eval(['[m,n] = size(D' num2str(F) ');']);
                Nback = n-Headings(F);
                if length(w) > 1
                    eval(['A =  D' num2str(F) '(w,Headings(F));']);
                    n0 = find(A <= 90 | A > 270);
                    n180 = find(A > 90 & A <= 270);
                    if length(n0) > length(n180)
                        wQ1 = find(A < 180 & A >= 0);
                        A(wQ1) = A(wQ1) + 360;
                    end
                    Ave = mean(A);
                    if Ave > 360,
                        Ave = Ave - 360;
                    end
                    AveDATA(end-Nback) = Ave; % place average heading in appropriate spot
                elseif length(w) == 1,
                    eval(['Ave =  D' num2str(F) '(w,Headings(F));']);
                    AveDATA(end-Nback) = Ave; % place average heading in appropriate spot
                end
            end
        end
        if length(AveDATA) == NumFields,
            fprintf(fidOUT,[CombinedFORMAT '\r\n'],AveDATA);
        end
    end
    
    fclose all;
end

toc
end

