function curr_statistics(theSite)
%
%   curr_statistics(theSite)
%
%   Function to create tables of the mean, standard deviation, min, and max
%   of current for each site of the Mass Bay Internal Wave Experiment.

%Soupy Alexander, 6/14/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.

startTime = [1998 08 06 12 00 00];
endTime = [1998 08 30 12 00 00];

dataDirectory = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';
plotDirectory = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT';

if strcmp(lower(theSite), 'b')
    hmFiles = {{'fsi1467-a.nc'}; ...
        {'fsi1468-a.nc'}; ...
        {'fsi1469-a.nc'}};
    hmDepths = [20];
    subFiles = {'vmcm0014-a.nc'; ...
        'vmcm0009-a.nc'; 'fsi1428-a.nc'};
    subDepths = [50 75 80]; 
    triFiles = {'5371v-a_d1.nc'; '5371v-a_d2.nc'};
    triDepths = [84 85];
    outFile = 'curr_stat_b.xls';
elseif strcmp(lower(theSite), 'c')
    subFiles = {'5344-a.cdf'};
    subDepths = [45];   
    outFile = 'curr_stat_c.xls';
elseif strcmp(lower(theSite), 'mblt-a')
    subFiles = {'5294-a_d1.nc'};
    subDepths = [21];   
    triFiles = {'5301Av-a_d1.nc'};
    triDepths = [30];
    outFile = 'curr_stat_e.xls';
end
   
fID = fopen(fullfile(plotDirectory, outFile), 'w');
fprintf(fID, ['Site ' upper(theSite) ', Current (cm/s)\n\n']);

