% simpleexamples shows a few examples of how these files can be used.

% It is probably good to have both simple and  simple/inputfiles in
% the matlab path.


% First, electrons.

% This should produce the same result as what is shown by the dashed
% lines in Figs. 8 and 9 of Löfgren and Gunell, Phys. Plasmas, vol. 4,
% 3469-3476 (1997).

[vph, UHPpole, lim, omega, f, k, v] = disprelelectron('ele3');

figure(101)
set(gcf,'paperpositionmode','auto')
clf
subplot(2,1,1)
plot(v,f,'k')
set(gca,'fontname','times','fontsize',14)
axis([-2 2 0 1.5])
xlabel('v/v_{d}','fontname','times','fontsize',18)
ylabel('Normalised f(v)','fontname','times','fontsize',18)

subplot(2,1,2)
rootpathplot(vph, UHPpole, lim)
set(gca,'fontname','times','fontsize',14)
axis([-2 2 -0.1 1])
grid on
xlabel('\Re(v_{ph}/v_{d})','fontname','times','fontsize',18)
ylabel('\Im(v_{ph}/v_{d})','fontname','times','fontsize',18)
hold on
plot(vph(1,:),'r') % LW2
plot(vph(2,:),'g') % LW1
plot(vph(7,:),'c') % SSCW
hold off

figure(102)
set(gcf,'paperpositionmode','auto')
clf
subplot(2,3,1)
plot(k,real(omega(2,:)),'g')
set(gca,'fontname','times','fontsize',14)
axis([0 2 -2 0])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Re(\omega)/\omega_{pe}','fontname','times','fontsize',16)
subplot(2,3,4)
plot(k,imag(omega(2,:)),'g')
set(gca,'fontname','times','fontsize',14)
axis([0 2 -0.2 0.5])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Im(\omega)/\omega_{pe}','fontname','times','fontsize',16)

subplot(2,3,2)
plot(k,real(omega(7,:)),'c')
set(gca,'fontname','times','fontsize',14)
axis([0 2 0 2])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Re(\omega)/\omega_{pe}','fontname','times','fontsize',16)
subplot(2,3,5)
plot(k,imag(omega(7,:)),'c')
set(gca,'fontname','times','fontsize',14)
axis([0 2 -0.2 0.5])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Im(\omega)/\omega_{pe}','fontname','times','fontsize',16)

subplot(2,3,3)
plot(k,real(omega(1,:)),'r')
set(gca,'fontname','times','fontsize',14)
axis([0 2 0 2])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Re(\omega)/\omega_{pe}','fontname','times','fontsize',16)
subplot(2,3,6)
plot(k,imag(omega(1,:)),'r')
set(gca,'fontname','times','fontsize',14)
axis([0 2 -0.2 0.5])
xlabel('kv_{d}/\omega_{pe}','fontname','times','fontsize',16)
ylabel('\Im(\omega)/\omega_{pe}','fontname','times','fontsize',16)


% Then ions

% This should be equivalent to the rightmost panels in Fig. 6 of the
% paper by Gunell and Skiff, Phys. Plasmas, vol. 8, 3550-3557 (2001).

[vph, UHPpole, lim, omega, f, k, v] = disprelion('ion2');

figure(201)
set(gcf,'paperpositionmode','auto')
clf
subplot(2,1,1)
plot(v,f,'k')
set(gca,'fontname','times','fontsize',14)
axis([-2 2 0 3])
xlabel('v/c_{s}','fontname','times','fontsize',18)

subplot(2,1,2)
rootpathplot(vph, UHPpole, lim)
set(gca,'fontname','times','fontsize',14)
axis([-2 2 0 1.5])
xlabel('\Re(v_{ph}/c_{s})','fontname','times','fontsize',18)
ylabel('\Im(v_{ph}/c_{s})','fontname','times','fontsize',18)
aa=axis;
hold on;plot(aa(2)*[-1+1i/2/pi 0 1+1i/2/pi],'k--');hold off

% Now, some fluctuation calculations

% This is meant to reproduce the distribution on the last line of Table I
% in the article by Gunell and Skiff, Phys. Plasmas, vol. 9, 2585-2592
% (2002), and the figures generated here should be equivalent to the
% lower panel of Fig. 7 and the dash-dotted line in Fig. 9 of that paper.
run('fluct1')

kplot     = 10.^(-3:0.01:1);
omegaplot = 10.^(-4:0.01:0.5);
k     = 10.^(0.5*(log10(kplot(1:end-1)) + log10(kplot(2:end))));
omega = 10.^(0.5*(log10(omegaplot(1:end-1)) + log10(omegaplot(2:end))));
[kmat omat]=meshgrid(k,omega);

Smat = S1Dofk(kmat, omat, ALLpole, ALLresid, UHPpole, UHPresid, ...
                   n_vec, nb_vec, mass_vec, me);
Smat(Smat<0)=NaN;

figure(301);
set(gcf,'paperpositionmode','auto','renderer','zbuffer')
clf
pp=log10(Smat);
ss=size(pp);
pp=[[pp zeros(ss(1),1)];zeros(1,ss(2)+1)];
surf(kplot,omegaplot,pp)
view(2)
shading flat
set(gca,'xscale','log','yscale','log','fontname','times','fontsize',14)
axis([min(kplot) max(kplot) min(omegaplot) max(omegaplot)])
grid off
hh=colorbar;
set(hh,'fontname','times','fontsize',14)
ylabel(hh,'lg(S/(m_{e}v_{te}^{2}/\omega_{pe}))', ...
       'fontname','times','fontsize',18)
xlabel('k\lambda_{De}','fontname','times','fontsize',18)
ylabel('\omega/\omega_{pe}','fontname','times','fontsize',18)


% integrate S to obtain PEE
omega1 = [0 10.^[-6:0.01:-1]];
Pnn=zeros(size(omega1));
for jj=1:length(omega1)
  oomega=omega1(jj);
  F = @(k) Q1Dofk(k, oomega, ALLpole, ALLresid, UHPpole, UHPresid, ...
                  n_vec, nb_vec, mass_vec, me);
  tmp=zeros(1,6);
  tmp(1)=(1/(2*pi^2))*integral(F,0.001,0.003,'RelTol',1e-9);
  tmp(2)=(1/(2*pi^2))*integral(F,0.003,0.01,'RelTol',1e-9);
  tmp(3)=(1/(2*pi^2))*integral(F,0.01,0.03,'RelTol',1e-9);
  tmp(4)=(1/(2*pi^2))*integral(F,0.03,0.1,'RelTol',1e-9);
  tmp(5)=(1/(2*pi^2))*integral(F,0.1,0.3,'RelTol',1e-9);
  tmp(6)=(1/(2*pi^2))*integral(F,0.3,1,'RelTol',1e-9);
  Pnn(jj)=sum(tmp);
end


figure(302)
set(gcf,'paperpositionmode','auto')
loglog(omega1,Pnn,'k')
set(gca,'fontname','times','fontsize',14)
xlabel('\omega/\omega_{pe}','fontname','times','fontsize',18)
ylabel('P_{nn}^{1D}(\omega)\lambda_{De}^{4}\omega_{pe}', ...
       'fontname','times','fontsize',18)
axis(10.^[-6 -1 -2 1])

