function temp_statistics(theSite)
%
%   temp_statistics(theSite)
%
%   Function to create tables of the mean, standard deviation, min, and max
%   of temperature 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 = 'temp_stat_a.xls';
elseif strcmp(lower(theSite), 'b')
    surfFiles = {'tpod4494-a.nc'; 'tpod4491-a.nc'; 'tpod4483-a.nc'; ...
        'tpod3662-a.nc'};
    surfDepths = [0 10 20 30];
    hmFiles = {{'sbe0035-a.nc'; 'sbe0052-a.nc'; 'sbe0045-a.nc'; ...
        'sbe0048-a.nc'; 'sbe0044-a.nc'; 'tpod8093-a.nc'}; ...
        {'fsi1467-a.nc'}; ...
        {'sbe0040-a.nc'; 'sbe0007-a.nc'; 'sbe0054-a.nc'; ...
        'sbe0046-a.nc'; 'sbe0050-a.nc'; 'tpod8094-a.nc'}; ...
        {'fsi1468-a.nc'}; ...
        {'mcat0683-a.nc'; []; 'mcat0686-a.nc'; ...
        'mcat0685-a.nc'; []; 'mcat0684-a.nc'}; ...
        {'sbe0049-a.nc'; 'sbe0051-a.nc'; []; ...
        'sbe0039-a.nc'; 'sbe0101-a.nc'; 'tpod8091-a.nc'}; ...
        {'fsi1469-a.nc'}; ...
        {'sbe0041-a.nc'; 'sbe0038-a.nc'; 'sbe0047-a.nc'; ...
        'sbe0103-a.nc'; []; 'tpod8090-a.nc'}; ...
        {'mcat0670-a.nc'}};
    hmDepths = [20 25 30 35 40 45];
    subFiles = {'vmcm0014-a.nc'; 'wadar62-a.nc'; '5383sc-a.cdf'; ...
        'vmcm0009-a.nc'; 'fsi1428-a.nc'; 'tide0046s-a.nc'};
    subDepths = [50 50 66 75 80 85]; 
    triFiles = {'5371tcp-a.nc'};
    triDepths = [83];
    outFile = 'temp_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 = 'temp_stat_c.xls';
elseif strcmp(lower(theSite), 'd/mblt-b')
    triFiles = {'5321sc-a_d1.nc'};
    triDepths = [21];
    outFile = 'temp_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 = 'temp_stat_e.xls';
end
   
fID = fopen(fullfile(plotDirectory, outFile), 'w');
fprintf(fID, ['Site ' upper(theSite) ', Temperature (°C)\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'}(:));
        theTemp = ncID{'T_20'}(:);
        if isempty(theTemp)
            theTemp = ncID{'temperature'}(:);
        end
        if isempty(theTemp)
            theTemp = ncID{'temp'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theTemp = theTemp(inTime);
        
        bads = find(theTemp > 1e3 | theTemp < -1e3);
        theTemp(bads) = NaN;
        
        numPoints = length(find(~isnan(theTemp)));
        theMean = nanmean(theTemp);
        theSD = nanstd(theTemp);
        theMin = nanmin(theTemp);
        theMax = nanmax(theTemp);
        startIn = min(find(~isnan(theTemp)));
        endIn = max(find(~isnan(theTemp)));
        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('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\tMean\tSt. Dev.\tMin\tMax\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'}(:));
                theTemp = ncID{'T_20'}(:);
                if isempty(theTemp)
                    theTemp = ncID{'temperature'}(:);
                end
                if isempty(theTemp)
                    theTemp = ncID{'temp'}(:);
                end

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

                bads = find(theTemp > 1e3 | theTemp < -1e3);
                theTemp(bads) = NaN;

                numPoints = length(find(~isnan(theTemp)));
                theMean = nanmean(theTemp);
                theSD = nanstd(theTemp);
                theMin = nanmin(theTemp);
                theMax = nanmax(theTemp);
                startIn = min(find(~isnan(theTemp)));
                endIn = max(find(~isnan(theTemp)));
                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(theMean, '%.2f') '\t' ...
                    num2str(theSD, '%.2f') '\t' ...
                    num2str(theMin, '%.2f') '\t' ...
                    num2str(theMax, '%.2f') '\n']);
            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\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'}(:));
        theTemp = ncID{'T_20'}(:);
        if isempty(theTemp)
            theTemp = ncID{'temperature'}(:);
        end
        if isempty(theTemp)
            theTemp = ncID{'temp'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theTemp = theTemp(inTime);

        bads = find(theTemp > 1e3 | theTemp < -1e3);
        theTemp(bads) = NaN;
        
        numPoints = length(find(~isnan(theTemp)));
        theMean = nanmean(theTemp);
        theSD = nanstd(theTemp);
        theMin = nanmin(theTemp);
        theMax = nanmax(theTemp);
        startIn = min(find(~isnan(theTemp)));
        endIn = max(find(~isnan(theTemp)));
        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'}(:));
        theTemp = ncID{'T_20'}(:);
        if isempty(theTemp)
            theTemp = ncID{'temperature'}(:);
        end
        if isempty(theTemp)
            theTemp = ncID{'temp'}(:);
        end
        
        inTime = find((theTime >= julian(startTime)) & ...
            (theTime <= julian(endTime)));
        theTime = theTime(inTime);
        theTemp = theTemp(inTime);

        bads = find(theTemp > 1e3 | theTemp < -1e3);
        theTemp(bads) = NaN;
        
        numPoints = length(find(~isnan(theTemp)));
        theMean = nanmean(theTemp);
        theSD = nanstd(theTemp);
        theMin = nanmin(theTemp);
        theMax = nanmax(theTemp);
        startIn = min(find(~isnan(theTemp)));
        endIn = max(find(~isnan(theTemp)));
        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)