function principal_axes 

%   principal_axes
%
%   Function to calculate the principal axes for the MBIWE and output a text file of the results.

%Soupy Alexander, 3/17/04

%This Matlab m-file was used to create portions of U.S. Geological Survey 
%Data Series 85. Although this program has been used by the USGS, no 
%warranty, expressed or implied, is made by the USGS or the United States 
%Government as to the accuracy and functioning of the program and related 
%program material nor shall the fact of distribution constitute any such 
%warranty, and no responsibility is assumed by the USGS in connection 
%therewith.

%Citation:  Butman, Bradford, Alexander, P. Soupy, Anderson, S.P., 
%Lightsom, F.L., Scotti, Alberto, and Beardsley, R.C., 2004, The 
%Massachusetts Bay Internal Wave Experiment, August 1998:  Data Report:  
%U.S. Geological Survey Data Series 85, 1 DVD-ROM.

dataDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';

theSites = {'B'; 'B'; 'B'; 'B'; 'B'; 'B'; ...
        'B'; 'B'; 'B'; 'C'; 'C'; 'D/MBLT-B'; 'D/MBLT-B'};
theFiles = {'fsi1467-a.nc'; 'fsi1468-a.nc'; 'fsi1469-a.nc'; ...
        'vmcm0014.epic'; 'vmcm0009.epic'; 'fsi1428-a.nc'; ...
        '5372adc-a.nc'; '5371v-a_d1.nc'; '5371v-a_d2.nc'; ...
        '5344-a.cdf'; '5352adc-a.nc';
        '5322adc-a.nc'; '5322adc-a.nc'};
theInstruments = {'FSI 1467'; 'FSI 1468'; 'FSI 1469'; ...
            'VMCM 0014'; 'VMCM 0009'; 'FSI 1428'; ...
            'ADCP 5372'; 'BASS 5371'; 'BASS 5371'; ...
            'VMCM 5344'; 'ADCP 5352'; ...
            'ADCP 5322'; 'ADCP 5322'};
theDepths = {'20 mbs'; '20 mbs'; '20 mbs'; ...
            '50 mbs'; '75 mbs'; '80 mbs'; ...
            '76 mbs'; '84 mbs'; '84.6 mbs'; ...
            '45 mbs'; '45 mbs'; '5 mbs'; '16 mbs'};

printDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\';

fID = fopen(fullfile(printDir, 'pr_axes.xls'), 'w');
fprintf(fID, 'Site\tInstrument\tDepth\tStart Time\tEnd Time\tSampling Interval\tMajor Axis\tMinor Axis\tOrientation\tEllipticity\n');


for indexFile = 1:length(theFiles)
    ncID = netcdf(fullfile(dataDir, theFiles{indexFile}), 'nowrite');
    theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    theDepth = ncID{'depth'}(:);
    if indexFile ~= 12
        hitDepth = find(theDepth == max(theDepth));
    else
        hitDepth = value2Index(theDepth, 5);
    end
    
    theU = ncID{'u_1205'}(:,hitDepth);
    theV = ncID{'v_1206'}(:,hitDepth);
    if isempty(theU)
        theU = ncID{'east'}(:);
        theV = ncID{'north'}(:);
    end
    bads = find(theU > 1e3 | theU < -1e3 | theV > 1e3 | theV < -1e3);
    theU(bads) = [];
    theV(bads) = [];   
    theTime(bads) = [];
    
    meanU = nanmean(theU);
    meanV = nanmean(theV);
    [majAx, majAz, minAx, minAz, theEllip] = pcaben(theU, theV);
    
    fprintf(fID, [theSites{indexFile} '\t' theInstruments{indexFile} '\t' ...
            theDepths{indexFile} '\t' datestr(julian2datenum(min(theTime)),0) '\t' ...
            datestr(julian2datenum(max(theTime)),0) '\t' ...
            num2str(median(diff(theTime))*24*60, '%.1f') '\t' num2str(majAx, '%.2f') '\t' ...
            num2str(minAx, '%.2f') '\t' num2str(majAz, '%.0f') '\t' ...
            num2str(theEllip, '%.2f') '\n']);
    
end

fclose('all')