%  plot_solution_xsection.m
%
% script to make plots from the structure solution_xsection
%
% revised: 14 August 2007  ROK UW-SSEC
% revised: 17 August 2007  ROK UW-SSEC  Include quicklook mode as an option (xhat_update plotted)
% revised: 22 August 2007  ROK UW-SSEC
% revised: 13 Sept   2011  ROK UW-SSEC  Customized for HS3 flight 9/9/2011
% revised: 27 Sept   2011  ROK UW-SSEC  
% revised: 30 Sept   2011  ROK UW-SSEC  
% revised: 03 Oct    2011  ROK UW-SSEC  Added smoothing and DS overlay
% revised: 04 Oct    2011  ROK UW-SSEC  
% revised: 05 Oct    2011  ROK UW-SSEC  
% revised: 20 Oct    2011  ROK UW-SSEC  
% revosed: 08 Nov  2011  ROK UW-SSEC  Automated plotting
% revised: 05 May   2012  ROK UW-SSEC  Fixed crossing UTC date boundary and added Latitude cross-section
% revised: 16 Sept  2012  ROK UW-SSEC  Added Longitude cross-section
% revised: 17 Sept  2012  ROK UW-SSEC  revised plot limits
%
%
% Load instrument Time Latitude Longitude if available
% load 20110909_flight_instrumentTime_Latitude_Longitude.mat
%
%  Plotting time limits
%time_plot_lo = 0;
%time_plot_hi = 24;
%
%chisq_cutoff = 1e10;    % no cutoff
chisq_cutoff = 2e6;    % quicklook cutoff
%chisq_cutoff = 2e5;    % 000 iteration cutoff
%chisq_cutoff =  4e4;   % 000 iteration cutoff
%chisq_cutoff = 9;      % optimal convergence criterion cutoff
%
disp('dir *.mat')
files = dir('solution*.mat');
% %%%%%%
%
for kk = 1:length(files)
    %
    %solution_fileName = input('Enter file name of file containing solution_xsection: ','s')
    solution_filename = files(kk).name;
    disp(['... loading ' solution_filename])
    load(solution_filename)
    %
    plot_mode = 0;   % 0 = xhat_update
    %
    disp('... extract variables and convert units ')
    %  xhat_update is the "next" retrieval 
    xall = solution_xsection.xhat_update;
    %
    dateUTC = solution_xsection.dateUTC;
    timeUTC = solution_xsection.timeUTC;
    dateUTC_string = num2str(dateUTC');
    YR = str2num(dateUTC_string(:,1:2))'+2000;
    MO = str2num(dateUTC_string(:,3:4))';
    DY = str2num(dateUTC_string(:,5:6))';
    dnum = datenum(YR,MO,DY+timeUTC/24);          % compute matlab datenumber
    disp([datestr(dnum(1),'yyyymmddTHHMMSS') ' - ' datestr(dnum(end),'yyyymmddTHHMMSS')])
    %
    if kk == 1
        start_dnum = dnum(1);
        end_dnum = dnum(end);
    else
        end_dnum = dnum(end);
    end
    %
    time_plot_lo = start_dnum;
    time_plot_hi = end_dnum;
    %
    latitude = solution_xsection.Latitude;
    longitude = solution_xsection.Longitude;
    %
%     %  If latitude/longitude is missing then extract it from the instrumentLatitude Longitude
%     for jj = 1: length(timeUTC)
%         ind = min(find(abs(flight_instrumentTime-timeUTC(jj))<0.5/60/60));
%         latitude(jj) = flight_instrumentLatitude(ind);
%         longitude(jj) = flight_instrumentLongitude(ind);
%     end
    %
%     % Find the time region of interest
%     Time_Period_Lo = 9.91/24;
%     Time_Period_Hi = 17.3/24;
%     Time_Period_Lo_dnum = datenum(2011,9,9+Time_Period_Lo);  % 50N
%     Time_Period_Hi_dnum = datenum(2011,9,9+Time_Period_Hi);  % 10N
    Time_Period_Lo_dnum = start_dnum;
    Time_Period_Hi_dnum = end_dnum;
    %
    index_time = find(dnum >= Time_Period_Lo_dnum & dnum <= Time_Period_Hi_dnum);
    %
    %%%%%%%%%%%%%%%%
    %
    % Convert units
    x_dim = solution_xsection.x_dim;
    v_levels = [1:x_dim(1)];
    p = solution_xsection.p;
    xhat_pressure = p(1:x_dim(1));
    xhat_Tskin = xall(end,:);
    xhat_Tdry = xall(1:x_dim(1),:);
    log_vmr_h2o = xall(x_dim(1)+1:x_dim(1)+x_dim(2),:);
    vmr_h2o = exp(log_vmr_h2o)*1.e6;           % ppmv
    Md=28.966;  %Molecular mass of dry air
    Mw=18.016;  %Molecular Mass of water
    xhat_h2omassmix = 0.001*(Mw/Md)*vmr_h2o;   % ppmv to g/kg
    xhat_h2orh = zeros(size(xhat_h2omassmix,1),size(xhat_h2omassmix,2));
    for k = 1:size(xhat_h2omassmix,2)
        pp = p(1:x_dim(1),1); t = xhat_Tdry(:,k); w = xhat_h2omassmix(:,k);
        [rh1,rh2] = mr2rh(pp,t,w);
        xhat_h2orh(:,k) = rh1;
    end
    log_vmr_ozone = xall(x_dim(1)+x_dim(2)+1:x_dim(1)+x_dim(2)+x_dim(3),:);
    xhat_vmr_ozone = exp(log_vmr_ozone)*1.e6;  % ppmv
    %
    %  xa is the mean of the a priori climatology
    xa = solution_xsection.xa;
    xa_pressure = p(1:x_dim(1));
    xa_Tskin = xa(end);
    xa_Tdry = xa(1:x_dim(1));
    log_vmr_h2o = xa(x_dim(1)+1:x_dim(1)+x_dim(2));
    vmr_h2o = exp(log_vmr_h2o)*1.e6;           % ppmv
    Md=28.966;  %Molecular mass of dry air
    Mw=18.016;  %Molecular Mass of water
    xa_h2omassmix = 0.001*(Mw/Md)*vmr_h2o;   % ppmv to g/kg
    pp = p(1:x_dim(1),1); t = xa_Tdry'; w = xa_h2omassmix';
    [rh1,rh2] = mr2rh(pp,t,w);
    xa_h2orh = rh1;
    log_vmr_ozone = xa(x_dim(1)+x_dim(2)+1:x_dim(1)+x_dim(2)+x_dim(3));
    xa_vmr_ozone = exp(log_vmr_ozone)*1.e6;  % ppmv
    %
    % Apply Time Smoother to S-HIS retrievals by layer (Trapezoid)
    %
    Tdry = xhat_Tdry;
    h2orh = xhat_h2orh;
    tmp_Tdry = Tdry;
    tmp_h2orh = h2orh;
    for k = 1:61
        for m = 2:size(Tdry,2)-1
            Tdry(k,m) = (0.5*tmp_Tdry(k,m-1)+tmp_Tdry(k,m)+0.5*tmp_Tdry(k,m+1))/2;
            h2orh(k,m) = (0.5*tmp_h2orh(k,m-1)+tmp_h2orh(k,m)+0.5*tmp_h2orh(k,m+1))/2;
        end
    end
    
    %
    % Assign NaNs when chi sq exceed cutoff and other QC checks
    Tdry(find(h2orh > 98)) = NaN;   % if RH > 98% then set to NaN
    h2orh(find(h2orh > 98)) = NaN;
    NaN_profile = repmat(NaN,61,1);
    reject_ind = find(solution_xsection.chi_sq_optimal > chisq_cutoff);
    for k = 1:length(reject_ind), Tdry(:,reject_ind(k)) = NaN_profile; end
    for k = 1:length(reject_ind), h2orh(:,reject_ind(k)) = NaN_profile; end
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot time-pressure cross-sections 
    figure(1)
    %
    %subplot(211)
    press = p(1:61);
    tmp = Tdry(:,index_time);
    pcolor(dnum(index_time),press,tmp);
    hold on
    shading flat  
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[273-80 273+20]) 
    colorbar
    set(gca,'XLim',[Time_Period_Lo_dnum Time_Period_Hi_dnum])