if exist('surfFiles', 'var')
    
    fprintf(fID, 'Surface Moorings\n');
    fprintf(fID, 'Depth (m)\tFile\tStart Time\tEnd Time\tSampling Interval (mins)\t# Points\tU Mean\tU STD\tV Mean\tV STD\tVec. Speed\tVec. Dir.\tSc. Mean\tSc. STD\tSc. Min\tSc. Max\n');
    
    for indexFile = 1:length(surfFiles)
        hitFile = surfFiles{indexFile};
        [junk, fileName, junk] = fileparts(hitFile);
        ncID = netcdf(fullfile(dataDirectory, hitFile), 'nowrite');
        theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
        theU = ncID{'u_1205'}(:);
        theV = ncID{'v_1206'}(:);
        if isempty(theU)
            theU = ncID{'east'}(:);
            theV = ncID{'north'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theU = theU(inTime);
        theV = theV(inTime);
        
        bads = find((theU + theV) > 1e3 | (theU + theV) < -1e3);
        theU(bads) = NaN;
        theV(bads) = NaN;
        
        numPoints = length(find(~isnan((theU + theV))));
        uMean = nanmean(theU);
        uSD = nanstd(theU);
        vMean = nanmean(theV);
        vSD = nanstd(theV);
        
        [vecTheta, vecSpeed] = cart2pol(uMean, vMean);
        vecDir = anglextocompass(vecTheta*180/pi);
        
        [scalTheta, scalSpeed] = cart2pol(theU, theV);
        scalDir = anglextocompass(scalTheta*180/pi);
        scalMean = nanmean(scalSpeed);
        scalSD = nanstd(scalSpeed);        
        scalMin = nanmin(scalSpeed);
        scalMax = nanmax(scalSpeed);
        startIn = min(find(~isnan(scalSpeed)));
        endIn = max(find(~isnan(scalSpeed)));
        startDate = datestr(julian2datenum(theTime(startIn)), 0);
        endDate = datestr(julian2datenum(theTime(endIn)), 0);
        sampleInt = median(diff(theTime))*24*60;
        
        fprintf(fID, [num2str(surfDepths(indexFile), '%.0f') '\t' ...
            fileName '\t' startDate '\t' endDate '\t' ...
            num2str(sampleInt, '%.2f') '\t' ...
            num2str(numPoints, '%.0f') '\t' ...
            num2str(uMean, '%.2f') '\t' ...
            num2str(uSD, '%.2f') '\t' ...    
            num2str(vMean, '%.2f') '\t' ...
            num2str(vSD, '%.2f') '\t' ...
            num2str(vecSpeed, '%.2f') '\t' ...
            num2str(vecDir, '%.0f') '\t' ...
            num2str(scalMean, '%.2f') '\t' ...
            num2str(scalSD, '%.2f') '\t' ...
            num2str(scalMin, '%.2f') '\t' ...
            num2str(scalMax, '%.2f') '\n']);
    end
    fprintf(fID, '\n');
end

if exist('hmFiles', 'var')
    
    fprintf(fID, 'Horizontal Array\n');
    fprintf(fID, 'Depth (m)\t');

    fprintf(fID, repmat('File\tStart Time\tEnd Time\tSampling Interval (mins)\t# Points\tU Mean\tU STD\tV Mean\tV STD\tVec. Speed\tVec. Dir.\tSc. Mean\tSc. STD\tSc. Min\tSc. Max\t', ...
        1, 1));
    fprintf(fID, '\n');
    
    for indexString = 1:length(hmFiles)
        thisString = hmFiles{indexString};
        for indexFile = 1:length(thisString)
            hitFile = thisString{indexFile};
            if isempty(hitFile)
                continue
            else   
                fprintf(fID, [num2str(hmDepths(indexFile), '%.0f') '\t']);
                [junk, fileName, junk] = fileparts(hitFile);
                ncID = netcdf(fullfile(dataDirectory, hitFile), 'nowrite');
        theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
        theU = ncID{'u_1205'}(:);
        theV = ncID{'v_1206'}(:);
        if isempty(theU)
            theU = ncID{'east'}(:);
            theV = ncID{'north'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theU = theU(inTime);
        theV = theV(inTime);
        
        bads = find((theU + theV) > 1e3 | (theU + theV) < -1e3);
        theU(bads) = NaN;
        theV(bads) = NaN;
        
        numPoints = length(find(~isnan((theU + theV))));
        uMean = nanmean(theU);
        uSD = nanstd(theU);
        vMean = nanmean(theV);
        vSD = nanstd(theV);
        
        [vecTheta, vecSpeed] = cart2pol(uMean, vMean);
        vecDir = anglextocompass(vecTheta*180/pi);
        
        [scalTheta, scalSpeed] = cart2pol(theU, theV);
        scalDir = anglextocompass(scalTheta*180/pi);
        scalMean = nanmean(scalSpeed);
        scalSD = nanstd(scalSpeed);        
        scalMin = nanmin(scalSpeed);
        scalMax = nanmax(scalSpeed);
        startIn = min(find(~isnan(scalSpeed)));
        endIn = max(find(~isnan(scalSpeed)));
        startDate = datestr(julian2datenum(theTime(startIn)), 0);
        endDate = datestr(julian2datenum(theTime(endIn)), 0);
        sampleInt = median(diff(theTime))*24*60;
        
        fprintf(fID, [fileName '\t' ...
            startDate '\t' endDate '\t' num2str(sampleInt, '%.2f') '\t' ...
            num2str(numPoints, '%.0f') '\t' ...
            num2str(uMean, '%.2f') '\t' ...
            num2str(uSD, '%.2f') '\t' ...    
            num2str(vMean, '%.2f') '\t' ...
            num2str(vSD, '%.2f') '\t' ...
            num2str(vecSpeed, '%.2f') '\t' ...
            num2str(vecDir, '%.0f') '\t' ...
            num2str(scalMean, '%.2f') '\t' ...
            num2str(scalSD, '%.2f') '\t' ...
            num2str(scalMin, '%.2f') '\t' ...
            num2str(scalMax, '%.2f') '\t']);
            end
        end
        fprintf(fID, '\n');
    end
    fprintf(fID, '\n');
end

if exist('subFiles', 'var')
    
    fprintf(fID, 'Subsurface Moorings\n');
    fprintf(fID, 'Depth (m)\tFile\tStart Time\tEnd Time\tSampling Interval (mins)\t# Points\tU Mean\tU STD\tV Mean\tV STD\tVec. Speed\tVec. Dir.\tSc. Mean\tSc. STD\tSc. Min\tSc. Max\n');
    
    for indexFile = 1:length(subFiles)
        hitFile = subFiles{indexFile};
        [junk, fileName, junk] = fileparts(hitFile);
        ncID = netcdf(fullfile(dataDirectory, hitFile), 'nowrite');
        theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
        theU = ncID{'u_1205'}(:);
        theV = ncID{'v_1206'}(:);
        if isempty(theU)
            theU = ncID{'east'}(:);
            theV = ncID{'north'}(:);
        end

        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theU = theU(inTime);
        theV = theV(inTime);

        bads = find((theU + theV) > 1e3 | (theU + theV) < -1e3);
        theU(bads) = NaN;
        theV(bads) = NaN;

        numPoints = length(find(~isnan((theU + theV))));
        uMean = nanmean(theU);
        uSD = nanstd(theU);
        vMean = nanmean(theV);
        vSD = nanstd(theV);
        
        [vecTheta, vecSpeed] = cart2pol(uMean, vMean);
        vecDir = anglextocompass(vecTheta*180/pi);
        
        [scalTheta, scalSpeed] = cart2pol(theU, theV);
        scalDir = anglextocompass(scalTheta*180/pi);
        scalMean = nanmean(scalSpeed);
        scalSD = nanstd(scalSpeed);        
        scalMin = nanmin(scalSpeed);
        scalMax = nanmax(scalSpeed);  
        startIn = min(find(~isnan(scalSpeed)));
        endIn = max(find(~isnan(scalSpeed)));
        startDate = datestr(julian2datenum(theTime(startIn)), 0);
        endDate = datestr(julian2datenum(theTime(endIn)), 0);
        sampleInt = median(diff(theTime))*24*60;
        
        fprintf(fID, [num2str(subDepths(indexFile), '%.0f') '\t' ...
            fileName '\t' startDate '\t' endDate '\t' ...
            num2str(sampleInt, '%.2f') '\t' ...
            num2str(numPoints, '%.0f') '\t' ...
            num2str(uMean, '%.2f') '\t' ...
            num2str(uSD, '%.2f') '\t' ...    
            num2str(vMean, '%.2f') '\t' ...
            num2str(vSD, '%.2f') '\t' ...
            num2str(vecSpeed, '%.2f') '\t' ...
            num2str(vecDir, '%.0f') '\t' ...
            num2str(scalMean, '%.2f') '\t' ...
            num2str(scalSD, '%.2f') '\t' ...
            num2str(scalMin, '%.2f') '\t' ...
            num2str(scalMax, '%.2f') '\n']);
    end
    fprintf(fID, '\n');
end

if exist('triFiles', 'var')
    
    fprintf(fID, 'Tripod\n');
    fprintf(fID, 'Depth (m)\tFile\tStart Time\tEnd Time\tSampling Interval (mins)\t# Points\tU Mean\tU STD\tV Mean\tV STD\tVec. Speed\tVec. Dir.\tSc. Mean\tSc. STD\tSc. Min\tSc. Max\n');
    
    for indexFile = 1:length(triFiles)
        hitFile = triFiles{indexFile};
        [junk, fileName, junk] = fileparts(hitFile);
        ncID = netcdf(fullfile(dataDirectory, hitFile), 'nowrite');
        theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
        theU = ncID{'u_1205'}(:);
        theV = ncID{'v_1206'}(:);
        if isempty(theU)
            theU = ncID{'east'}(:);
            theV = ncID{'north'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theU = theU(inTime);
        theV = theV(inTime);

        bads = find((theU + theV) > 1e3 | (theU + theV) < -1e3);
        theU(bads) = NaN;
        theV(bads) = NaN;
        
        numPoints = length(find(~isnan((theU + theV))));
        uMean = nanmean(theU);
        uSD = nanstd(theU);
        vMean = nanmean(theV);
        vSD = nanstd(theV);
        
        [vecTheta, vecSpeed] = cart2pol(uMean, vMean);
        vecDir = anglextocompass(vecTheta*180/pi);
        
        [scalTheta, scalSpeed] = cart2pol(theU, theV);
        scalDir = anglextocompass(scalTheta*180/pi);
        scalMean = nanmean(scalSpeed);
        scalSD = nanstd(scalSpeed);        
        scalMin = nanmin(scalSpeed);
        scalMax = nanmax(scalSpeed);  
        startIn = min(find(~isnan(scalSpeed)));
        endIn = max(find(~isnan(scalSpeed)));
        startDate = datestr(julian2datenum(theTime(startIn)), 0);
        endDate = datestr(julian2datenum(theTime(endIn)), 0);
        sampleInt = median(diff(theTime))*24*60;
        
        fprintf(fID, [num2str(triDepths(indexFile), '%.0f') '\t' ...
            fileName  '\t' startDate '\t' endDate '\t' ...
            num2str(sampleInt, '%.2f') '\t' ...
            num2str(numPoints, '%.0f') '\t' ...
            num2str(uMean, '%.2f') '\t' ...
            num2str(uSD, '%.2f') '\t' ...    
            num2str(vMean, '%.2f') '\t' ...
            num2str(vSD, '%.2f') '\t' ...
            num2str(vecSpeed, '%.2f') '\t' ...
            num2str(vecDir, '%.0f') '\t' ...
            num2str(scalMean, '%.2f') '\t' ...
            num2str(scalSD, '%.2f') '\t' ...
            num2str(scalMin, '%.2f') '\t' ...
            num2str(scalMax, '%.2f') '\n']);
    end
    fprintf(fID, '\n');
end
fclose(fID)