close all 
nCombination = 67; % Number of the combination as indicated by Table S2 (1-80)
fUactH = 0.45; % Activation energy of H in eV in the range 0.20, 0.21, ..., 0.5
nRMSslope = 9; % Surface roughness RMS slope in DEG (9 or 20)
fvYLim_OH = [0 1e17]; % Limits of the y-axis showing the OH column density
fvYLim_H = [0 3e16]; % Limits of the y-axis showing the H column density
nRegion = 11; % Region index (1-18)

% Conversions from column density to ppm
fColDens2PPM_OH = 2.1708e-16;
fColDens2PPM_H = 1.2769e-17;

%% Load data
fileModel = load(sprintf('RMSslope_%02d/Model_Result_Combi_%03d/Model_Uact%.3f_Region_%02d_Combi_%03d.mat',nRMSslope,nCombination,fUactH,nRegion,nCombination));
fileParameter = load(sprintf('RMSslope_%02d/Model_Result_Combi_%03d/Model_Uact%.3f_Combi_%03d_Parameters.mat',nRMSslope,nCombination,fUactH,nCombination));
fileData = load(sprintf('RMSslope_%02d/Region_%02d',nRMSslope,nRegion));

fvMedian = 1e17*fileParameter.fvX(1,end)*fileData.fvMedian;
fmQuantiles = 1e17*fileParameter.fvX(1,end)*fileData.fmQuantiles;
fvLocalTime = fileData.fvLocalTime;
  
%% Display OH and IBD
f = figure;
hold on;
if(any(isnan(fileModel.fvColDens_OH)))
  hText = text(12,mean(fvYLim_OH),'The model did not converge for this set of parameters.','horizontalalignment','center','verticalalignment','middle');
  ylabel('OH column density [m^{-2}]')
  ylabel('OH mass fraction [ppm]')
  xlabel('Local time')
  title(sprintf('Region %d, parameter combination %d',nRegion,nCombination))
  temp = regexp(sprintf('%02d:00,',6:4:18),',','split');
  set(gca,'ylim',fvYLim_OH,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
else
  [hAx,hPlot1,hPlot2] = plotyy(fileModel.fvModelTime,fileModel.fvColDens_OH+1e17*fileParameter.fvX(1,nRegion),fileModel.fvModelTime,fColDens2PPM_OH*(fileModel.fvColDens_OH+1e17*fileParameter.fvX(1,nRegion)));
  set(hPlot1,'linestyle','-','color','b')
  set(hPlot2,'linestyle','-','color','b')
  plot(hAx(1),fvLocalTime,fvMedian,'xr')
  plot(hAx(1),fvLocalTime,fmQuantiles(2,:),'^r')
  plot(hAx(1),fvLocalTime,fmQuantiles(1,:),'vr')
  for nTime = 1:length(fvLocalTime)
      plot(hAx(1),[fvLocalTime(nTime);fvLocalTime(nTime)],fmQuantiles(:,nTime),'-r')
  end
  %hLegend = legend(hAx(1),'Model','IBD','location','northoutside','orientation','horizontal');
  ylabel(hAx(1),'OH column density [m^{-2}]')
  ylabel(hAx(2),'OH mass fraction [ppm]')
  xlabel(hAx(1),'Local time')
  title(sprintf('Region %d, parameter combination %d',nRegion,nCombination))
  temp = regexp(sprintf('%02d:00,',6:4:18),',','split');
  set(hAx(1),'ylim',fvYLim_OH,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
  set(hAx(2),'ylim',fColDens2PPM_OH*fvYLim_OH,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
end

%% Display H
f = figure;
hold on;    
if(any(isnan(fileModel.fvColDens_H)))
  hText = text(12,mean(fvYLim_H),'The model did not converge for this set of parameters.','horizontalalignment','center','verticalalignment','middle');
  ylabel('H column density [m^{-2}]')
  ylabel('H mass fraction [ppm]')
  xlabel('Local time')
  title(sprintf('Region %d, parameter combination %d',nRegion,nCombination))
  temp = regexp(sprintf('%02d:00,',6:4:18),',','split');
  set(gca,'ylim',fvYLim_H,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
else
  [hAx,hPlot1,hPlot2] = plotyy(fileModel.fvModelTime,fileModel.fvColDens_H,fileModel.fvModelTime,fColDens2PPM_H*fileModel.fvColDens_H);
  set(hPlot1,'linestyle','-','color','b')
  set(hPlot2,'linestyle','-','color','b')
  ylabel(hAx(1),'H column density [m^{-2}]')
  ylabel(hAx(2),'H mass fraction [ppm]')
  xlabel(hAx(1),'Local time')
  title(sprintf('Region %d, parameter combination %d',nRegion,nCombination))
  temp = regexp(sprintf('%02d:00,',6:4:18),',','split');
  set(hAx(1),'ylim',fvYLim_H,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
  set(hAx(2),'ylim',fColDens2PPM_H*fvYLim_H,'ytickmode','auto','xlim',[6 18],'xtick',6:4:18,'xticklabel',temp(1:end-1),'ycolor','k');
end