function ctd_bvfq_profile
%
%   ctd_bvfq_profile
%
%   Function to create profile 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.

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

targDepths = [0:4:100];

%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)));

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]);

colList = {'b'; 'r'; 'g'; 'c'; 'y'; 'm'; 'k'; ...
    'b:'; 'r:'; 'g:'; 'c:'; 'y:'; 'm:'; 'k:'; ...
    'b-.'; 'r-.'; 'g-.'; 'c-.'; 'y-.'; 'm-.'; 'k-.'};

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 = 1:length(theRun)
        figure, orient tall, hold on
        eval(['target = ' theRun{indexRun} ';']);
        clear castList
        for index = 1:length(target)
            clear newTarg
            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) = [];
            for indexDepth = 1:length(targDepths)
                depthSpace = median(diff(targDepths));
                inDepth = ...
                    find((hitDepths > (targDepths(indexDepth)-depthSpace/2)) & ...
                    (hitDepths <= (targDepths(indexDepth)+depthSpace/2)));
                if ~isempty(inDepth)
                    newTarg(indexDepth) = nanmean(theTarg(inDepth));
                else
                    newTarg(indexDepth) = NaN;
                end
            end
            plot(newTarg(:), targDepths(:), colList{index})
            castList{index} = [theCasts{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']};

        axis ij
        ylabel('Depth (m)')
        xlabel('BV Frequency')
        ylim([0 100])
        xlim([-2e-3 8e-3])
        box on
        l = legend(castList, 'location', 'northeastoutside');

        t = title(theTitle);
        currentPos = get(t, 'position');
        currentPos(2) = -4;
        set(t, 'position', currentPos)
        box('on')
        print(gcf, '-dpdf', ['argo_' theNum '_' plotCom{indexVar} '_p'])
        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)
            figure, hold on, orient tall
            eval(['target = ' theRun{indexRun} ';']);
            clear castList
            for index = 1:length(target)
                clear newTarg
                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) = [];
                for indexDepth = 1:length(targDepths)
                    depthSpace = median(diff(targDepths));
                    inDepth = ...
                        find((hitDepths > (targDepths(indexDepth)-depthSpace/2)) & ...
                        (hitDepths <= (targDepths(indexDepth)+depthSpace/2)));
                    if ~isempty(inDepth)
                        newTarg(indexDepth) = nanmean(theTarg(inDepth));
                    else
                        newTarg(indexDepth) = NaN;
                    end
                end
                plot(newTarg(:), targDepths(:), colList{index})
                castList{index} = [theCasts{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']};
            
        l = legend(castList, 'location', 'northeastoutside');
                
                axis ij
                ylabel('Depth (m)')
                xlabel('BV Frequency')
                ylim([0 100])
                xlim([-2e-3 8e-3])
                box on

                t = title(theTitle);
                currentPos = get(t, 'position');
                currentPos(2) = -4;
                set(t, 'position', currentPos)
                print(gcf, '-dpdf', ['candrew_' theNum '_' plotCom{indexVar} '_p'])
                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)
            figure, hold on, orient tall
            eval(['target = ' theRun{indexRun} ';']);
            clear castList
            for index = 1:length(target)
                clear newTarg
                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) = [];
                for indexDepth = 1:length(targDepths)
                    depthSpace = median(diff(targDepths));
                    inDepth = ...
                        find((hitDepths > (targDepths(indexDepth)-depthSpace/2)) & ...
                        (hitDepths <= (targDepths(indexDepth)+depthSpace/2)));
                    if ~isempty(inDepth)
                        newTarg(indexDepth) = nanmean(theTarg(inDepth));
                    else
                        newTarg(indexDepth) = NaN;
                    end
                end
                plot(newTarg(:), targDepths(:), colList{index})   
                castList{index} = [theCasts{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']};
            
                axis ij
                ylabel('Depth (m)')
                xlabel('BV Frequency')
                ylim([0 100])
                xlim([-2e-3 8e-3])
                box on
        l = legend(castList, 'location', 'northeastoutside');
                
            t = title(theTitle);
            currentPos = get(t, 'position');
            currentPos(2) = -4;
            set(t, 'position', currentPos)
            print(gcf, '-dpdf', ['candrew_' theNum '_' plotCom{indexVar} '_p'])
            close
        end
    end
end