%     set(gca,'XLim',[time_plot_lo time_plot_hi])
    datetick('x','HH:MM','keeplimits')
    tmp_press = repmat(press,size(Tdry(:,index_time)',1),1);
    tmp_time = repmat(dnum(index_time)',1,61);
    [cs,h] = contour(tmp_time,tmp_press,Tdry(:,index_time)',linspace(273-80,273+20,100/5+1),'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Time (UTC hours) ' datestr(start_dnum,'dd-mmm-yyyy')])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Temperature (K)   ')
    %
    figure(2)
%     %subplot(212)
    tmp = h2orh(:,index_time);
    pcolor(dnum(index_time),p(1:61),tmp);
    hold on
    shading flat; 
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[0 100])
    colorbar
    set(gca,'XLim',[Time_Period_Lo_dnum Time_Period_Hi_dnum])
%     set(gca,'XLim',[time_plot_lo time_plot_hi])
    datetick('x','HH:MM','keeplimits')

    tmp_press = repmat(press,size(h2orh(:,index_time)',1),1);
    tmp_time = repmat(dnum(index_time)',1,61);
    [cs,h] = contour(tmp_time,tmp_press,h2orh(:,index_time)',[10 30 50 70 90],'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Time (UTC hours) ' datestr(start_dnum,'dd-mmm-yyyy')])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Relative Humidity (%) ')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot Latitude-Pressure cross-sections     
    %
    figure(3)
    %
    %subplot(211)
    press = p(1:61);
    tmp = Tdry(:,index_time);
    pcolor(latitude(index_time),press,tmp);
    hold on
    shading flat  
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[273-80 273+20]) 
    colorbar
    %set(gca,'XLim',[10 50])
    tmp_press = repmat(press,size(Tdry(:,index_time)',1),1);
    tmp_latitude = repmat(latitude(index_time)',1,61);
    [cs,h] = contour(tmp_latitude,tmp_press,Tdry(:,index_time)',linspace(273-80,273+20,100/5+1),'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Latitude (deg) ' datestr(start_dnum,'dd-mmm-yyyy') ' '])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Temperature (K)   ')
    %
    figure(4)
%     %subplot(212)
    tmp = h2orh(:,index_time);
    pcolor(latitude(index_time),p(1:61),tmp);
    hold on
    shading flat; 
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[0 100])
    colorbar
    %set(gca,'XLim',[10 50])

    tmp_press = repmat(press,size(h2orh(:,index_time)',1),1);
    tmp_latitude = repmat(latitude(index_time)',1,61);
    [cs,h] = contour(tmp_latitude,tmp_press,h2orh(:,index_time)',[10 30 50 70 90],'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Latitude (deg) ' datestr(start_dnum,'dd-mmm-yyyy') ' '])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Relative Humidity (%)  ')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot Longitude-Pressure cross-sections     
    %
    figure(5)
    %
    %subplot(211)
    press = p(1:61);
    tmp = Tdry(:,index_time);
    pcolor(longitude(index_time),press,tmp);
    hold on
    shading flat  
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[273-80 273+20]) 
    colorbar
    %set(gca,'XLim',[10 50])
    tmp_press = repmat(press,size(Tdry(:,index_time)',1),1);
    tmp_longitude = repmat(longitude(index_time)',1,61);
    [cs,h] = contour(tmp_longitude,tmp_press,Tdry(:,index_time)',linspace(273-80,273+20,100/5+1),'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Longitude (deg) ' datestr(start_dnum,'dd-mmm-yyyy') ' '])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Temperature (K)   ')
    %
    figure(6)
%     %subplot(212)
    tmp = h2orh(:,index_time);
    pcolor(longitude(index_time),p(1:61),tmp);
    hold on
    shading flat; 
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[0 100])
    colorbar
    %set(gca,'XLim',[10 50])

    tmp_press = repmat(press,size(h2orh(:,index_time)',1),1);
    tmp_longitude = repmat(longitude(index_time)',1,61);
    [cs,h] = contour(tmp_longitude,tmp_press,h2orh(:,index_time)',[10 30 50 70 90],'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Longitude (deg) ' datestr(start_dnum,'dd-mmm-yyyy') ' '])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Relative Humidity (%)  ')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

     figure(6)
%     %subplot(212)
    tmp = h2orh(:,index_time);
    pcolor(longitude(index_time),p(1:61),tmp);
    hold on
    shading flat; 
    set(gca,'YLim',[50 1000])
    set(gca,'YScale','log')
    set(gca,'YTick',[50 100 200 300 400 500 600 700 800 900 1000])
    set(gca,'YDir','reverse')
    set(gca,'CLim',[0 100])
    colorbar
    %set(gca,'XLim',[10 50])

    tmp_press = repmat(press,size(h2orh(:,index_time)',1),1);
    tmp_longitude = repmat(longitude(index_time)',1,61);
    [cs,h] = contour(tmp_longitude,tmp_press,h2orh(:,index_time)',[10 30 50 70 90],'k');
%     clabel(cs,h);
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Longitude (deg) ' datestr(start_dnum,'dd-mmm-yyyy') ' '])
    ylabel('Pressure (mb) ')
    title('S-HIS Retrieved Relative Humidity (%)  ')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    figure(7)
%     %subplot(212)
    plot(longitude(index_time),latitude(index_time),'k.');
    hold on
    grid on
    set(gca,'FontSize',12)
    set(gca,'FontWeight','Bold')
    xlabel(['Longitude (deg) ' ])
    ylabel(['Latitude (deg) '])
    title(['S-HIS Flight Track  ' datestr(start_dnum,'dd-mmm-yyyy') ' ' ])
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

    %
end
%
figure(1)
saveas(gcf,'SHIS_Temperature_Time_Pressure_xsection','fig')
saveas(gcf,'SHIS_Temperature_Time_Pressure_xsection','png')
saveas(gcf,'SHIS_Temperature_Time_Pressure_xsection','epsc')
figure(2)
saveas(gcf,'SHIS_RH_Time_Pressure_xsection','fig')
saveas(gcf,'SHIS_RH_Time_Pressure_xsection','png')
saveas(gcf,'SHIS_RH_Time_Pressure_xsection','epsc')
figure(3)
saveas(gcf,'SHIS_Temperature_Latitude_Pressure_xsection','fig')
saveas(gcf,'SHIS_Temperature_Latitude_Pressure_xsection','png')
saveas(gcf,'SHIS_Temperature_Latitude_Pressure_xsection','epsc')
figure(4)
saveas(gcf,'SHIS_RH_Latitude_Pressure_xsection','fig')
saveas(gcf,'SHIS_RH_Latitude_Pressure_xsection','png')
saveas(gcf,'SHIS_RH_Latitude_Pressure_xsection','epsc')
figure(5)
saveas(gcf,'SHIS_Temperature_Longitude_Pressure_xsection','fig')
saveas(gcf,'SHIS_Temperature_Longitude_Pressure_xsection','png')
saveas(gcf,'SHIS_Temperature_Longitude_Pressure_xsection','epsc')
figure(6)
saveas(gcf,'SHIS_RH_Longitude_Pressure_xsection','fig')
saveas(gcf,'SHIS_RH_Longitude_Pressure_xsection','png')
saveas(gcf,'SHIS_RH_Longitude_Pressure_xsection','epsc')
figure(7)
saveas(gcf,'SHIS_Latitude_Longitude','fig')
saveas(gcf,'SHIS_Latitude_Longitude','png')
saveas(gcf,'SHIS_Latitude_Longitude','epsc')
