function scatter_plots

%   scatter_plots
%
%   Function to create scatter plots of east/north current for the MBIWE.

%Soupy Alexander, 3/16/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'; 'C'; 'D/MBLT-B'};
theFiles = {{{'fsi1467-a.nc'; 'fsi1468-a.nc'; 'fsi1469-a.nc'}; ...
        {'vmcm0014-a.nc'; 'vmcm0009-a.nc'; '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'}}};
theTitles = {{{'FSI 1467, 20 mbs'; 'FSI 1468, 20 mbs'; 'FSI 1469, 20 mbs'}; ...
            {'VMCM 0014, 50 mbs'; 'VMCM 0009, 75 mbs'; 'FSI 1428, 80 mbs'}; ...
            {'ADCP 5372, 76 mbs'; 'BASS 5371, 84 mbs'; 'BASS 5371, 84.6 mbs'}}; ...
        {{'VMCM 5344, 45 mbs'}; {'ADCP 5352, 45 mbs'}}; ...
        {{'ADCP 5322, 5 mbs'; 'ADCP 5322, 16 mbs'}}};
superTitles = {{'Site B, Horizontal Array'; 'Site B, Subsurface'; 'Site B, Tripod'}; ...
        {'Site C, Subsurface'; 'Site C, Tripod'}; {'Site D/MBLT-B, ADCP'}};

printDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\';
printNames = {'scat_b_horiz'; 'scat_b_sub'; 'scat_b_tri'; ...
        'scat_c_sub'; 'scat_c_tri'; 'scat_d_adcp'};
   
%Calculate the principal ellipses and means.
theta=2*pi*(1:64)/64;
theta=[0 theta];
    
for indexSite = 1 %1:length(theSites)
    siteFiles = theFiles{indexSite};
    siteTitles = theTitles{indexSite};
    siteSuper = superTitles{indexSite};
    
    for indexSet = 2%:length(siteFiles)
        setFiles = siteFiles{indexSet};
        setTitles = siteTitles{indexSet};
        setSuper = siteSuper{indexSet};
        
        figure, orient tall
        
        for indexFile = 1:length(setFiles)
            ncID = netcdf(fullfile(dataDir, setFiles{indexFile}), 'nowrite');
            theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
            theDepth = ncID{'depth'}(:);
            if indexSite == 3 & indexFile ==1
                hitDepth = value2Index(theDepth, 5);
            else
                hitDepth = find(theDepth == max(theDepth));
            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);
            [x1, y1] = cmgspd2uv(majAx, majAz);
            [x2, y2] = cmgspd2uv(minAx, minAz);
            xx = majAx*cos(theta);
            yy = minAx*sin(theta);
            theAngle = -1*majAz*pi/180+pi/2;
            xxx = xx.*cos(theAngle) - yy.*sin(theAngle);
            yyy = xx.*sin(theAngle) + yy.*cos(theAngle);
            
            subplot(length(setFiles), 2, indexFile*2-1)
            thePosition = get(gca, 'position');
            plot(theU, theV, '.')
            hold on
            plot([0 meanU], [0 meanV], 'r', 'linewidth', 2);
            plot(xxx+meanU, yyy+meanV, 'r', 'linewidth', 2);     
            set(gca, 'fontsize', 8)
            xlabel('Eastward Component (cm/s)')
            ylabel('Northward Component (cm/s)')
            axis equal, axis tight   
            axis([-50 50 -50 50])
            labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
            l = line([-50 50], [0 0]);
            l2 = line([0 0], [-50 50]);
            set(l,'color','k', 'linestyle','-.');
            set(l2,'color','k', 'linestyle','-.');
            set(gca,'xtick',labels);
            set(gca,'ytick',labels);
            title({setTitles{indexFile}; 'Basic Sampling Interval'})
            set(gca, 'position', thePosition)
            
            %Hourly averaged
            timeStep = median(diff(theTime));
            newTime = min(theTime):timeStep:max(theTime);     
            inHour = round((1/24)/median(diff(theTime)));
            leftOver = rem(length(newTime), inHour);
            newTime(end-leftOver+1:end) = [];
            newU = smart_interp(theTime, theU, newTime, 0.5);
            newV = smart_interp(theTime, theV, newTime, 0.5);

            newU2 = reshape(newU, inHour, length(newU)/inHour);
            newV2 = reshape(newV, inHour, length(newV)/inHour);
            newTime2 = reshape(newTime, inHour, length(newTime)/inHour);
            
            hourTime = mean(newTime2);
            hourU = mean(newU2);
            hourV = mean(newV2);
            
            %Low-pass
            [lowpassU, lowpassTime] = plfilt(hourU, hourTime, 6);
            [lowpassV] = plfilt(hourV, hourTime, 6);
            meanU = nanmean(lowpassU);
            meanV = nanmean(lowpassV);
            [majAx, majAz, minAx, minAz, theEllip] = pcaben(lowpassU, lowpassV);
            [x1, y1] = cmgspd2uv(majAx, majAz);
            [x2, y2] = cmgspd2uv(minAx, minAz);
            xx = majAx*cos(theta);
            yy = minAx*sin(theta);
            theAngle = -1*majAz*pi/180+pi/2;
            xxx = xx.*cos(theAngle) - yy.*sin(theAngle);
            yyy = xx.*sin(theAngle) + yy.*cos(theAngle);
            
            subplot(length(setFiles), 2, indexFile*2)
            thePosition = get(gca, 'position');
            plot(lowpassU, lowpassV, '.')
            hold on
            plot([0 meanU], [0 meanV], 'r', 'linewidth', 2);
            plot(xxx+meanU, yyy+meanV, 'r', 'linewidth', 2);     
            set(gca, 'fontsize', 8)
            xlabel('Eastward Component (cm/s)')
            ylabel('Northward Component (cm/s)')
            axis equal, axis tight   
            axis([-25 25 -25 25])
            labels = [-25 -20 -15 -10 -5 0 5 10 15 20 25];
            l = line([-25 25], [0 0]);
            l2 = line([0 0], [-25 25]);
            set(l,'color','k', 'linestyle','-.');
            set(l2,'color','k', 'linestyle','-.');
            set(gca,'xtick',labels);
            set(gca,'ytick',labels);    
            title({setTitles{indexFile}; 'Low-Passed'})
            set(gca, 'position', thePosition)
        end
        suptitle(setSuper)
        print('-dpdf', fullfile(printDir, printNames{indexSet}))
    end
end
            
close all