function ctd_bvfq
%
%   ctd_bvfq
%
%   Function to create plots of MBIWE brunt-vaisala frequency.

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


theExt = '_f';
%theExt = '_i';

titleStart = {'Brunt-Vaisala Frequency'};
theLimits = {[-1e-4:1e-4:4e-3]};
theLimits2 = {[-1e-4:5e-4:4e-3]};
plotCom = {'bvfq'};
thePick = [5e-4];

stationNames = {'A'; 'B'; 'C'; 'D'};
stationLats = [42.3315 42.2932 42.2470 42.1620];
stationLons = [-70.3913 -70.4574 -70.5530 -70.6400];

%Loading in/extracting the data
[ctdData, ctdData2, ctdData3] = mbay_iwave_ctd;
theCasts = fieldnames(ctdData);
for index = 1:length(theCasts)
    eval(['hitCTD = ctdData.' theCasts{index} ';']);
    theLons(index) = -1*hitCTD(1,3);
    theLats(index) = hitCTD(1,2);
    theTimes(index) = hitCTD(1,1);
    theDepths(index) = max(hitCTD(:,4));
end

%Load in the pressure time data
presFile = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\OTHER\tide_times.xls';
[presNums, presText] = xlsread(presFile);
presText(1:4,:) = [];
tideTimes = datenum2julian(datenum(presText(:,1))); 

pivotLat = theLats(27);
pivotLon = theLons(27);

depthGrid = 5:1:95;

%Loading in the bathymetry/creating the mask for depth
[bathyData] = xlsread('C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES\ctddepth.xls');
bathyLats = bathyData(:,4);
bathyLons = bathyData(:,3);
theBathy = abs(bathyData(:,5));        
for indexBathy = 1:length(theBathy)
    bathy_dist(indexBathy) = latlon2dist(pivotLat, pivotLon, ...
        bathyLats(indexBathy), bathyLons(indexBathy));
end
bathy_km =  bathy_dist./1000;
theX2s = [min(bathy_km); min(bathy_km); bathy_km(:); ...
        max(bathy_km); max(bathy_km); min(bathy_km)];
theXs = [-0.8; -0.8; bathy_km(:); 24.99; 24.99; -0.8];
theY2s = [99; theBathy(find(bathy_km == min(bathy_km))); theBathy(:); ...
        theBathy(find(bathy_km == max(bathy_km))); 99; 99];
theYs = [99; (theBathy(find(bathy_km == min(bathy_km)))-2); theBathy(:)-2; ...
        (theBathy(find(bathy_km == max(bathy_km)))-2); 99; 99];

section_1 = fliplr([1:9]);
section_2 = fliplr([10:18]);
section_3 = fliplr([19:27]);
section_4 = fliplr([41 40 39 42 37 36 35 34 33 32 31 30 29 28]);
section_5 = fliplr([52 51 50 49 48 47 46 45]);
section_6 = fliplr([53:55]);
section_7 = fliplr([61 60 59 58 57 56 55]);
section_8 = fliplr([63 64 65 66 67]);

