% aa16appendix draws a figure for the appendix in order that we may
% explain the distribution function model. 
%
% HG 2016-09-09

kTe=7*elc;            % Electron temperature
cs=sqrt(kTe/(18*mp)); % This is the ion acoustic speed, which is used
                      % as reference.
aa16case3blue;

% velocity vector normalised to cs
v=[-30000:100:-1100 -1000:10:-110 -100:100 110:10:1000 1100:50000]/1000;

% Total distribution function
ftot=sum(ALLresid.'*ones(size(v))./ ...
       (ones(size(ALLpole)).'*v - ALLpole.'*ones(size(v))));
ftot=real(ftot);

% One for each part
po=[1:4];
f1=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
f1=real(f1);

po=[5:12];
f2=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
f2=real(f2);

po=[13:16];
f3=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
f3=real(f3);

po=[17:22];
f4=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
f4=real(f4);


aa16case1blue;
po=[9:12];
M=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
M=real(M);
po=[13:18];
T=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
T=real(T);
po=[9:18];
f=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
f=real(f);

aa16illustration;
po=[1:4];
M=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
M=real(M);M=M/max(M);
po=[5:10];
T=sum(ALLresid(po).'*ones(size(v))./ ...
       (ones(size(ALLpole(po))).'*v - ALLpole(po).'*ones(size(v))));
T=real(T);T=T/max(T);


ha=[];
figure(78)
clf
set(gcf,'paperpositionmode','auto','position',[150 641 660 330])
ha(1)=axes('position',[0.12 0.17 0.38 0.815]);
ha(2)=axes('position',[0.61 0.17 0.38 0.815]);
semilogy(ha(1),v,f1,'b-.',v,f2,'r--',v,f3,'g--',v,f4,'m-.')
set(ha(1),'nextplot','add')
semilogy(ha(1),v,ftot,'k-','linewidth',2)
set(ha(1),'nextplot','replace')
axis(ha(1),[-4 7 5e-7 11])
set(ha(1),'fontname','times','fontsize',14)
xlabel(ha(1),'v_{\chi}/c_{s}','fontname','times','fontsize',18)
ylabel(ha(1),'f(v_{\chi})','fontname','times','fontsize',18)
hl1=legend(ha(1),'f_{1}','f_{2}','f_{3}','f_{4}','f_{tot}');
set(hl1,'pos',[0.3986 0.46 0.1144 0.5379])
set(ha(1),'ytick',10.^[-10:2:10])

plot(ha(2),v,M,'k--',v,T,'r--',v,M.*T,'b-')
axis(ha(2),[-0.5 6.5 0 1.01])
set(ha(2),'fontname','times','fontsize',14)
xlabel(ha(2),'v_{\chi}/c_{s}','fontname','times','fontsize',18)
hyl=ylabel(ha(2),'f(v_{\chi})','fontname','times','fontsize',18);
set(hyl,'pos',[-1.158 0.4975 1])
hl2=legend(ha(2),'M','T','f');
set(hl2,'pos',[0.89 0.735 0.1068 0.2644])

% add panel labels
for ii=1:length(ha)
  aa=axis(ha(ii));
  if strcmp(get(ha(ii),'yscale'),'linear')
    text(aa(1)+0.05*diff(aa(1:2)), ...
         aa(3)+0.9*diff(aa(3:4)),['(' char(96+ii) ')'], ...
         'fontname','times','fontsize',18,'parent',ha(ii))
  else
    text(aa(1)+0.05*diff(aa(1:2)), ...
         aa(3)*10^(0.9*diff(log10(aa(3:4)))),['(' char(96+ii) ')'], ...
         'fontname','times','fontsize',18,'parent',ha(ii))    
  end
end


if 3==4
  ccc=pwd;
  cd ~
  print -f78 -depsc2 fig-appendix.eps
  cd(ccc)
end



