proc_ES_ebw_jd

 Non-relativistic electrostatic dispersion processing for high-frequency waves

 Note: xalphaphi_lin is normalized to (omega/clum)

   INPUTS:
       
       - xa: (wpe/w)^2 [1,nx]
       - xq: w/wce [1,nx]
       - xbeta: sqrt(Te/mc2) [1,nx]
       - xNpar: parallel index of refraction [1,nx]
       - xNperp_in: guess for root [1,nx] or [1,1] (with progressive guess)
       - opt_disp: option for dispersion relation: 
           - (0): exact ES NR 
           - (1): Z function expansion large argument 2nd order 
           - (2): Z function expansion large argument 1st order + Nperp >> Npar

       - nmax: number of harmonic numbers to be considered

 J. Decker 14/12/05 (joan.decker@cea.fr)

0001 function [mask,xNperp,xepol_pmz,xphi_xyz,xalphaphi_lin_nr,...
0002         xalphaphi_lin_wr,xalphaphi_lin_fr,xepol_xyz,xepol_pyk,xphi_pmz,xphi_pyk] = ...
0003     proc_ES_ebw_jd(xa,xq,xbeta,xNpar,xNperp,opt_disp_in,nmax_in)
0004 %
0005 % Non-relativistic electrostatic dispersion processing for high-frequency waves
0006 %
0007 % Note: xalphaphi_lin is normalized to (omega/clum)
0008 %
0009 %   INPUTS:
0010 %
0011 %       - xa: (wpe/w)^2 [1,nx]
0012 %       - xq: w/wce [1,nx]
0013 %       - xbeta: sqrt(Te/mc2) [1,nx]
0014 %       - xNpar: parallel index of refraction [1,nx]
0015 %       - xNperp_in: guess for root [1,nx] or [1,1] (with progressive guess)
0016 %       - opt_disp: option for dispersion relation:
0017 %           - (0): exact ES NR
0018 %           - (1): Z function expansion large argument 2nd order
0019 %           - (2): Z function expansion large argument 1st order + Nperp >> Npar
0020 %
0021 %       - nmax: number of harmonic numbers to be considered
0022 %
0023 % J. Decker 14/12/05 (joan.decker@cea.fr)
0024 %
0025 if nargin < 7
0026     nmax_in = 50;
0027 end
0028 if nargin < 6
0029     opt_disp_in = 0;
0030 end
0031 if nargin < 5
0032     error('An initial guess must be provided')
0033 end
0034 %
0035 opt_disp = opt_disp_in;
0036 nmax = nmax_in;
0037 %
0038 nx = max([length(xa),length(xq),length(xbeta),length(xNpar),length(xNperp)]);
0039 %
0040 if length(xa) ~= nx,
0041     error('Argument dimension mismatch')
0042 end
0043 %
0044 if length(xq) ~= nx,
0045     error('Argument dimension mismatch')
0046 end
0047 %
0048 if length(xbeta) ~= nx,
0049     error('Argument dimension mismatch')
0050 end
0051 %
0052 if length(xNpar) ~= nx,
0053     error('Argument dimension mismatch')
0054 end
0055 %
0056 if length(xNperp) ~= nx,
0057     error('Argument dimension mismatch')
0058 end
0059 %
0060 xa = xa(:).';
0061 xq = xq(:).';
0062 xbeta = xbeta(:).';
0063 xNpar = xNpar(:).';
0064 xNperp = xNperp(:).';
0065 %
0066 % --------------- Initialization -----------------
0067 %
0068 xalphaphi_lin_nr = NaN*ones(1,nx);
0069 xalphaphi_lin_wr = NaN*ones(1,nx);
0070 xalphaphi_lin_fr = NaN*ones(1,nx);
0071 %
0072 xepol_xyz = NaN*ones(3,nx);
0073 xepol_pmz = NaN*ones(3,nx);
0074 xepol_pyk = NaN*ones(3,nx);
0075 %
0076 xphi_xyz = NaN*ones(3,nx);
0077 xphi_pmz = NaN*ones(3,nx);
0078 xphi_pyk = NaN*ones(3,nx);
0079 %
0080 xlpar = (xNpar.*xbeta.*xq).^2;
0081 %
0082 xNperp(real(xNperp) < 1e-2) = NaN;
0083 xlperp = (xNperp.*xbeta.*xq).^2;
0084 %
0085 mask = ~isnan(xNperp); 
0086 nx = length(mask);
0087 %
0088 xa = xa(mask);
0089 xq = xq(mask);
0090 xbeta = xbeta(mask);
0091 xNpar = xNpar(mask);
0092 xNperp = xNperp(mask);
0093 %
0094 % ---------------- Polarization ------------------
0095 %
0096 xepol_pyk(1,mask) = 0;
0097 xepol_pyk(2,mask) = 0;
0098 xepol_pyk(3,mask) = 1;
0099 %
0100 xNtot = sqrt(xNperp.^2 + xNpar.^2);
0101 %
0102 xepol_xyz(1,mask) = xNperp./xNtot;
0103 xepol_xyz(2,mask) = 0./xNtot;
0104 xepol_xyz(3,mask) = xNpar./xNtot;
0105 %
0106 xepol_pmz(1,mask) = xNperp./xNtot/sqrt(2);
0107 xepol_pmz(2,mask) = xNperp./xNtot/sqrt(2);
0108 xepol_pmz(3,mask) = xNpar./xNtot;
0109 %
0110 % ------------- Energy flow density --------------
0111 %
0112 n = -(nmax+1):nmax+1;
0113 nn = 2*nmax+1;
0114 %
0115 XBessel = besseli(n.',xlperp,1).'; % takes the Bessel function on lambda with orders 1:N
0116 XGamma = XBessel(2:end-1,:);
0117 XdGammadlperp = (XBessel(1:end-2,:) + XBessel(3:end,:))/2 - XGamma;
0118 %
0119 xdlperpdNperp = 2*xNperp.*(xbeta.*xq).^2;
0120 xdlpardNpar = 2*xNpar.*(xbeta.*xq).^2;
0121 %
0122 Xn = repmat(n(2:end-1).',[1,nx]);
0123 %
0124 Xq = repmat(xq,[nn,1]);
0125 Xlperp = repmat(xlperp,[nn,1]);
0126 Xlpar = repmat(xlpar,[nn,1]);
0127 %
0128 if opt_disp == 0,
0129     %
0130     Xzeta = (Xq-Xn)./sqrt(2*Xlpar);
0131     %
0132     [XZ,XZp] = Zlasmadisp_jd(Xzeta);
0133     %
0134     xdYdNperp = xdlperpdNperp.*xa.*xq.^2./(xlperp + xlpar).*...
0135         sum((XdGammadlperp - XGamma./(Xlperp + Xlpar)).*(1 + Xn./Xq + Xq.*real(XZ)./sqrt(2*Xlpar)),1);
0136     xdYdNpar = - xdlpardNpar.*xa.*xq.^2./(xlperp + xlpar).*...
0137         sum(XGamma./(Xlperp + Xlpar).*(1 + Xn./Xq + Xq.*real(XZ)./sqrt(2*Xlpar)) + XGamma.*Xq.*(real(XZ) + Xzeta.*real(XZp))./(2*Xlpar).^(3/2),1);
0138 elseif opt_disp == 1,
0139     %
0140     xdYdNperp = xdlperpdNperp.*xa./(xlperp + xlpar).*...
0141         sum((XGamma./(Xlperp + Xlpar) - XdGammadlperp).*(Xn.^2./(1 - Xn./Xq) + Xlpar./(1 - Xn./Xq).^3),1);
0142     %
0143     xdYdNpar = xdlpardNpar.*xa./(xlperp + xlpar).*...
0144         sum(XGamma./(Xlperp + Xlpar).*(Xn.^2./(1 - Xn./Xq) + Xlpar./(1 - Xn./Xq).^3) - XGamma./(1 - Xn./Xq).^3,1);
0145     %
0146 elseif opt_disp == 2,
0147     %
0148     xdYdNperp = xdlperpdNperp.*xa./xlperp.*...
0149         sum((XGamma./Xlperp - XdGammadlperp).*Xn.^2./(1 - Xn./Xq),1);
0150     %
0151     xdYdNpar = zeros(1,nx);
0152     %
0153 end
0154 %
0155 xphi_xyz(1,mask) = - xdYdNperp/2;
0156 xphi_xyz(2,mask) = 0;
0157 xphi_xyz(3,mask) = - xdYdNpar/2;
0158 %
0159 % ----------- Absorption Coefficient -------------
0160 %
0161 if nargout > 4,
0162     %
0163     [xalphaphi_lin_nr(mask),xalphaphi_lin_wr(mask),xalphaphi_lin_fr(mask)] = alphaphi_ebw_jd(xa,xq,xbeta,xNpar,xNperp,opt_disp,nmax);
0164     %
0165 end
0166