function density_grad(theSite)
%
%   density_grad(theSite)
%
%   Function to calculate Sigma-T (and the difference between Sigma-T at
%   measured depths) for the MBIWE.

%Soupy Alexander, 6/18/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';

if strcmp(lower(theSite), 'a')
    theFiles = {'5391ts-a.nc'; '5392ts-a.nc'; '5393ts-a.nc'; ...
        '5394ts-a.nc'; '5395ts-a.nc'};
elseif strcmp(lower(theSite), 'c')
    theFiles = {'5331mc-a.cdf'; '5341mc-a.cdf'; '5342mc-a.cdf'; ...
        '5343sc-a.cdf'; '5351tcp-a.nc'};
end

timeRange = julian(1998, 8, 1, 00):1/24:julian(1998, 10, 1, 00); 

for indexFile = 1:length(theFiles)
    ncID = netcdf(fullfile(dataDir, theFiles{indexFile}), 'nowrite');
    
    theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    theTemp = ncID{'T_20'}(:);
    theSal = ncID{'S_40'}(:);
    theDepths(indexFile) = ncID{'depth'}(:);
    theDepthList{indexFile} = [num2str(ncID{'depth'}(:), '%.0f') ' mbs'];
    
    if isempty(theTemp)
        theTemp = ncID{'temperature'}(:);
        theSal = ncID{'salinity'}(:);
    end
    
    bads = find((theTemp + theSal) > 1e3 | (theTemp + theSal) < -1e3);
    theTime(bads) = [];
    theTemp(bads) = [];
    theSal(bads) = [];
    theDensity = sw_dens0(theSal, theTemp);
    
    for indexTime = 1:length(timeRange)
        inTime = find((theTime > (timeRange(indexTime)- 0.5/24)) & ...
            (theTime <= (timeRange(indexTime) + 0.5/24)));
        if isempty(inTime)
            hourDensity(indexTime) = NaN;
        else
            hourDensity(indexTime) = nanmean(theDensity(inTime));
        end
    end
    
    allDensity{indexFile} = hourDensity;
    theStarts(indexFile) = min(find(~isnan(hourDensity)));
    theEnds(indexFile) = max(find(~isnan(hourDensity)));
end

for indexFile = 1:(length(theFiles) - 1)
    firstDepth = theDepths(indexFile);
    secondDepth = theDepths(indexFile + 1);
    theStr = [num2str(secondDepth, '%.0f') '-' ...
        num2str(firstDepth, '%.0f') ' mbs'];
    theDepthList2{indexFile} = theStr;
end

colList = {'r'; 'k'; 'g'; 'm'; 'b'};

%Hourly Average
figure, orient landscape, hold on

for indexFile = 1:length(theFiles)
    
    hitDensity = allDensity{indexFile}-1000;
    plot(julian2datenum(timeRange), hitDensity, colList{indexFile})
    
end

legend(theDepthList, 'location', 'northeastoutside')
xlim([julian2datenum(min(timeRange(theStarts))) ...
    julian2datenum(max(timeRange(theEnds)))])
ylim([22.0 25.5])
set(gca, 'xtick', julian2datenum(min(timeRange)):4:...
    julian2datenum(max(timeRange)))
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('sigma-t (kg/m^3)')
title(['Site ' upper(theSite) ', Hourly Averaged Sigma-t'])
box('on')

%Hourly Average (Difference)
figure, orient landscape, hold on

for indexFile = 1:(length(theFiles)-1)
    
    hitDensity = allDensity{indexFile}-1000;
    hitDensity2 = allDensity{indexFile+1} - 1000;
    diffDensity = hitDensity2 - hitDensity;
    plot(julian2datenum(timeRange), diffDensity, colList{indexFile})
    
end

legend(theDepthList2, 'location', 'northeastoutside')
xlim([min(julian2datenum(timeRange(theStarts))) ...
    max(julian2datenum(timeRange(theEnds)))])
ylim([0 2.5])
set(gca, 'xtick', julian2datenum(min(timeRange)):4:...
    julian2datenum(max(timeRange)))
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('sigma-t (kg/m^3)')
title(['Site ' upper(theSite) ', Hourly Averaged Density Difference'])
box('on')

%Low-passed
figure, orient landscape, hold on

for indexFile = 1:length(theFiles)
    
    hitDensity = allDensity{indexFile}-1000;
    [lpDensity, lpTime] = plfilt(hitDensity, timeRange, 6);
    plot(julian2datenum(lpTime), lpDensity, colList{indexFile})
    
end

legend(theDepthList, 'location', 'northeastoutside')
xlim([julian2datenum(min(timeRange(theStarts))) ...
    julian2datenum(max(timeRange(theEnds)))])
ylim([22 25.5])
set(gca, 'xtick', julian2datenum(min(timeRange)):4:...
    julian2datenum(max(timeRange)))
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('sigma-t (kg/m^3)')
title(['Site ' upper(theSite) ', Low-passed Sigma-t'])
box('on')

%Low-passed difference
figure, orient landscape, hold on

for indexFile = 1:(length(theFiles)-1)
    
    hitDensity = allDensity{indexFile}-1000;
    [lpDensity, lpTime] = plfilt(hitDensity, timeRange, 6);
    hitDensity2 = allDensity{indexFile+1} - 1000;
    [lpDensity2, lpTime2] = plfilt(hitDensity2, timeRange, 6);
    diffDensity = lpDensity2 - lpDensity;
    plot(julian2datenum(lpTime), diffDensity, colList{indexFile})
    
end

legend(theDepthList2, 'location', 'northeastoutside')
xlim([julian2datenum(min(timeRange(theStarts))) ...
    julian2datenum(max(timeRange(theEnds)))])
ylim([0 2.5])
set(gca, 'xtick', julian2datenum(min(timeRange)):4:...
    julian2datenum(max(timeRange)))
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('sigma-t (kg/m^3)')
title(['Site ' upper(theSite) ', Low-passed Density Difference'])
box('on')