% aa16fig6 is a program intended to produce figure showing
% dispersion relations for one case where the solar wind ions are
% included in order to show that these ions do not matter. The figure is
% intended for use in a paper in Astronomy and Astrophysics in 2016. 
%
% The case used is what is called case 1 blue curve in the other files
% and corresponds to the ICA sweep obtained at 08:00 on 20 January 2015.
%
% HG 2016-09-09


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

% Case 1: 08:00 on 20 January 2015
nprot  = 5.334635e4;
nalpha = 6.76553e4;


% blue curve
kTe=7*elc;         % Electron temperature
cs=sqrt(kTe/mH2O); % 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/(mH2O*eps0));fpH2O=omega_pH2O/2/pi;

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

vvv=[-3:0.005:6]*1e5/cs;

aa16case1blue;
fionblue=sum(ALLresid.'*ones(size(vvv))./ ...
             (ones(size(ALLpole)).'*vvv - ALLpole.'*ones(size(vvv))));
fionblue=real(fionblue);
vblue=vvv;

aa16WithSWprotons;
fionp=sum(ALLresid.'*ones(size(vvv))./ ...
          (ones(size(ALLpole)).'*vvv - ALLpole.'*ones(size(vvv))));
fionp=real(fionp);
vp=vvv;

aa16WithSWalphas;
fiona=sum(ALLresid.'*ones(size(vvv))./ ...
          (ones(size(ALLpole)).'*vvv - ALLpole.'*ones(size(vvv))));
fiona=real(fiona);
va=vvv;


if ~exist('matfiles/aa16WithSWall.mat')
  % 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('aa16WithSWall',2.^(8.3),vvv,op);
  % Then we update that guess manually here
  Aa=vphinit(2);
  vphinit(2)  = vphinit(16);

  vphinit(3)  = UHPpole(3)+(vphinit(1)-UHPpole(1));
  vphinit(31) = UHPpole(3)+(vphinit(29)-UHPpole(1));
  vphinit(4)  = UHPpole(4)+(vphinit(2)-UHPpole(2));
  vphinit(32) = UHPpole(4)+(vphinit(30)-UHPpole(2));
  
  Ab=vphinit(6);
  vphinit(33)  = vphinit(44);
  vphinit(6)  = vphinit(36);
  vphinit(34)  = 2*UHPpole(6)-vphinit(6);

  Ac=vphinit(8);Ad=vphinit(36);Ae=vphinit(37);
  vphinit(7)   = UHPpole(7) + (vphinit(5)-UHPpole(5));
  vphinit(35)  = 2*UHPpole(7)-vphinit(7);
  vphinit(8)   = UHPpole(8) + abs(vphinit(5)-UHPpole(5));
  vphinit(36)  = 2*UHPpole(8)-vphinit(8);
  vphinit(9)   = UHPpole(9) + (vphinit(6)-UHPpole(6));
  vphinit(37)  = 2*UHPpole(9)-vphinit(9);

  Af=vphinit(10);Ag=vphinit(11);Ah=vphinit(39);
  vphinit(11)  = vphinit(10);
  vphinit(39)  = 2*UHPpole(11)-vphinit(11);
  vphinit(10)  = UHPpole(10)-conj(vphinit(11)-UHPpole(11));
  vphinit(38)  = 2*UHPpole(10)-vphinit(10);

  Ai=vphinit(12);Aj=vphinit(13);
  vphinit(12)  = Ah;
  vphinit(40)  = 2*UHPpole(12)-vphinit(12);
  vphinit(13)  = Ag;
  vphinit(41)  = 2*UHPpole(13)-vphinit(13);
  vphinit(14)  = UHPpole(14)-conj(vphinit(12)-UHPpole(12));
  vphinit(42)  = 2*UHPpole(14)-vphinit(14);

  vphinit(16)  = Aj;
  vphinit(44)  = 2*UHPpole(16)-vphinit(16);
  vphinit(15)  = UHPpole(15)-conj(vphinit(16)-UHPpole(16));
  vphinit(43)  = 2*UHPpole(15)-vphinit(15);

  Ak=vphinit(17);Al=vphinit(18);Am=vphinit(19);An=vphinit(20);
  Ao=vphinit(45);Ap=vphinit(46);Aq=vphinit(48);
  vphinit(17)  = UHPpole(17) + (vphinit(15)-UHPpole(15));
  vphinit(45)  = 2*UHPpole(17)-vphinit(17);
  vphinit(18)  = UHPpole(18)-conj(vphinit(17)-UHPpole(17));
  vphinit(46)  = 2*UHPpole(18)-vphinit(18);
  vphinit(19)  = UHPpole(19) + (vphinit(17)-UHPpole(17));
  vphinit(47)  = 2*UHPpole(19)-vphinit(19);
  vphinit(20)  = UHPpole(20) + (vphinit(18)-UHPpole(18));
  vphinit(48)  = 2*UHPpole(20)-vphinit(20);

  vphinit(21)  = UHPpole(21) + (vphinit(17)-UHPpole(17));
  vphinit(49)  = 2*UHPpole(21)-vphinit(21);
  vphinit(22)  = UHPpole(22)-conj(vphinit(21)-UHPpole(21));
  vphinit(50)  = 2*UHPpole(22)-vphinit(22);
  vphinit(23)  = UHPpole(23) +abs(vphinit(21)-UHPpole(21));
  vphinit(51)  = 2*UHPpole(23)-vphinit(23);
  
  vphinit(53)  = 2*UHPpole(25)-vphinit(25);
  vphinit(54)  = 2*UHPpole(26)-vphinit(26);
  vphinit(55)  = Aa;
  
  
  [vph,UHPpole,lim,omega,f,k,v,lim2] = ...
      disprelelectronnew('aa16WithSWall',2.^[-6:0.01:8.3],vvv,op,vphinit);
  save matfiles/aa16WithSWall.mat vph ...
      UHPpole lim omega f k v lim2 op vvv
