function sal_statistics(theSite)
%
%   sal_statistics(theSite)
%
%   Function to create tables of the mean, standard deviation, min, and max
%   of salinity 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), 'a')
    surfFiles = {'5391ts-a.nc'; '5392ts-a.nc'; '5393ts-a.nc'; ...
        '5394ts-a.nc'; '5395ts-a.nc'};
    surfDepths = [10 20 30 40 50];
    outFile = 'sal_stat_a.xls';
elseif strcmp(lower(theSite), 'b')
    subFiles = {'5383sc-a.cdf'};
    subDepths = [66]; 
    triFiles = {'5371tcp-a.nc'};
    triDepths = [83];
    outFile = 'sal_stat_b.xls';
elseif strcmp(lower(theSite), 'c')
    surfFiles = {'5331mc-a.cdf'};
    surfDepths = [10];
    subFiles = {'5341mc-a.cdf'; '5342mc-a.cdf'; '5343sc-a.cdf'};
    subDepths = [21 31 40];   
    triFiles = {'5351tcp-a.nc'};
    triDepths = [49];
    outFile = 'sal_stat_c.xls';
elseif strcmp(lower(theSite), 'd/mblt-b')
    triFiles = {'5321sc-a_d1.nc'};
    triDepths = [21];
    outFile = 'sal_stat_d.xls';
elseif strcmp(lower(theSite), 'mblt-a')
    subFiles = {'5291mc-a_d1.nc'; '5293sc-a_d1.nc'};
    subDepths = [5 21];   
    triFiles = {'5301Atcp-a.nc'};
    triDepths = [30];
    outFile = 'sal_stat_e.xls';
end
   
fID = fopen(fullfile(plotDirectory, outFile), 'w');
fprintf(fID, ['Site ' upper(theSite) ', Salinity (ppt)\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\tMean\tSt. Dev.\tMin\tMax\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'}(:));
        theSal = ncID{'S_40'}(:);
        if isempty(theSal)
            theSal = ncID{'salinity'}(:);
        end
        if isempty(theSal)
            theSal = ncID{'sal'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theSal = theSal(inTime);
        
        bads = find(theSal > 1e3 | theSal < -1e3);
        theSal(bads) = NaN;
        
        numPoints = length(find(~isnan(theSal)));
        theMean = nanmean(theSal);
        theSD = nanstd(theSal);
        theMin = nanmin(theSal);
        theMax = nanmax(theSal);
        startIn = min(find(~isnan(theSal)));
        endIn = max(find(~isnan(theSal)));
        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(theMean, '%.2f') '\t' ...
            num2str(theSD, '%.2f') '\t' ...
            num2str(theMin, '%.2f') '\t' ...
            num2str(theMax, '%.2f') '\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\tMean\tSt. Dev.\tMin\tMax\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'}(:));
        theSal = ncID{'S_40'}(:);
        if isempty(theSal)
            theSal = ncID{'salinity'}(:);
        end
        if isempty(theSal)
            theSal = ncID{'sal'}(:);
        end

        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theSal = theSal(inTime);
        
        bads = find(theSal > 1e3 | theSal < -1e3);
        theSal(bads) = NaN;
        
        numPoints = length(find(~isnan(theSal)));
        theMean = nanmean(theSal);
        theSD = nanstd(theSal);
        theMin = nanmin(theSal);
        theMax = nanmax(theSal);
        startIn = min(find(~isnan(theSal)));
        endIn = max(find(~isnan(theSal)));
        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(theMean, '%.2f') '\t' ...
            num2str(theSD, '%.2f') '\t' ...
            num2str(theMin, '%.2f') '\t' ...
            num2str(theMax, '%.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\tMean\tSt. Dev.\tMin\tMax\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'}(:));
        theSal = ncID{'S_40'}(:);
        if isempty(theSal)
            theSal = ncID{'salinity'}(:);
        end
        if isempty(theSal)
            theSal = ncID{'sal'}(:);
        end

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

        bads = find(theSal > 1e3 | theSal < -1e3);
        theSal(bads) = NaN;
        
        numPoints = length(find(~isnan(theSal)));
        theMean = nanmean(theSal);
        theSD = nanstd(theSal);
        theMin = nanmin(theSal);
        theMax = nanmax(theSal);  
        startIn = min(find(~isnan(theSal)));
        endIn = max(find(~isnan(theSal)));
        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(theMean, '%.2f') '\t' ...
            num2str(theSD, '%.2f') '\t' ...
            num2str(theMin, '%.2f') '\t' ...
            num2str(theMax, '%.2f') '\n']);
    end
    fprintf(fID, '\n');
end
fclose(fID)