for indexVar = 1
    theStep = mean(diff(theLimits{indexVar}));
    theCasts = fieldnames(ctdData);
    for index = 1:length(theCasts)
        eval(['hitCTD = ctdData.' theCasts{index} ';']);
        theLons(index) = -1*abs(hitCTD(1,3));
        theLats(index) = hitCTD(1,2);
        theTimes(index) = hitCTD(1,1);
        theDepths(index) = max(hitCTD(:,4));
    end
    theRun = {'section_1'; 'section_2'; 'section_3'; 'section_4'; 'section_5'; ...
            'section_6'; 'section_7'; 'section_8'};
    for indexRun = 8%1:length(theRun)
        eval(['target = ' theRun{indexRun} ';']); 
        clear theDist theGrid
        for index = 1:length(target)
            indexCast = target(index);
            eval(['hitCTD = ctdData.' theCasts{indexCast} ';']);
            [junk, theGoods] = unique(hitCTD(:,4));
            bads = find(hitCTD(theGoods,4) < 5);
            theGoods(bads) = [];
            theSalinity = hitCTD(theGoods, 6);
            thePressure = hitCTD(theGoods, 7);
            theTemperature = hitCTD(theGoods, 8);
            theLat = hitCTD(1,2);
            theTarg = sw_bfrq(theSalinity,theTemperature,thePressure,theLat);
            %theTarg(find(theTarg <= 0)) = NaN;
            hitDepths = filter(ones(1,2)/2,1,hitCTD(theGoods, 4));
            hitDepths(1) = [];
            theInterp = interp1(hitDepths, theTarg, depthGrid);
            if indexRun ~=8
                lastNonNan = max(find(~isnan(theInterp)));
                theInterp(lastNonNan:end) = theInterp(max(find(~isnan(theInterp))));
            end
            theGrid(:,index) = theInterp(:);
            theDist(index) = latlon2dist(pivotLat, pivotLon, theLats(indexCast), theLons(indexCast));
        end 
        sectionNum = theRun{indexRun};
        theDash = strfind(sectionNum, '_');
        theNum = sectionNum(theDash+1:end);
         theDiffs = tideTimes-min(theTimes(target));
        thisEvent = max(find(theDiffs < 0));
        startAfter = julian2datenum(min(theTimes(target)) - tideTimes(thisEvent));
        endAfter = julian2datenum(max(theTimes(target)) - tideTimes(thisEvent));
        
        
        theTitle = {[titleStart{indexVar} ', Argo Maine Transect ' theNum ': ' ...
                datestr(julian2datenum(min(theTimes(target))),1) ' ' ...
                datestr(julian2datenum(min(theTimes(target))),15) ' to ' ...
                datestr(julian2datenum(max(theTimes(target))),1) ' '...
                datestr(julian2datenum(max(theTimes(target))),15)]; ...
                [datestr(startAfter, 15) ' after low water to ' ...
                datestr(endAfter, 15) ' after low water']};
        dist_km = theDist./1000;
        figure, orient landscape
        if ~strcmp(sectionNum, 'section_8')
            if strcmp(theExt, '_f')
                [c, h] = contourf(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar});
                shading flat
            else
                [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits2{indexVar}, 'k');
            end 
            %clabel(c, h, [min(theLimits{indexVar}):thePick(indexVar):max(theLimits{indexVar})], ...
            %    'labelspacing', 500)  
            p = patch(theXs, theYs, 'w');
            set(p, 'linestyle', 'none')
            p2 = patch(theX2s, theY2s, [0.7 0.7 0.7]);
            set(p2, 'linestyle', 'none');
            for index = 1:length(target)
                l = line([dist_km(index) dist_km(index)], [0 100]);
                set(l, 'color', 'k')
                t = text(dist_km(index)+0.25, 90, [num2str(target(index), '%02.0f')]);
                set(t, 'color', 'k')
            end
            axis([-1 25 0 100])
            xlabel('Distance to Ref. Point (Station 27, SW Corner)') 
            hold on
            plot(bathy_km, theBathy, 'k', 'linewidth', 2)
        else
            if strcmp(theExt, '_f')
                [c, h] = contourf(theTimes(section_8)', depthGrid(:), theGrid, theLimits{indexVar});
                shading flat
            else
                [c,h] = contour(theTimes(section_8)', depthGrid(:), theGrid, theLimits2{indexVar}, ...
                    'k');
            end 
            %clabel(c, h, [min(theLimits{indexVar}):thePick(indexVar):max(theLimits{indexVar})], ...
             %   'labelspacing', 500)
            for index = 1:length(target)
                l = line([theTimes(section_8((index))) theTimes(section_8((index)))], [0 100]);
                set(l, 'color', 'k')
                t = text(theTimes(section_8((index)))+0.001, 90, [num2str(target(index), '%02.0f')]);
                set(t, 'color', 'k')
            end
            ylim([0 100])
            xlim([julian(1998, 9, 2, 13) julian([1998 9 2 15 30 00])])  
            xlabel('Time')
            datetick('x', 15, 'keeplimits')
        end
        for indexStation = 1:length(stationNames)
            st_dist(indexStation) = latlon2dist(pivotLat, pivotLon, ...
                stationLats(indexStation), stationLons(indexStation));
            st_km(indexStation) = st_dist(indexStation)/1000;
            t = text(st_km(indexStation), -2, stationNames{indexStation});
            set(t, 'color', 'k')
        end       
        w = text(0, 92, '<----West');
        set(w, 'horizontala', 'left');
        e = text(24.5, 92, 'East----->');
        set(e, 'horizontala', 'right');
        axis ij
        set(gca, 'clim', [min(theLimits{indexVar}) max(theLimits{indexVar})])
        if strcmp(theExt, '_f')
            colorbar
        end
        t = title(theTitle);
        currentPos = get(t, 'position');
        currentPos(2) = -4;
        set(t, 'position', currentPos)
        ylabel('Depth (m)')
        box('on')
        print(gcf, '-dpdf', ['argo_' theNum '_' plotCom{indexVar} theExt])
        close
    end
    
    if indexVar ~=3
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         %Christopher Andrew 1
%         candrew_1 = fliplr([1:14]);
%         
%         theCasts = fieldnames(ctdData2);
%         for index = 1:length(theCasts)
%             eval(['hitCTD = ctdData2.' theCasts{index} ';']);
%             theLons(index) = -1*abs(hitCTD(1,3));
%             theLats(index) = hitCTD(1,2);
%             theTimes(index) = hitCTD(1,1);
%             theDepths(index) = max(hitCTD(:,4));
%         end
%         
%         theRun = {'candrew_1'};
%         
%         for indexRun = 1:length(theRun)
%             depthGrid = 1:1:95;
%             eval(['target = ' theRun{indexRun} ';']);
%             clear theDist theGrid
%             for index = 1:length(target)
%                 indexCast = target(index);
%                 eval(['hitCTD = ctdData2.' theCasts{indexCast} ';']);
%                 [junk, theGoods] = unique(hitCTD(:,4));
%                 bads = find(hitCTD(theGoods,4) < 5);
%                 theGoods(bads) = [];
%                 theSalinity = hitCTD(theGoods, 6);
%                 thePressure = hitCTD(theGoods, 7);
%                 theTemperature = hitCTD(theGoods, 8);
%                 theLat = hitCTD(1,2);
%                 theTarg = sw_bfrq(theSalinity,theTemperature,thePressure,theLat);
%                 %theTarg(find(theTarg <= 0)) = NaN;
%                 hitDepths = filter(ones(1,2)/2,1,hitCTD(theGoods, 4));
%                 hitDepths(1) = [];
%                 theInterp = interp1(hitDepths, theTarg, depthGrid);
%                 if indexRun ~=8
%                     lastNonNan = max(find(~isnan(theInterp)));
%                     theInterp(lastNonNan:end) = theInterp(max(find(~isnan(theInterp))));
%                 end
%                 theGrid(:,index) = theInterp(:);
%                 theDist(index) = latlon2dist(pivotLat, pivotLon, theLats(indexCast), theLons(indexCast));
%             end 
%             sectionNum = theRun{indexRun};
%             theDash = strfind(sectionNum, '_');
%             theNum = sectionNum(theDash+1:end);
%         theDiffs = tideTimes-min(theTimes(target));
%         thisEvent = max(find(theDiffs < 0));
%         startAfter = julian2datenum(min(theTimes(target)) - tideTimes(thisEvent));
%         endAfter = julian2datenum(max(theTimes(target)) - tideTimes(thisEvent));
%         
% 
%             theTitle = {[titleStart{indexVar} ', Christopher Andrew Transect ' theNum ': ' ...
%                 datestr(julian2datenum(min(theTimes(target))),1) ' ' ...
%                 datestr(julian2datenum(min(theTimes(target))),15) ' to ' ...
%                 datestr(julian2datenum(max(theTimes(target))),1) ' ' ...
%                 datestr(julian2datenum(max(theTimes(target))),15)]; ...
%                 [datestr(startAfter, 15) ' after low water to ' ...
%                 datestr(endAfter, 15) ' after low water']};
%             dist_km = theDist./1000;
%             figure, orient landscape
%             if strcmp(theExt, '_f')
%                 [c, h] = contourf(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar});
%                 shading flat
%             else
%                 [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits2{indexVar}, 'k');
%             end
%             %clabel(c, h, [min(theLimits{indexVar}):thePick(indexVar):max(theLimits{indexVar})], ...
%             %    'labelspacing', 500)     
%             p = patch(theXs, theYs, 'w');
%             set(p, 'linestyle', 'none')    
%             p2 = patch(theX2s, theY2s, [0.7 0.7 0.7]);
%             set(p2, 'linestyle', 'none');
%             for index = 1:length(target)
%                 l = line([dist_km(index) dist_km(index)], [0 100]);
%                 set(l, 'color', 'k')
%                 t = text(dist_km(index)+0.25, 90, [num2str(target(index), '%02.0f')]);
%                 set(t, 'color', 'k')
%             end  
%             for indexStation = 1:length(stationNames)
%                 st_dist = latlon2dist(pivotLat, pivotLon, ...
%                     stationLats(indexStation), stationLons(indexStation));
%                 st_km = st_dist/1000;
%                 t = text(st_km, -2, stationNames{indexStation});
%                 set(t, 'color', 'k')
%             end    
%             w = text(-0.5, 92, '<----West');
%             set(w, 'horizontala', 'left');
%             e = text(24.5, 92, 'East----->');
%             set(e, 'horizontala', 'right');
%             hold on
%             plot(bathy_km, theBathy, 'k', 'linewidth', 2)
%             axis ij
%             axis([-1 25 0 100])
%             set(gca, 'clim', [min(theLimits{indexVar}) max(theLimits{indexVar})])
%             if strcmp(theExt, '_f')
%                 colorbar
%             end
%             t = title(theTitle);
%             currentPos = get(t, 'position');
%             currentPos(2) = -4;
%             set(t, 'position', currentPos)
%             xlabel('Distance to Ref. Point (Station 27, SW Corner)')
%             ylabel('Depth (m)')
%             print(gcf, '-dpdf', ['candrew_' theNum '_' plotCom{indexVar} theExt])
%             close
%         end
%         
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         %Christopher Andrew 2
%         candrew_2 = fliplr([1:14]);
%         
%         theCasts = fieldnames(ctdData3);
%         for index = 1:length(theCasts)
%             eval(['hitCTD = ctdData3.' theCasts{index} ';']);
%             theLons(index) = -1*abs(hitCTD(1,3));
%             theLats(index) = hitCTD(1,2);
%             theTimes(index) = hitCTD(1,1);
%             theDepths(index) = max(hitCTD(:,4));
%         end
%         
%         theRun = {'candrew_2'};
%         
%         for indexRun = 1:length(theRun)
%             depthGrid = 1:1:95;
%             eval(['target = ' theRun{indexRun} ';']);
%             clear theDist theGrid
%             for index = 1:length(target)
%                 indexCast = target(index);
%                 eval(['hitCTD = ctdData3.' theCasts{indexCast} ';']);
%                 [junk, theGoods] = unique(hitCTD(:,4));
%                 bads = find(hitCTD(theGoods,4) < 5);
%                 theGoods(bads) = [];
%                 theSalinity = hitCTD(theGoods, 6);
%                 thePressure = hitCTD(theGoods, 7);
%                 theTemperature = hitCTD(theGoods, 8);
%                 theLat = hitCTD(1,2);
%                 theTarg = sw_bfrq(theSalinity,theTemperature,thePressure,theLat);
%                 %theTarg(find(theTarg <= 0)) = NaN;
%                 hitDepths = filter(ones(1,2)/2,1,hitCTD(theGoods, 4));
%                 hitDepths(1) = [];
%                 theInterp = interp1(hitDepths, theTarg, depthGrid);
%                 if indexRun ~=8
%                     lastNonNan = max(find(~isnan(theInterp)));
%                     theInterp(lastNonNan:end) = theInterp(max(find(~isnan(theInterp))));
%                 end
%                 theGrid(:,index) = theInterp(:);
%                 theDist(index) = latlon2dist(pivotLat, pivotLon, theLats(indexCast), theLons(indexCast));
%             end
%             sectionNum = theRun{indexRun};
%             theDash = strfind(sectionNum, '_');
%             theNum = sectionNum(theDash+1:end);
% 
%         theDiffs = tideTimes-min(theTimes(target));
%         thisEvent = max(find(theDiffs < 0));
%         startAfter = julian2datenum(min(theTimes(target)) - tideTimes(thisEvent));
%         endAfter = julian2datenum(max(theTimes(target)) - tideTimes(thisEvent));
%         
% 
%             theTitle = {[titleStart{indexVar} ', Christopher Andrew Transect ' theNum ': ' ...
%                 datestr(julian2datenum(min(theTimes(target))),1) ' ' ...
%                 datestr(julian2datenum(min(theTimes(target))),15) ' to ' ...
%                 datestr(julian2datenum(max(theTimes(target))),1) ' ' ...
%                 datestr(julian2datenum(max(theTimes(target))),15)]; ...
%                 [datestr(startAfter, 15) ' after low water to ' ...
%                 datestr(endAfter, 15) ' after low water']};
%             dist_km = theDist./1000;
%             figure, orient landscape
%             if strcmp(theExt, '_f')
%                 [c, h] = contourf(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar});
%                 shading flat
%             else
%                 [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits2{indexVar}, 'k');
%             end
%             %clabel(c, h, [min(theLimits{indexVar}):thePick(indexVar):max(theLimits{indexVar})], ...
%             %    'labelspacing', 500)   
%             p = patch(theXs, theYs, 'w');
%             set(p, 'linestyle', 'none')      
%             p2 = patch(theX2s, theY2s, [0.7 0.7 0.7]);
%             set(p2, 'linestyle', 'none');
%             for index = 1:length(target)
%                 l = line([dist_km(index) dist_km(index)], [0 100]);
%                 set(l, 'color', 'k')
%                 t = text(dist_km(index)+0.25, 90, [num2str(target(index), '%02.0f')]);
%                 set(t, 'color', 'k')
%             end   
%             for indexStation = 1:length(stationNames)
%                 st_dist = latlon2dist(pivotLat, pivotLon, ...
%                     stationLats(indexStation), stationLons(indexStation));
%                 st_km = st_dist/1000;
%                 t = text(st_km, -2, stationNames{indexStation});
%                 set(t, 'color', 'k')
%             end
%             w = text(-0.5, 92, '<----West');
%             set(w, 'horizontala', 'left');
%             e = text(24.5, 92, 'East----->');
%             set(e, 'horizontala', 'right');
%             hold on
%             plot(bathy_km, theBathy, 'k', 'linewidth', 2)
%             axis ij
%             axis([-1 25 0 100])
%             set(gca, 'clim', [min(theLimits{indexVar}) max(theLimits{indexVar})])
%             if strcmp(theExt, '_f')
%                 colorbar
%             end
%             t = title(theTitle);
%             currentPos = get(t, 'position');
%             currentPos(2) = -4;
%             set(t, 'position', currentPos)
%             xlabel('Distance to Ref. Point (Station 27, SW Corner)')
%             ylabel('Depth (m)')
%             print(gcf, '-dpdf', ['candrew_' theNum '_' plotCom{indexVar} theExt])
%             close
%         end
    end
end