% aa16figs7to9 is a program intended to produce figures showing
% dispersion relations for use in a paper in Astronomy and Astrophysics
% in 2016. 
%
% This file is used to compute updated dispersion relations that fit
% Langmuir probe results better than the first guesses did.
%
% HG 2016-09-09

mass=18*mp;        % assume water ions are 18 times heavier than protons

% Case 1: 08:00 on 20 January 2015
hb=[];
figure(800)
clf
set(gcf,'paperpositionmode','auto','renderer','zbuffer')
set(gcf,'position',[5 253 560 440])
hb(1)=axes('position',[0.14 0.565 0.83 0.43]);
hb(2)=axes('position',[0.14 0.13 0.83 0.43]);

% red curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
ne=6.963607702e+05;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));
if ~exist('matfiles/aa16case1redwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run it first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1redwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  vphinit(2)  = 2*UHPpole(2)-vphinit(23);

  Aa=vphinit(3);Ab=vphinit(4);
  vphinit(3)  = 2*UHPpole(3)-vphinit(24);
  vphinit(4)  = Aa;
  vphinit(25)  = 2*UHPpole(4)-vphinit(4);

  Ac=vphinit(6);Ad=vphinit(9);Ae=vphinit(10);Af=vphinit(11);
  Ag=vphinit(30);Ah=vphinit(31);Ai=vphinit(32);
  vphinit(26) = vphinit(27);
  vphinit(6)  = vphinit(29);
  vphinit(27)  = 2*UHPpole(6)-vphinit(6);
  vphinit(28)  = 2*UHPpole(7)-vphinit(7);
  vphinit(29)  = 2*UHPpole(8)-vphinit(8);
  vphinit(9)   = UHPpole(9)-conj(vphinit(7)-UHPpole(7));
  vphinit(30)  = 2*UHPpole(9)-vphinit(9);
  vphinit(10)  = 0.01875766132741+0.00384778498889i;
  vphinit(31)  = 2*UHPpole(10)-vphinit(10);
  vphinit(11)  = vphinit(35);
  vphinit(32)  = 2*UHPpole(11)-vphinit(11);

  Aj=vphinit(12);Ak=vphinit(13);Al=vphinit(34);
  vphinit(12)  = Ai;
  vphinit(13)  = vphinit(37);
  vphinit(34)  = 2*UHPpole(13)-vphinit(13);
  vphinit(14)  = UHPpole(14)-conj(vphinit(12)-UHPpole(12));
  vphinit(35)  = 2*UHPpole(14)-vphinit(14);

  vphinit(36)  = 2*UHPpole(15)-vphinit(15);
  vphinit(37)  = 2*UHPpole(16)-vphinit(16);

  vphinit(20)  = vphinit(41);
  vphinit(39)  = 2*UHPpole(18)-vphinit(18);
  vphinit(40)  = 2*UHPpole(19)-vphinit(19);
  vphinit(41)  = 2*UHPpole(20)-vphinit(20);
  
  % Here we compute the dispersion relation. It may take a while.
  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1redwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case1redwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case1redwe.mat
end

set(hb(1:2),'nextplot','add')
plot(hb(1),k/lambdaDe,real(omega(2,:))*fpe,'r')
plot(hb(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'r')
set(hb,'box','on')


% blue curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
ne=1.652659692e+08;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;

lambdaDe=sqrt(eps0*kTe/(ne*elc^2));


if ~exist('matfiles/aa16case1bluewe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  % now with electrons
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run the original distribution first for the highest k value to get an
  % initial guess. 
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1bluewe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(2);
  vphinit(2)  = vphinit(4);

  vphinit(3)  = UHPpole(3)+(vphinit(1)-UHPpole(1));
  vphinit(24) = UHPpole(3)+(vphinit(22)-UHPpole(1));
  vphinit(4)  = UHPpole(4)+(vphinit(2)-UHPpole(2));
  vphinit(25) = UHPpole(4)+(vphinit(23)-UHPpole(2));
  
  Ab=vphinit(6);
  vphinit(26)  = 2*UHPpole(5)-vphinit(5);
  vphinit(6)  = vphinit(29);
  vphinit(27)  = 2*UHPpole(6)-vphinit(6);

  Ac=vphinit(8);Ad=vphinit(29);Ae=vphinit(30);
  vphinit(7)   = UHPpole(7) + (vphinit(5)-UHPpole(5));
  vphinit(28)  = 2*UHPpole(7)-vphinit(7);
  vphinit(29)  = 2*UHPpole(8)-vphinit(8);
  vphinit(9)   = UHPpole(9) + (vphinit(6)-UHPpole(6));
  vphinit(30)  = 2*UHPpole(9)-vphinit(9);

  Ag=vphinit(11);Ah=vphinit(32);
  vphinit(11)  = vphinit(12);
  vphinit(32)  = vphinit(34);
  vphinit(10)  = UHPpole(10)-conj(vphinit(11)-UHPpole(11));
  vphinit(31)  = 2*UHPpole(10)-vphinit(10);

  Ai=vphinit(12);Aj=vphinit(13);Ak=vphinit(34);
  vphinit(12)  = Ah;
  vphinit(33)  = 2*UHPpole(12)-vphinit(12);
  vphinit(13)  = Ag;
  vphinit(34)  = 2*UHPpole(13)-vphinit(13);
  vphinit(14)  = UHPpole(14)-conj(vphinit(12)-UHPpole(12));
  vphinit(35)  = 2*UHPpole(14)-vphinit(14);

  vphinit(16)  = Aj;
  vphinit(37)  = 2*UHPpole(16)-vphinit(16);
  vphinit(15)  = UHPpole(15)-conj(vphinit(16)-UHPpole(16));
  vphinit(36)  = 2*UHPpole(15)-vphinit(15);

  vphinit(20)  = vphinit(41);
  vphinit(39)  = 2*UHPpole(18)-vphinit(18);
  vphinit(40)  = 2*UHPpole(19)-vphinit(19);
  vphinit(41)  = 2*UHPpole(20)-vphinit(20);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1bluewe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case1bluewe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case1bluewe.mat
end

plot(hb(1),k/lambdaDe,real(omega(2,:))*fpe,'b')

plot(hb(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'b')


% green curve
kTe=10.2*elc;      % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
clear nbeam nbulk ntot; ne=8.31845e7;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;

lambdaDe=sqrt(eps0*kTe/(ne*elc^2));


if ~exist('matfiles/aa16case1greenwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  % now with electrons
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run the original distribution first for the highest k value to get an
  % initial guess. 
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1greenwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(2);
  vphinit(2)  = 2*UHPpole(2)-vphinit(23);

  vphinit(24)  = 2*UHPpole(3)-vphinit(3);
  vphinit(4)   = UHPpole(4)-conj(vphinit(3)-UHPpole(3));
  vphinit(25)  = 2*UHPpole(4)-vphinit(4);
  
  Ab=vphinit(27);
  vphinit(26)  = 2*UHPpole(5)-vphinit(5);
  vphinit(27)  = vphinit(10);
  vphinit(6)   = UHPpole(6)-conj(vphinit(5)-UHPpole(5));
  vphinit(27)  = 2*UHPpole(6)-vphinit(6);
  
  Ac=vphinit(8);Ad=vphinit(9);
  Ae=vphinit(29);Af=vphinit(30);
  vphinit(7)   = UHPpole(7) + (vphinit(5)-UHPpole(5));
  vphinit(28)  = 2*UHPpole(7)-vphinit(7);
  vphinit(8)   = UHPpole(8) + (vphinit(13)-UHPpole(16));
  vphinit(29)  = 2*UHPpole(8)-vphinit(8);
  vphinit(9)   = UHPpole(9) + (vphinit(6)-UHPpole(6));
  vphinit(30)  = 2*UHPpole(9)-vphinit(9);

  Ag=vphinit(11);Ah=vphinit(32);
  vphinit(11)  = vphinit(34);
  vphinit(32)  = 2*UHPpole(11)-vphinit(11);
  vphinit(10)  = UHPpole(10)-conj(vphinit(11)-UHPpole(11));
  vphinit(31)  = 2*UHPpole(10)-vphinit(10);

  Ai=vphinit(12);Aj=vphinit(13);Ak=vphinit(34);
  vphinit(12)  = Ah;
  vphinit(33)  = 2*UHPpole(12)-vphinit(12);
  vphinit(13)  = Ag;
  vphinit(34)  = 2*UHPpole(13)-vphinit(13);
  vphinit(14)  = UHPpole(14)-conj(vphinit(12)-UHPpole(12));
  vphinit(35)  = 2*UHPpole(14)-vphinit(14);

  vphinit(16)  = Aj;
  vphinit(37)  = 2*UHPpole(16)-vphinit(16);
  vphinit(15)  = UHPpole(15)-conj(vphinit(16)-UHPpole(16));
  vphinit(36)  = 2*UHPpole(15)-vphinit(15);

  vphinit(20)  = vphinit(41);
  vphinit(39)  = 2*UHPpole(18)-vphinit(18);
  vphinit(40)  = 2*UHPpole(19)-vphinit(19);
  vphinit(41)  = 2*UHPpole(20)-vphinit(20);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case1greenwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case1greenwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case1greenwe.mat
end


plot(hb(1),k/lambdaDe,real(omega(2,:))*fpe,'g')
axis(hb(1),[0 pi/lambdaDe 0 550])

plot(hb(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'g')
axis(hb(2),[0 pi/lambdaDe -6 90])
set(hb(1:2),'nextplot','replace')

aa=axis(hb(2));
xlabel(hb(2),'k   [m^{-1}]','fontname','times','fontsize',18)
ylabel(hb(1),'f   [Hz]','fontname','times','fontsize',18)
ylabel(hb(2),'\gamma   [s^{-1}]','fontname','times','fontsize',18)
set(hb(2),'nextplot','add')
plot(hb(2),aa(1:2),[0 0],'k--')
set(hb(1),'xticklabel',[],'ytick',[0:100:1000])
set(hb,'fontname','times','fontsize',14)



% add panel labels
for ii=1:2
  aa=axis(hb(ii));
  if strcmp(get(hb(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',hb(ii))
  else
    text(aa(1)*10^(0.05*diff(log10(aa(1:2)))), ...
         aa(3)*10^(0.9*diff(log10(aa(3:4)))),['(' char(96+ii) ')'], ...
         'fontname','times','fontsize',18,'parent',hb(ii))    
  end
end



% Case 2: 16:24 on 20 January 2015
hf=[];
figure(1624)
clf
set(gcf,'paperpositionmode','auto','renderer','zbuffer')
set(gcf,'position',[413 253 560 440])
hf(1)=axes('position',[0.14 0.565 0.83 0.43]);
hf(2)=axes('position',[0.14 0.13 0.83 0.43]);

% red curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which may also be
                   % used as reference.
ne=7.3519376e6;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));

if ~exist('matfiles/aa16case2redwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run it first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2redwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(4);
  vphinit(21) = 2*UHPpole(2)-vphinit(2);
  vphinit(4)  = 2*UHPpole(4)-vphinit(23);

  vphinit(6)  = Aa; Ab=vphinit(8);
  vphinit(25) = 2*UHPpole(6)-vphinit(6);
  vphinit(8) = 2*UHPpole(8)-vphinit(27);

  Ac=vphinit(9);
  vphinit(9)  = 2*UHPpole(9)-vphinit(28);
  vphinit(29) = 2*UHPpole(10)-vphinit(10);

  vphinit(31) = 2*UHPpole(12)-vphinit(12);
  vphinit(11) = UHPpole(11) -conj(vphinit(12)-UHPpole(12));
  vphinit(30) = 2*UHPpole(11)-vphinit(11);

  vphinit(32) = 2*UHPpole(13)-vphinit(13);
  vphinit(14) = UHPpole(14) -conj(vphinit(13)-UHPpole(13));
  vphinit(33) = 2*UHPpole(14)-vphinit(14);

  vphinit(35) = 2*UHPpole(16)-vphinit(16);
  vphinit(36) = 2*UHPpole(17)-vphinit(17);
  vphinit(37) = 2*UHPpole(18)-vphinit(18);
  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2redwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case2redwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case2redwe.mat
end

set(hf(1:2),'nextplot','add')
plot(hf(1),k/lambdaDe,real(omega(2,:))*fpe,'r')
plot(hf(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'r')
set(hf,'box','on')

% blue curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which may also be
                   % used as reference.
ne=2.6576e+08;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));
if ~exist('matfiles/aa16case2bluewe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run it first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2bluewe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(4);Ab=vphinit(17);
  vphinit(4)  = vphinit(5);
  vphinit(17) = 2*UHPpole(4)-vphinit(4);

  vphinit(6)  = 0.022239501617 + 0.002563881126683i;
  vphinit(19) = 2*UHPpole(6)-vphinit(6);
  vphinit(5)  = UHPpole(5) + conj(vphinit(6)-UHPpole(6));
  vphinit(18) = 2*UHPpole(5)-vphinit(5);

  vphinit(8)  = 0.04632493553291 + 0.00427328799294i;
  vphinit(21) = 2*UHPpole(8)-vphinit(8);
  vphinit(7)  = UHPpole(7) + conj(vphinit(8)-UHPpole(8));
  vphinit(20) = 2*UHPpole(7)-vphinit(7);

  vphinit(23) = 2*UHPpole(10)-vphinit(10);
  vphinit(24) = 2*UHPpole(11)-vphinit(11);
  vphinit(12) = vphinit(25);
  vphinit(25) = 2*UHPpole(12)-vphinit(12);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2bluewe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case2bluewe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case2bluewe.mat
end
plot(hf(1),k/lambdaDe,real(omega(2,:))*fpe,'b')
axis(hf(1),[0 pi/lambdaDe 0 1700])

plot(hf(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'b')
axis(hf(2),[0 pi/lambdaDe -41 550])



% green curve
kTe=7.9*elc;       % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which may also be
                   % used as reference.
ne=6.66497e7;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));

if ~exist('matfiles/aa16case2greenwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run it first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2greenwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(4);
  vphinit(16) = 2*UHPpole(3)-vphinit(3);
  vphinit(4) = 2*UHPpole(4)-vphinit(17);

  vphinit(6)  = vphinit(5);
  vphinit(19) = 2*UHPpole(6)-vphinit(6);
  vphinit(5)  = UHPpole(5) + conj(vphinit(6)-UHPpole(6));
  vphinit(18) = 2*UHPpole(5)-vphinit(5);

  vphinit(7)  = UHPpole(7) + (vphinit(5)-UHPpole(5));
  vphinit(20) = 2*UHPpole(7)-vphinit(7);
  vphinit(8)  = UHPpole(8) + conj(vphinit(7)-UHPpole(7));
  vphinit(21) = 2*UHPpole(8)-vphinit(8);

  vphinit(23) = 2*UHPpole(10)-vphinit(10);
  vphinit(24) = 2*UHPpole(11)-vphinit(11);
  vphinit(12) = vphinit(25);
  vphinit(25) = 2*UHPpole(12)-vphinit(12);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case2greenwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case2greenwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case2greenwe.mat
end

plot(hf(1),k/lambdaDe,real(omega(2,:))*fpe,'g')
axis(hf(1),[0 pi/lambdaDe 0 1050])

plot(hf(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'g')
axis(hf(2),[0 pi/lambdaDe -20 550])
set(hf(1:2),'nextplot','replace')

aa=axis(hf(2));
xlabel(hf(2),'k   [m^{-1}]','fontname','times','fontsize',18)
ylabel(hf(1),'f   [Hz]','fontname','times','fontsize',18)
ylabel(hf(2),'\gamma   [s^{-1}]','fontname','times','fontsize',18)
set(hf(2),'nextplot','add')
plot(hf(2),aa(1:2),[0 0],'k--')
set(hf(1),'xticklabel',[],'ytick',[0:200:1000])
set(hf(2),'ytick',[-1000:100:1000])
set(hf,'fontname','times','fontsize',14)



% add panel labels
for ii=1:2
  aa=axis(hf(ii));
  if strcmp(get(hf(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',hf(ii))
  else
    text(aa(1)*10^(0.05*diff(log10(aa(1:2)))), ...
         aa(3)*10^(0.9*diff(log10(aa(3:4)))),['(' char(96+ii) ')'], ...
         'fontname','times','fontsize',18,'parent',hf(ii))    
  end
end



% Case 3: 20:46 on 20 January 2015
he=[];
figure(2046)
clf
set(gcf,'paperpositionmode','auto','renderer','zbuffer')
set(gcf,'position',[797 253 560 440])
he(1)=axes('position',[0.14 0.565 0.83 0.43]);
he(2)=axes('position',[0.14 0.13 0.83 0.43]);

% red curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which also is
                   % used as reference.
ne=3.947577e+05;
omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));
if ~exist('matfiles/aa16case3redwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run it first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3redwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  vphinit(2)  = 2*UHPpole(2)-vphinit(18);

  Aa=vphinit(3);Ab=vphinit(4);
  vphinit(3)  = 2*UHPpole(3)-vphinit(19);
  vphinit(4)  = vphinit(6);
  vphinit(20) = 2*UHPpole(4)-vphinit(4);
  vphinit(6)  = vphinit(5);
  vphinit(22) = 2*UHPpole(6)-vphinit(6);
  vphinit(5)  = UHPpole(5)-conj(vphinit(6)-UHPpole(6));
  vphinit(21) = 2*UHPpole(5)-vphinit(5);

  Ac = vphinit(8);
  vphinit(7)  = vphinit(24);
  vphinit(8)  = vphinit(25);
  vphinit(24) = 2*UHPpole(8)-vphinit(8);
  
  Ad = vphinit(9);Ae = vphinit(10);
  vphinit(9) = Ac;
  vphinit(25) = 2*UHPpole(9)-vphinit(9);
  vphinit(10)  = UHPpole(10)-conj(vphinit(9)-UHPpole(9));
  vphinit(26) = 2*UHPpole(10)-vphinit(10);
  vphinit(27)  = Ae;

  vphinit(29)  = 2*UHPpole(13)-vphinit(13);
  vphinit(30)  = 2*UHPpole(14)-vphinit(14);
  vphinit(31)  = 2*UHPpole(15)-vphinit(15);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3redwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case3redwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case3redwe.mat
end

set(he(1:2),'nextplot','add')
plot(he(1),k/lambdaDe,real(omega(2,:))*fpe,'r')
plot(he(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'r')
set(he,'box','on')



% blue curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which also is
                   % used as reference.

ne=2.100761188e+08;

omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));
if ~exist('matfiles/aa16case3bluewe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run the old one first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3bluewe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(2);
  vphinit(2)  = 2*UHPpole(2)-vphinit(18);

  vphinit(3)  = UHPpole(3)+(vphinit(1)-UHPpole(1));
  vphinit(19) = 2*UHPpole(3)-vphinit(3);
  vphinit(4)  = UHPpole(4)+(vphinit(2)-UHPpole(2));
  vphinit(20) = 2*UHPpole(4)-vphinit(4);
  vphinit(5)  = UHPpole(5)+(vphinit(1)-UHPpole(1));
  vphinit(21) = 2*UHPpole(5)-vphinit(5);
  vphinit(6)  = UHPpole(6)+(vphinit(2)-UHPpole(2));
  vphinit(22) = 2*UHPpole(6)-vphinit(6);

  vphinit(7)  = UHPpole(7)+(vphinit(1)-UHPpole(1));
  vphinit(23) = 2*UHPpole(7)-vphinit(7);
  vphinit(8)  = UHPpole(8)+(vphinit(2)-UHPpole(2));
  vphinit(24) = 2*UHPpole(8)-vphinit(8);

  vphinit(9)  = UHPpole(9)+(vphinit(1)-UHPpole(1));
  vphinit(25) = 2*UHPpole(9)-vphinit(9);
  vphinit(10) = UHPpole(10)+(vphinit(2)-UHPpole(2));
  vphinit(26) = 2*UHPpole(10)-vphinit(10);

  vphinit(11) = UHPpole(11)-abs(vphinit(1)-UHPpole(1));
  vphinit(27) = 2*UHPpole(11)-vphinit(11);
 
  vphinit(29) = 2*UHPpole(13)-vphinit(13);
  vphinit(30) = 2*UHPpole(14)-vphinit(14);
  vphinit(31) = 2*UHPpole(15)-vphinit(15);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3bluewe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case3bluewe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case3bluewe.mat
end

plot(he(1),k/lambdaDe,real(omega(2,:))*fpe,'b')
axis(he(1),[0 pi/lambdaDe 0 1700])

plot(he(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'b')
axis(he(2),[0 pi/lambdaDe -35 550])



% green curve
kTe=12.9*elc;      % Electron temperature
cs=sqrt(kTe/mass); % This is the ion acoustic speed, which is used
                   % as reference.
vte=sqrt(kTe/me);  % This is the electron thermal speed, which also is
                   % used as reference.

ne=6.87621e7;

omega_pe=sqrt(ne*elc^2/(me*eps0));fpe=omega_pe/2/pi;
omega_pH2O=sqrt(ne*elc^2/(mass*eps0));fpH2O=omega_pH2O/2/pi;
lambdaDe=sqrt(eps0*kTe/(ne*elc^2));
if ~exist('matfiles/aa16case3greenwe.mat')
  vvv=[-3:0.005:7]*1e4/cs;
  op=optimset('MaxFunEvals',10000,'tolx',1e-12);
  % Run the old one first for the highest k value to get an initial guess
  [vphinit,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3greenwe',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(2);
  vphinit(2)  = vphinit(9);

  vphinit(6)  = vphinit(3);
  vphinit(22) = 2*UHPpole(6)-vphinit(6);  
  vphinit(5)  = UHPpole(5)-conj(vphinit(6)-UHPpole(6));
  vphinit(21) = 2*UHPpole(5)-vphinit(5);  
  vphinit(3)  = UHPpole(3)+(vphinit(5)-UHPpole(5));
  vphinit(19) = 2*UHPpole(3)-vphinit(3);
  vphinit(4)  = UHPpole(4)+(vphinit(6)-UHPpole(6));
  vphinit(20) = 2*UHPpole(4)-vphinit(4);

  Ab=vphinit(7);
  vphinit(7)  = vphinit(27);
  vphinit(23) = 2*UHPpole(7)-vphinit(7);
  vphinit(8)  = Ab;
  vphinit(24) = 2*UHPpole(8)-vphinit(8);

  vphinit(9)  = UHPpole(9) + (vphinit(3)-UHPpole(3));
  vphinit(25) = 2*UHPpole(9)-vphinit(9);
  vphinit(10) = UHPpole(10) + (vphinit(4)-UHPpole(4));
  vphinit(26) = 2*UHPpole(10)-vphinit(10);
  vphinit(11) = UHPpole(11) + abs(vphinit(4)-UHPpole(4));
  vphinit(27) = 2*UHPpole(11)-vphinit(11);

  vphinit(29) = 2*UHPpole(13)-vphinit(13);
  vphinit(30) = 2*UHPpole(14)-vphinit(14);
  vphinit(31) = 2*UHPpole(15)-vphinit(15);

  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16case3greenwe',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16case3greenwe.mat vph UHPpole lim omega f k v lim2 op
else
  load matfiles/aa16case3greenwe.mat
end

plot(he(1),k/lambdaDe,real(omega(2,:))*fpe,'g')
axis(he(1),[0 1.23*pi/lambdaDe 0 650])

plot(he(2),k/lambdaDe,imag(omega(2,:))*omega_pe,'g')
axis(he(2),[0 1.23*pi/lambdaDe 0 90])
set(he(1:2),'nextplot','replace')

aa=axis(hf(2));
xlabel(he(2),'k   [m^{-1}]','fontname','times','fontsize',18)
ylabel(he(1),'f   [Hz]','fontname','times','fontsize',18)
ylabel(he(2),'\gamma   [s^{-1}]','fontname','times','fontsize',18)
set(he(2),'nextplot','add')
plot(he(2),aa(1:2),[0 0],'k--')
set(he(1),'xticklabel',[],'ytick',[0:100:1000])
set(he(2),'ytick',[-100:20:200])
set(he,'fontname','times','fontsize',14)

% Root paths are not shown in the paper, but they may be instructive. 
% Here are those for the last computed dispersion relation.
figure(30)
plot(UHPpole,'ko')
hold on
plot(vph.','k')
plot(vph(2,:),'r')
hold off


if 3==4
  ccc=pwd;
  cd ~
  print -f800 -depsc2 fig-disp0800.eps
  print -f1624 -depsc2 fig-disp1624.eps
  print -f2046 -depsc2 fig-disp2046.eps
  cd(ccc)
end

