function ctd_plots
%
%   ctd_plots
%
%   Function to create plots of MBIWE CTD data contoured (filled or not)
%   with the bathymetry overlaid.

%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 = {'Temperature'; 'Salinity'; 'Attenuation'; 'Sigma Theta'};
rowIndex= [8 6 9 5];
theLimits = {[4:0.2:16]; [30.5:0.05:32]; [0.5:0.05:2]; [20:0.05:25.5]};
plotCom = {'temp'; 'sal'; 'attn'; 'sigma_theta'};
thePick = {[5:5:15]; [30.5:0.5:32]; [0.5:0.2:2]; [20:0.5:25.5]};

addpath('C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\MFILES')

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:3,:) = [];
tideTimes = datenum2julian(datenum(presText(:,1))); 

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

depthGrid = 1: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 = 4%1:4
    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)
        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));    
            theTarg = hitCTD(theGoods, rowIndex(indexVar));
            theTarg(find(theTarg <= 0)) = NaN;
            theInterp = interp1(hitCTD(theGoods,4), 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});
            else
                [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar}, 'k');
            end 
            clabel(c, h, thePick{indexVar}, ...
                'labelspacing', 288)
            shading flat    
            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});
            else
                [c,h] = contour(theTimes(section_8)', depthGrid(:), theGrid, theLimits{indexVar}, ...
                    'k');
            end 
            clabel(c, h, thePick{indexVar}, ...
                'labelspacing', 288)
            shading flat
            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));    
                theTarg = hitCTD(theGoods, rowIndex(indexVar));
                theTarg(find(theTarg <= 0)) = NaN;
                theInterp = interp1(hitCTD(theGoods,4), theTarg, depthGrid);
                lastNonNan = max(find(~isnan(theInterp)));
                theInterp(lastNonNan:end) = theInterp(max(find(~isnan(theInterp))));
                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});
            else
                [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar}, 'k');
            end
            clabel(c, h, thePick{indexVar}, ...
                'labelspacing', 288)  
            shading flat    
            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));    
                theTarg = hitCTD(theGoods, rowIndex(indexVar));
                theTarg(find(theTarg <= 0)) = NaN;
                theInterp = interp1(hitCTD(theGoods,4), theTarg, depthGrid);
                lastNonNan = max(find(~isnan(theInterp)));
                theInterp(lastNonNan:end) = theInterp(max(find(~isnan(theInterp))));
                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});
            else
                [c,h] = contour(dist_km(:)', depthGrid(:), theGrid, theLimits{indexVar}, 'k');
            end
            clabel(c, h, thePick{indexVar}, ...
                'labelspacing', 288)
            shading flat    
            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