% This script generates all figures from the paper and aves them 
% as png files in the "figures" folder.

close all;
clear all;

save_images=1;
save_dir='../figures/';

unit_str_BR='Gauss';
unit_str_FLUX='Mx';
unit_ste_EUV='log10(I)';
Rsun=6.957e10;
Rsun2=(Rsun^2);
r1au=1.496e13;
gauss2nt=1e5;
screen_width    = 2560;
screen_height   = 1352;
screen_x_offset = 1;
screen_y_offset = 89;
fontsz=25;
lw=2;
load('psi_pals.mat');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  FIGURE 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Observed synchronic EUV
euv_file_name='../data/EUV/CR2098_S_20100708_OBS_EUV.hdf';
[MAP,tvec,pvec]=rdhdf_2d(euv_file_name);
MAP=double(MAP);
tvec=double(tvec);
pvec=double(pvec);
sinlat_vec=cos(tvec);

figure('Position',[507 498 1456 900],'Color','w');
h_map=pcolor(pvec,sinlat_vec,MAP);
set(gca,'FontSize',fontsz,'FontName','Arial');
shading interp;
cmaptmp=psi_pal_eit195;
cmaptmp(1,:)=[0 0 0.5];
colormap(cmaptmp);
cbh=colorbar;
caxis([[-0.000001 3]]);
xlabel('\phi','FontSize',fontsz,'FontName','Arial');
cblab_h=cblabel(unit_ste_EUV,'FontSize',fontsz,'FontName','Arial');
axis(gca,[0 2*pi -1 1]);
set(gca,'YTick',[-1 -0.5 0 0.5 1],'YTickLabel',{'-1','-1/2','0','1/2','1'}); 
ylabel('Sin-Latitude','FontSize',fontsz,'FontName','Arial');
grid on;
set(gca,'layer','top');
set(gca,'XTick',[0 pi/2 pi 3*pi/2 2*pi],'XTickLabel',{'0','90','180','270','360'},'FontName','symbol');
set(cbh,'FontSize',fontsz,'FontName','Arial');
set(cbh,'YTick',[0 1 2 3],'YTickLabel',{'1','2','3','4'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure1a_observed_synchronic_euv.png'],'-r96');
end
%
% Observed synchronic CHM
%
chm_file_name='../data/CHM/CR2098_S_20100708_OBS_CHM.hdf';
[CHM,~,~]=rdhdf_2d(chm_file_name);
CHM=double(CHM);
hold on;
[~,h_cont]=contour(pvec,sinlat_vec,CHM>0.5,[0 0],'m','LineWidth',lw);
hold off;
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure1a_observed_synchronic_euv_w_chm.png'],'-r96');
end
delete(h_cont);
CHM(MAP<-99)=1.01;
set(h_map,'CData',CHM);
caxis([0 1.01]);
cmaptmp=flipud(gray);
cmaptmp(end,:)=[0 0 0.5];
colormap(cmaptmp);
set(cblab_h,'String',' ');
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure1b_observed_synchronic_chm.png'],'-r96');
end
%
% ADAPT Model TMHD OPENFLUX
%
file_name='../data/OFM/CR2098_S_20100707_ADAPT_FAR_MHD_OFM.hdf';
[CHM,~,~]=rdhdf_2d(file_name);
CHM=double(CHM);
set(h_map,'CData',CHM);
colormap(flipud(gray));
caxis([0 1]);
set(cblab_h,'String',' ');
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});
hold on;
[~,h_cont]=contour(gca,pvec,sinlat_vec,MAP<-99,[1 1],'LineColor',[0 1 1],'LineWidth',lw+1,'LineStyle','--');
hold off;
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure1c_model_tmhd_adapt_far_open_flux.png'],'-r96');
end
delete(h_cont)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  FIGURE 3 (figure 2 below)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Observed Synoptic CHM
%
file_name='../data/CHM/CR2098_D_OBS_CHM.hdf';
[CHM,~,~]=rdhdf_2d(file_name);
CHM=double(CHM);
set(h_map,'CData',CHM);
colormap(flipud(gray));
caxis([0 1]);
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure3b_observed_synoptic_chm.png'],'-r96');
end
%
% Observed Synoptic EUV
%
file_name='../data/EUV/CR2098_D_OBS_EUV.hdf';
[MAP,~,~]=rdhdf_2d(file_name);
MAP=double(MAP);
set(h_map,'CData',MAP);
colormap(psi_pal_eit195)
caxis([0 3]);
set(cblab_h,'String',unit_ste_EUV);
set(cbh,'YTick',[0 1 2 3],'YTickLabel',{'1','2','3','4'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure3a_observed_synoptic_euv.png'],'-r96');
end
%
% Observed Synoptic Magnetogram
%
file_name='../data/BR0/CR2098_D_VSM_BR0.hdf';
[MAP,~,~]=rdhdf_2d(file_name);
MAP=double(MAP);
set(h_map,'CData',MAP);
colormap(psi_pal_bluered)
caxis([-20 20]);
set(cblab_h,'String',unit_str_BR);
set(cbh,'YTick',[-20 -10 0 10 20],'YTickLabel',[-20 -10 0 10 20]);
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure3c_observed_synoptic_magnetogram_shine.png'],'-r96');
end
%
% Observed Synoptic Magnetogram (SHINE) masked with Obs Synoptic CHM
%
MAP=CHM.*MAP;
set(h_map,'CData',MAP);
hold on;
[~,h_cont]=contour(pvec,sinlat_vec,CHM>0.5,[0 0],'k','LineWidth',1);
hold off;
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure3d_obs_syn_mag_shine_masked_w_obs_syn_chm.png'],'-r96');
end
delete(h_cont);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  FIGURE 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Model Synthetic Synoptic EUV from TMHD
%
file_name='../data/EUV/CR2098_D_VSM_PFIT_MHD_EUVSYNTH.hdf';
[MAP,~,~]=rdhdf_2d(file_name);
MAP=double(MAP);
set(h_map,'CData',MAP);
colormap(psi_pal_eit195)
caxis([0 3]);
set(cblab_h,'String',unit_ste_EUV);
set(cbh,'YTick',[0 1 2 3],'YTickLabel',{'1','2','3','4'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure4a_model_tmhd_synthetic_euv.png'],'-r96');
end
%
% Model Synthetic Synoptic CHD from TMHD
%
file_name='../data/CHM/CR2098_D_VSM_PFIT_MHD_CHMSYNTH.hdf';
[CHM,~,~]=rdhdf_2d(file_name);
CHM=double(CHM);
set(h_map,'CData',CHM);
colormap(flipud(gray));
caxis([0 1]);
set(cblab_h,'String',' ');
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure4b_model_tmhd_synthetic_chd_minsymdiff.png'],'-r96');
end
%
% Model Magnetogram masked with synthetic EUV CHD
%
file_name='../data/BR0/CR2098_D_VSM_PFIT_BR0.hdf';
[MAP,~,~]=rdhdf_2d(file_name);
MAP=double(MAP);
MAP=MAP.*CHM;
set(h_map,'CData',MAP);
colormap(psi_pal_bluered)
caxis([-20 20]);
set(cblab_h,'String',unit_str_BR);
set(cbh,'YTick',[-20 -10 0 10 20],'YTickLabel',[-20 -10 0 10 20]);
hold on;
[~,h_cont]=contour(pvec,sinlat_vec,CHM>0.5,[0 0],'k','LineWidth',1);
hold off;
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure4c_model_tmhd_synthetic_chd_masking_br_from_sim.png'],'-r96');
end
delete(h_cont);
%
% Model Open Flux from TMHD
%
file_name='../data/OFM/CR2098_D_VSM_PFIT_MHD_OFM.hdf';
[CHM,~,~]=rdhdf_2d(file_name);
CHM=double(CHM);
set(h_map,'CData',CHM);
colormap(flipud(gray));
caxis([0 1]);
set(cblab_h,'String',' ');
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure4d_model_vsm_tmhd_open_flux.png'],'-r96');
end
%
% Model Magnetogram masked with model openflux
%
file_name='../data/BR0/CR2098_D_VSM_PFIT_BR0.hdf';
[MAP,~,~]=rdhdf_2d(file_name);
MAP=double(MAP);
MAP=MAP.*CHM;
set(h_map,'CData',MAP);
colormap(psi_pal_bluered)
caxis([-20 20]);
set(cblab_h,'String',unit_str_BR);
set(cbh,'YTick',[-20 -10 0 10 20],'YTickLabel',[-20 -10 0 10 20]);
hold on;
[~,h_cont]=contour(pvec,sinlat_vec,CHM>0.5,[0 0],'k','LineWidth',1);
hold off;
drawnow;
if(save_images==1)
  screen2png(gcf,[save_dir 'figure4e_model_vsm_tmhd_openflux_masking_br_from_sim.png'],'-r96');
end
delete(h_cont);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  FIGURE 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%S:
%
% VSM OFM of PFSS for various Rss values.
%
colormap(flipud(gray));
caxis([0 1]);
set(cblab_h,'String',' ');
set(cbh,'YTick',[0 1],'YTickLabel',{'No CH','CH'});

rss_list={'25','20','14','13'};
fig_list={'a','b','c','d'};
n_rss=4;

for j=1:n_rss
  file_name=['../data/OFM/CR2098_D_VSM_PFSS' rss_list{j} '_OFM.hdf'];
  [CHM,~,~]=rdhdf_2d(file_name);
  CHM=double(CHM);
  set(h_map,'CData',CHM);
  drawnow;
  if(save_images==1)
    screen2png(gcf,[save_dir 'figure2' fig_list{j} '_model_vsm_pfss_openflux_multi_rss.png'],'-r96');
  end
end