else
  load matfiles/aa16WithSWall.mat
end

hb=[];
figure(801)
clf
set(gcf,'paperpositionmode','auto','renderer','zbuffer')
set(gcf,'position',[133 410 560 560])
hb(1)=axes('position',[0.12 0.115 0.35 0.88]);
hb(2)=axes('position',[0.60 0.555 0.39 0.44]);
hb(3)=axes('position',[0.60 0.115 0.39 0.435]);

semilogy(hb(1),vblue*cs*1e-3,fionblue*ne/cs,'b', ...
         va*cs*1e-3,fiona*nalpha/cs,'k--')
set(hb(1),'nextplot','add')
semilogy(hb(1),vp*cs*1e-3,fionp*nprot/cs,'-.','color',[1 127/255 0])
set(hb(1),'nextplot','replace')
axis(hb(1),[-200 500 1e-5 9e5])

loglog(hb(2),k/lambdaDe,real(omega(2,:))*fpe,'b', ...
     k/lambdaDe,real(omega(18,:))*fpe,'k--')
set(hb(2),'nextplot','add')
loglog(hb(2),k/lambdaDe,real(omega(22,:))*fpe,'-.','color',[1 127/255 0])
set(hb(2),'nextplot','replace')
axis(hb(2),[0.01 1.1*pi/lambdaDe 8 4e4])

loglog(hb(3),k/lambdaDe,imag(omega(2,:))*omega_pe,'b', ...
     k/lambdaDe,imag(omega(18,:))*omega_pe,'k--')
set(hb(3),'nextplot','add')
loglog(hb(3),k/lambdaDe,imag(omega(22,:))*omega_pe,'-.','color',[1 127/255 0])
set(hb(3),'nextplot','replace')
axis(hb(3),[0.01 1.1*pi/lambdaDe 1e-1 4e4])

xlabel(hb(1),'v   [km/s]','fontname','times','fontsize',18)
ylabel(hb(1),'f(v)   [m^{-4}s]','fontname','times','fontsize',18)
xlabel(hb(3),'k   [m^{-1}]','fontname','times','fontsize',18)
ylabel(hb(2),'f   [Hz]','fontname','times','fontsize',18)
ylabel(hb(3),'\gamma   [s^{-1}]','fontname','times','fontsize',18)
set(hb(2),'xticklabel',[])
set(hb(1),'ytick',[10.^[-10:10]])
set(hb(2),'ytick',[10.^[-1:5]])
set(hb(3),'ytick',[10.^[-1:5]])
set(hb(2),'xtick',[10.^[-2:5]])
set(hb(3),'xtick',[10.^[-2:5]])
set(hb,'fontname','times','fontsize',14)

hl=legend(hb(3),'H_{2}O^{+}','He^{2+}','H^{+}');
set(hl,'position',[0.31 0.747 0.1723 0.2483])

% add panel labels
for ii=1:length(hb)
  aa=axis(hb(ii));
  if strcmp(get(hb(ii),'xscale'),'linear')
    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',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


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

