function plot_2spt_function_pdf(spt_names, fnameout)    

%% read and plot in 2 spt files, then save to 1 page pdf
% example input:
% spt_names{1} = 'zci20151026sgeeac.101';
% spt_names{2} = 'C:\Users\Jacob_2\Documents\Research\EM27\examine_spt_files\fits\jaxa_20160622_fits\zjm20160622s0e00c_co4233.0003';
% fnameout = 'spt_out.pdf';
% figure(1); clf;

%% data input check
if (~isequal(size(spt_names), [2,1]) && ~isequal(size(spt_names), [1,2])) || ~iscell(spt_names);
  error('spt_name should be 2x1 cell')
end
% spt_name = flipud(spt_name);

if ~isempty(spt_names{2});
  n_spts = 2; %number of spt files
else %only one plot
  n_spts = 1;
end
  
hsp = NaN*zeros(n_spts*2,1); %pre-allocate subplot handles
hl = NaN*zeros(n_spts,1); %pre-allocate legend handles
for i2=1:n_spts;
  %% read in file, determine how many species to plot
  spt = read_spt_fast(spt_names{i2}); %read in .spt file

  fn = fieldnames(spt);
  nfg = length(fn)-12; % Number of fitted gases
  col = lines(nfg);
  if isequal(fn{nfg+3},'solar');
    lgd_lst = fn(2:end-9);
  else %no solar gas fitting
    nfg = nfg-1;
    lgd_lst = fn(2:end-10);
  end
  maingas = fn{4};
  ctr_wn = floor(median(spt.Freq)); %center wavenumber
  
  %% plot residuals
  hold on;
  hsp(i2*2-1) = subplot(n_spts*2,1,i2*2-1);

  hold on;
  resid_pct = (spt.Tm-spt.Tc)*100; %calculate the residuals in percent
  plot(spt.Freq,resid_pct,'k','linewidth',1.5)
  grid on; box on
  ylabel('%','fontsize',12,'fontweight','bold')
  title([spt.filename ',  ' num2str(spt.zenith) '^\circ SZA,  ' num2str(spt.rms) '% RMS,  ' num2str(spt.tot_cols(1)) ' molec\cdotcm^{-2}'],'fontsize',12,'fontweight','bold')
  set(hsp(i2*2-1),'XTickLabel','','xlim',[spt.Freq_min spt.Freq_max],'fontsize',12,'fontweight','bold');

  %% find the largest residual, for "smart" assigning of ylimits on residuals
  resid_lim = abs(resid_pct);
  if length(resid_lim)>100;
    resid_lim = sort(resid_lim); 
    resid_lim = resid_lim(round(0.995*length(resid_lim))); %takes a given percentile as maximum
  else
    resid_lim = max(resid_lim);
  end

  if resid_lim>1;
    resid_lim=ceil(resid_lim);
  else %less than 1
    resid_lim=ceil(resid_lim*10)/10;
  end

  ylim([-1 1]*resid_lim)
  text(spt.Freq_max+(spt.Freq_max-spt.Freq_min)*0.02, 0, {maingas, ctr_wn},'fontsize',12','fontweight','bold')


  %% plot spectral window and species
  hsp(i2*2) = subplot(n_spts*2,1,i2*2);
  hold on;

  plot(spt.Freq,spt.Tm,'k.', 'Markersize', 12) %plot measured points
  plot(spt.Freq,spt.Tc,'k', 'linewidth', 1.5) %plot calculated points
  plot(spt.Freq,spt.(fn{1+3}),'-','color',col(1,:),'linewidth',2)
  for i1=2:nfg
      plot(spt.Freq,spt.(fn{i1+3}),'-','color',col(i1,:),'linewidth',1.5)
  end

  y_min = min([spt.Tm; spt.Tc]);
  y_min = floor(y_min*20)/20; %floor to nearest 0.05
  % y_max = max([spt.Tm; spt.Tc]);
  % if y_max>1; y_max=1.01; else y_max =1;end %set upper limit
  y_max = 1.02;
  ylim([y_min, y_max])
  ylabel('Transmittance')

  set(hsp(i2*2),'xlim',[spt.Freq_min spt.Freq_max],'fontsize',12,'fontweight','bold')
  hl(i2)=legend(lgd_lst,3);

  set(hl(i2),'Box','off','Location','eastOutside','Interpreter','none')
  % set(hl,'Interpreter','none')
  % pos = get(gca,'Position');
  grid on; box on;
  xlabel('Wavenumber (cm^{-1})')
  linkaxes(hsp((i2*2-1):i2*2),'x') %make xlimits the same for both
  
  
  set(hsp(i2*2-1), 'Pos', [0.10, 0.33 + (n_spts-i2)*0.5, 0.8, 0.10])
  set(hsp(i2*2), 'Pos', [0.10, 0.08 + (n_spts-i2)*0.5, 0.8, 0.25])
  
  sp_xlims = get(hsp(i2*2-1), 'Position');
  lgd_pos = get(hl(i2), 'Position'); %legend handle position
  set(hl(i2), 'Pos', [sp_xlims(1)+sp_xlims(3)+0.005, lgd_pos(2:4)]); %move legend closer to axes
end

%% position subplots
% set(hsp(1), 'Pos', [0.10, 0.73, 0.8, 0.23])
% set(hsp(2), 'Pos', [0.10, 0.14, 0.8, 0.58])



%% save figure
figh = gcf; %get current figure handle
figh.Units = 'inches';
figh.PaperUnits = 'inches';
figh.OuterPosition = [1 1 8 10];
figh.PaperPosition = [0 0 8 10];
figh.PaperPositionMode = 'manual';


print('-dpdf', fnameout, '-r300')


