% gfd_ds_slice_plot.m
% contours interpolated cross-sections of distance and height
% data to be contoured is created in gfd_ds_read_grid
% Original code by Ian Renfrew 1997
% Modified by Stephen Outten 2007, and IAR in Feb 2008
% set flight, leg and then variable using ivar

close all
iprint = 1;             % print to a file

%------------------------------------------
% set up which dataset to use, and which leg
%------------------------------------------
iflight = 268;
if (iflight == 268)
    dir_name =      'F:\data\gfd\FAAM\B268\';
    outdir = 'F:\data\gfd\B268_RTJ\dropsonde_pics\sections\';
    input_leg = 3;                                  % *** input_leg = 1, 2, 3
elseif (iflight == 276)
    dir_name =      'C:\data\gfd\FAAM\B276\';
    outdir =  'C:\data\gfd\B276_BW\dropsonde_pics\';
    input_leg = 14;
end;
outfile1 = ['ds_section_B' num2str(iflight) '_' num2str(input_leg)];

% set variables and load data
if input_leg==1
    load ([dir_name '210207_leg1.mat'])
    title_text1=(' - northern leg');
elseif input_leg==2
    load ([dir_name '210207_leg2.mat'])
    title_text1=(' - along jet');
elseif input_leg==3
    load ([dir_name '210207_leg3.mat'])
    title_text1=(' - southern leg');
elseif input_leg==14
    load ([dir_name '050307_leg1.mat'])
    title_text1=(' - dropsonde leg');
end

%---------------------------------------------------------------
% set variable to be contoured using ivar  *************
% LOOP over ivar is possible
% 1=theta, 2=thetae, 3=temperature, 4=rh, 5 = q, 6 = ws, 7 = wd
% 8 = u, 9 = v, 10 = pressure
% cl = contour levels; clab = contour labels; caxis = color limits
%--------------------------------------------------------------

for ivar = 21:22
    if (ivar == 1)
        zvar = slice_theta;
        cl = 250:2:300; clab = 250:2:300;
        caxis_var = [260 300];
        title_text=['\theta (K)' title_text1];
        outfile = [outfile1  '_theta'];
    elseif (ivar == 2)
        zvar = slice_thetae;
        cl = 250:2:300; clab = 250:2:300;
        caxis_var = [270 310];
        title_text=['\theta_{e}  (K)' title_text1];
        outfile = [outfile1  '_thetae'];
    elseif (ivar == 3)
        zvar = slice_t;
        cl = -30:2:10; clab = -30:4:10;
        caxis_var = [-30 10];
        title_text=['T (^{o}C)' title_text1];
        outfile = [outfile1  '_t'];
    elseif (ivar == 4)
        zvar = slice_rh;
        cl = 0:10:100; clab = 0:20:100;
        caxis_var = [0 100];
        title_text=['RH (%)' title_text1];
        outfile = [outfile1  '_rh'];
    elseif (ivar == 5)
        zvar = slice_q;
        cl = 0:0.5:5; clab = 0:1:5;
        caxis_var = [0 5];
        title_text=['q (g kg^{-1})' title_text1];
        outfile = [outfile1  '_q'];
    elseif (ivar == 6)
        zvar = slice_ws;
        cl = 0:2:50; clab = 16:4:50;
        caxis_var = [0 50];
        title_text=['ws (m s^{-1})' title_text1];
        outfile = [outfile1  '_ws'];
    elseif (ivar == 7)
        zvar = slice_wd;
        cl = 0:10:360; clab = 20:20:360;
        caxis_var = [0 90];
        title_text=['wd (deg.)' title_text1];
        outfile = [outfile1  '_wd'];
    elseif (ivar == 8)
        zvar = slice_u;
        cl = -40:4:40; clab = -30:4:30;
        caxis_var = [-30 30];
        title_text=['u (m s^{-1})' title_text1];
        outfile = [outfile1  '_u'];
    elseif (ivar == 9)
        zvar = slice_v;
        cl = -50:4:50; clab = -50:4:50;
        caxis_var = [-40 40];
        title_text=['v (m s^{-1})' title_text1];
        outfile = [outfile1  '_v'];
    elseif (ivar == 10)
        zvar = slice_p;
        cl = 500:10:1100; clab = 500:50:1100;
        caxis_var = [500 1100];
        title_text=['pressure (hPa)' title_text1];
        outfile = [outfile1  '_p'];
    elseif (ivar == 21)
        zvar = slice_vg2;
        cl = -50:4:50; clab = -50:4:50;
        caxis_var = [-40 40];
        title_text=['v_g (m s^{-1})' title_text1];
        outfile = [outfile1  '_vg'];
     elseif (ivar == 22)
        zvar = slice_v2;
        cl = -50:4:50; clab = -50:4:50;
        caxis_var = [-40 40];
        title_text=['v_{s} (m s^{-1})' title_text1];
        outfile = [outfile1  '_vs'];
    end;
    
    %----------------------------------------------------------------------
    % do contouring and plot annotations
    %----------------------------------------------------------------------
    figure(ivar)
    if (ivar < 20)
        [c_cont,h_cont] = contourf(distance, zi, zvar, cl, 'k');
    elseif (ivar > 20)
        [c_cont,h_cont] = contourf(distance2, zi, zvar, cl, 'k');
    end;
    hold on
    caxis([caxis_var(1) caxis_var(2)]);
    clabel(c_cont,h_cont,clab, 'fontsize', 10,'label',250); 
    set(gca, 'fontsize', 12)
    set(gca, 'XTick', 0:20:300)
    set(gca, 'YTick', 0:1000:6000)
    ylabel('height (m)')
    xlabel('distance (km)')
    set(gca,'PlotBoxAspectRatio', [max(distance)/70 z_max/4000 1]); % makes consistent aspect ratio
    pos = get(gca, 'Position');
    pos = [pos(1) pos(2) pos(3)*(max(distance)/250) pos(4)];   % use for a consistent width
    set(gca, 'Position', pos);
    orient portrait
    title(title_text)
    
    %------------------------------------------------------
    % Alternative, use plot to add markers at ds locations
    %------------------------------------------------------
    if input_leg<3
        ds_tagx = max(leg1_distance) - leg1_distance;   %flipped leg
    else
        ds_tagx = leg1_distance;
    end;
    if max(leg1_distance) > max(distance)     % substitute end ds for edge of plot
        ds_tagx(end) = max(distance);
    end;
    ds_tagy = z_max*ones(size(leg1_distance));
    ds_tag = plot(ds_tagx, ds_tagy,'v','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor','k');
    
    %-------------------------------------------------
    % print figure to file
    %-------------------------------------------------
    if (iprint == 1)
        print('-dpng', [outdir outfile '.png']);
        %       print('-depsc', [outdir outfile '.eps']);
        disp(['printing ' outdir outfile]);
    end
    clear zvar cl clab caxis_var title_text outfile
    
end; % of ivar loop