kineticdisp_jd

PURPOSE ^

SYNOPSIS ^

function [Nperp,e_xyz,e_pmz,e_pyk,D,X,X_s,X_sn,alphaphi_lin,phiT_xyz,phiP_xyz,phiT_pmz,phiP_pmz,phiT_pyk,phiP_pyk]= kineticdisp_jd(sa,sy,sbeta,Npar,Nperp_in,herm_mode,display_mode,nmax,NZ,tol,method)

DESCRIPTION ^

 This function solves the dispersion relation and the wave equation
 in the linear kinetic plasma model. It calculates the perpendicular index
 of refraction, the polarization, and the energy flow density as a 
 function of the wave frequency and the parallel index of refraction. 

 - (x,y,z): We assume that the magnetic field is in the z direction and 
         the perpendicular index of refraction in the x direction. 
 - (p,m,z): p and m refer to the left hand and right hand rotating fields 
         respectively.
 - (p,y,k): p and y are the transverse components of the polarization, 
         with respect to the wave vector. k is the longitudinal component.

   INPUTS: 

       - sa: ratio (wps/w)^2 for each species [Ns,1]
       - sy: ratio wcs/w for each species [Ns,1] 
           >>> Warning: sy < 0 for negatively charged species
       - sbeta: ratio (Ts/mc2)^(1/2) for each species [Ns,1]
       - Npar: parallel index of refraction [1,1]       
       - Nperp_in: perpendicular index of refraction (initial guess) [1,1] 
           >>> default: 1    
       - herm_mode: mode for the calculation: (0) full dielectric tensor 
                                              (1) hermitian part 
           >>> default: 0
           >>> In the case herm_mode = 1, the power flow density and
               the density of power absorbed are also calculated
       - display_mode: (0) no display (1) partial (2) figures  
           >>> default: 0
       - nmax: maximum harmonic number
           >>> default: 10
       - NZ: number of grid points in Z function integration calculation 
           >>> default: 1000
       - tol: tolerance on the value of the dispersion relation 
           >>> default: 1e-5
       - method_type: 'newton', 'segments'
           >>> default: 'newton'

   OUTPUTS:

       - Nperp: perpendicular index of refraction [1,1]
       - e_xyz: polarization vector in the (x,y,z) coordinates [3,1]
       - e_pmz: polarization vector in the (p,m,z) coordinates [3,1] 
       - e_pyk: polarization vector in the (p,y,k) coordinates [3,1]  
       - D: kinetic plasma dispersion tensor elements [3,3]
       - X: kinetic plasma susceptibility tensor elements [3,3]
       - X_s: kinetic plasma susceptibility tensor elements for species s 
           [3,3,Ns]
       - X_sn: kinetic plasma susceptibility tensor elements for species s 
           and harmonic n [3,3,Ns,Nn]

   if herm_mode = 1

       - alphaphi_lin: density of power absorbed [1,1]
           normalized to the (e0w/2)|E|^2
       - phiT_xyz: kinetic energy flow density,(x,y,z) coordinates [3,1]
       - phiP_xyz: Poynting energy flow density,(x,y,z) coordinates [3,1] 
       - phiT_pmz: kinetic energy flow density,(p,m,z) coordinates [3,1] 
       - phiP_pmz: Poynting energy flow density,(p,m,z) coordinates [3,1] 
       - phiT_pyk: kinetic energy flow density,(p,y,k) coordinates [3,1] 
       - phiP_pyk: Poynting energy flow density,(p,y,k) coordinates [3,1] 
           flows normalized to the vaccuum wave value (e0c/2)|E|^2

 created 05/02/2003
 by Joan Decker <joan.decker@cea.fr> and 
    Yves Peysson <yves.peysson@cea.fr>

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Nperp,e_xyz,e_pmz,e_pyk,D,X,X_s,X_sn,...
0002           alphaphi_lin,phiT_xyz,phiP_xyz,phiT_pmz,phiP_pmz,phiT_pyk,phiP_pyk]...
0003     = kineticdisp_jd(sa,sy,sbeta,Npar,Nperp_in,herm_mode,display_mode,nmax,NZ,tol,method)
0004 %
0005 % This function solves the dispersion relation and the wave equation
0006 % in the linear kinetic plasma model. It calculates the perpendicular index
0007 % of refraction, the polarization, and the energy flow density as a
0008 % function of the wave frequency and the parallel index of refraction.
0009 %
0010 % - (x,y,z): We assume that the magnetic field is in the z direction and
0011 %         the perpendicular index of refraction in the x direction.
0012 % - (p,m,z): p and m refer to the left hand and right hand rotating fields
0013 %         respectively.
0014 % - (p,y,k): p and y are the transverse components of the polarization,
0015 %         with respect to the wave vector. k is the longitudinal component.
0016 %
0017 %   INPUTS:
0018 %
0019 %       - sa: ratio (wps/w)^2 for each species [Ns,1]
0020 %       - sy: ratio wcs/w for each species [Ns,1]
0021 %           >>> Warning: sy < 0 for negatively charged species
0022 %       - sbeta: ratio (Ts/mc2)^(1/2) for each species [Ns,1]
0023 %       - Npar: parallel index of refraction [1,1]
0024 %       - Nperp_in: perpendicular index of refraction (initial guess) [1,1]
0025 %           >>> default: 1
0026 %       - herm_mode: mode for the calculation: (0) full dielectric tensor
0027 %                                              (1) hermitian part
0028 %           >>> default: 0
0029 %           >>> In the case herm_mode = 1, the power flow density and
0030 %               the density of power absorbed are also calculated
0031 %       - display_mode: (0) no display (1) partial (2) figures
0032 %           >>> default: 0
0033 %       - nmax: maximum harmonic number
0034 %           >>> default: 10
0035 %       - NZ: number of grid points in Z function integration calculation
0036 %           >>> default: 1000
0037 %       - tol: tolerance on the value of the dispersion relation
0038 %           >>> default: 1e-5
0039 %       - method_type: 'newton', 'segments'
0040 %           >>> default: 'newton'
0041 %
0042 %   OUTPUTS:
0043 %
0044 %       - Nperp: perpendicular index of refraction [1,1]
0045 %       - e_xyz: polarization vector in the (x,y,z) coordinates [3,1]
0046 %       - e_pmz: polarization vector in the (p,m,z) coordinates [3,1]
0047 %       - e_pyk: polarization vector in the (p,y,k) coordinates [3,1]
0048 %       - D: kinetic plasma dispersion tensor elements [3,3]
0049 %       - X: kinetic plasma susceptibility tensor elements [3,3]
0050 %       - X_s: kinetic plasma susceptibility tensor elements for species s
0051 %           [3,3,Ns]
0052 %       - X_sn: kinetic plasma susceptibility tensor elements for species s
0053 %           and harmonic n [3,3,Ns,Nn]
0054 %
0055 %   if herm_mode = 1
0056 %
0057 %       - alphaphi_lin: density of power absorbed [1,1]
0058 %           normalized to the (e0w/2)|E|^2
0059 %       - phiT_xyz: kinetic energy flow density,(x,y,z) coordinates [3,1]
0060 %       - phiP_xyz: Poynting energy flow density,(x,y,z) coordinates [3,1]
0061 %       - phiT_pmz: kinetic energy flow density,(p,m,z) coordinates [3,1]
0062 %       - phiP_pmz: Poynting energy flow density,(p,m,z) coordinates [3,1]
0063 %       - phiT_pyk: kinetic energy flow density,(p,y,k) coordinates [3,1]
0064 %       - phiP_pyk: Poynting energy flow density,(p,y,k) coordinates [3,1]
0065 %           flows normalized to the vaccuum wave value (e0c/2)|E|^2
0066 %
0067 % created 05/02/2003
0068 % by Joan Decker <joan.decker@cea.fr> and
0069 %    Yves Peysson <yves.peysson@cea.fr>
0070 %
0071 if nargin < 11,
0072    method = 'newton';
0073 end
0074 if nargin < 10,
0075    tol = 1e-5;
0076 end
0077 if nargin < 9,
0078    NZ = 1000;
0079 end
0080 if nargin < 8,
0081    nmax = 10;
0082 end
0083 if nargin < 7,
0084    display_mode = 0;
0085 end
0086 if nargin < 6,
0087    herm_mode = 0;
0088 end
0089 if nargin < 5,
0090    Nperp_in = 1;
0091 end
0092 if nargin < 4,
0093    error('Not enough arguments');
0094 end
0095 %
0096 % Internal parameters
0097 %
0098 nit_disp = 50;%maximum number of iterations in root finding procedure
0099 %
0100 proxcoef = 1e-4;%proximity parameter for local derivative calculations
0101 %(newton's method only)
0102 %
0103 mode = 0;%mode = 1 when root is localized
0104 varfac = 100;%reduction factor in exponential range expansion
0105 %
0106 % Initializations
0107 %
0108 sa = sa(:);
0109 sy = sy(:);
0110 sbeta = sbeta(:);
0111 %
0112 Ns = length(sa);
0113 %
0114 if length(sy) ~= Ns || length(sbeta) ~= Ns,
0115     error('Argument size not compatible')
0116 end
0117 %
0118 sy(sy.^2 == 1) = sy(sy.^2 == 1) + eps;%to avoid singularity
0119 %
0120 ss = sign(sy);
0121 spar = sign(Npar);
0122 %
0123 n = (-nmax:nmax);
0124 Nn = 2*nmax+1;
0125 n0 = nmax+1;
0126 %
0127 % Calculation of the plasma dispersion function
0128 %
0129 zeta_sn = (1 - abs(sy)*n)./(sqrt(2)*abs(Npar)*(sbeta*ones(1,Nn)));
0130 %
0131 [Z_sn,Zp_sn,Zpp_sn] = Zlasmadisp_jd(zeta_sn,display_mode,NZ);
0132 %
0133 %   Dielectric tensor elements
0134 %
0135 X_sn = zeros(3,3,Ns,Nn);
0136 Y_sn = zeros(3,3,Ns,Nn);
0137 %
0138 Msa_sn = shiftdim(repmat(sa,[1,Nn,3,3]),2);
0139 Mss_sn = reshape(repmat(ss,[1,Nn]),[1,1,Ns,Nn]);
0140 Mn_sn = reshape(repmat(n,[Ns,1]),[1,1,Ns,Nn]);
0141 %
0142 Mzeta_sn = reshape(zeta_sn,[1,1,Ns,Nn]);
0143 Mzeta0_sn = repmat(Mzeta_sn(:,:,:,n0),[3,3,1,Nn]);
0144 MZ_sn = reshape(Z_sn,[1,1,Ns,Nn]);
0145 MZp_sn = reshape(Zp_sn,[1,1,Ns,Nn]);
0146 %
0147 X_s = zeros(3,3,Ns);
0148 %
0149 X = zeros(3,3);
0150 D = zeros(3,3);
0151 %
0152 I = eye(3);
0153 %
0154 Drel = NaN*ones(1,2*nit_disp+1);
0155 Nperp = NaN*ones(1,2*(nit_disp+1));
0156 %
0157 Nperp(1) = Nperp_in;
0158 %
0159 for it_disp = 1:2*nit_disp+1,
0160     %
0161     N = [Nperp(it_disp),0,Npar].';
0162     %
0163     slambda = (Nperp(it_disp)*sbeta./sy).^2;
0164     Gamma_sn = besseli(n,slambda,1).*(exp(-i*imag(slambda))*ones(1,Nn));
0165     Gammap_sn = ((besseli(n+1,slambda,1) + besseli(n-1,slambda,1))/2 - besseli(n,slambda,1)).*(exp(-i*imag(slambda))*ones(1,Nn));
0166     %
0167     MGamma_sn = reshape(Gamma_sn,[1,1,Ns,Nn]);
0168     MGammap_sn = reshape(Gammap_sn,[1,1,Ns,Nn]);
0169     %
0170     Mlambda_sn = reshape(repmat(slambda,[1,Nn]),[1,1,Ns,Nn]);
0171     %
0172     Y_sn(1,1,:,:) = Mn_sn.^2.*MGamma_sn.*MZ_sn./Mlambda_sn;%Yxx
0173     Y_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*MGammap_sn.*MZ_sn;%Yxy
0174     Y_sn(1,3,:,:) = -spar.*Mn_sn.*MGamma_sn.*MZp_sn./sqrt(2*Mlambda_sn);%Yxz
0175     Y_sn(2,1,:,:) = -Y_sn(1,2,:,:);%Yyx
0176     Y_sn(2,2,:,:) = (Mn_sn.^2.*MGamma_sn./Mlambda_sn - 2*Mlambda_sn.*MGammap_sn).*MZ_sn;%Yyy
0177     Y_sn(2,3,:,:) = i*Mss_sn*spar.*sqrt(Mlambda_sn/2).*MGammap_sn.*MZp_sn;%Yyz
0178     Y_sn(3,1,:,:) = Y_sn(1,3,:,:);%Yzx
0179     Y_sn(3,2,:,:) = -Y_sn(2,3,:,:);%Yzy
0180     Y_sn(3,3,:,:) = -MGamma_sn.*Mzeta_sn.*MZp_sn;%Yzz
0181     %
0182     X_sn = Msa_sn.*Mzeta0_sn.*Y_sn;%Xsn
0183     %
0184     X_s = sum(X_sn,4);%Xs
0185     %
0186     X = sum(X_s,3);
0187     %
0188     D = N*N.' - (N.'*N)*I + I + X;
0189     %
0190     if herm_mode == 1,
0191         D = (D + D')/2;
0192     end
0193     %
0194     Drel(it_disp) = det(D); 
0195     %
0196     if strcmp(method,'newton'),
0197         if rem(it_disp,2) == 1 && abs(Drel(it_disp)) <= tol,
0198             break
0199         elseif rem(it_disp,2) == 1,
0200             if it_disp == 1,
0201                 Nperp(it_disp + 1) = (1 - proxcoef)*Nperp(it_disp);
0202             else
0203                 Nperp(it_disp + 1) = (1 - proxcoef)*Nperp(it_disp) + proxcoef*Nperp(it_disp - 2);
0204             end
0205         else
0206             dDreldNperp = (Drel(it_disp) - Drel(it_disp - 1))/(Nperp(it_disp) - Nperp(it_disp - 1));
0207             if dDreldNperp == 0 || isnan(dDreldNperp)
0208                 disp('WARNING : The calculation did not converge')
0209                 Nperp(it_disp + 1) = Nperp(it_disp - 1);
0210                 tol = Inf;
0211             else
0212                 Nperp(it_disp + 1) = Nperp(it_disp - 1) - Drel(it_disp - 1)/dDreldNperp; 
0213             end
0214         end
0215     elseif strcmp(method,'segments'),
0216         if abs(Drel(it_disp)) <= tol,
0217             break
0218         end
0219         if mode == 0,
0220             if it_disp > 1 && Drel(it_disp)*Drel(it_disp - 1) < 0,
0221                 mode = 1;
0222             elseif rem(it_disp,2) == 1,
0223                 Nperp(it_disp + 1) = Nperp(1)*exp((it_disp + 1)/2/varfac);
0224             else
0225                 Nperp(it_disp + 1) = Nperp(1)*exp(-it_disp/2/varfac);
0226             end
0227         end
0228         if mode == 1,
0229             if Drel(it_disp)*Drel(it_disp - 1) < 0,
0230                 it_ref = it_disp - 1;
0231             end
0232             Nperp(it_disp + 1) = (Nperp(it_ref) + Nperp(it_disp))/2;
0233         end
0234     else
0235         error('root finding method not recognized')
0236     end
0237 end
0238 %
0239 if display_mode >= 2,
0240 %
0241     figure(1),
0242     plot(1:(it_disp + 1)/2,real(Nperp(1:2:it_disp)),'b+',1:(it_disp + 1)/2,imag(Nperp(1:2:it_disp)),'r+')
0243 %
0244     figure(2),
0245     plot(1:(it_disp + 1)/2,real(Drel(1:2:it_disp)),'b+',1:(it_disp + 1)/2,imag(Drel(1:2:it_disp)),'r+')
0246 %
0247     figure(3),
0248     plot(1:(it_disp + 1)/2,abs(Drel(1:2:it_disp)),'b+')
0249     set(gca,'yscale','log')
0250 end
0251 %
0252 if it_disp == 2*nit_disp+1,    
0253     disp('WARNING : The root was not found');
0254 %
0255     Nperp = NaN;
0256     e_xyz = NaN*ones(3,1);
0257     e_pmz = NaN*ones(3,1);
0258     e_pyk = NaN*ones(3,1);
0259     D = NaN*ones(3,3);
0260     X = NaN*ones(3,3);
0261     X_s = NaN*ones(3,3,Ns);
0262     X_sn = NaN*ones(3,3,Ns,Nn);
0263     alphaphi_lin = NaN;
0264     phiT_xyz = NaN*ones(3,1);
0265     phiP_xyz = NaN*ones(3,1);
0266     phiT_pmz = NaN*ones(3,1);
0267     phiP_pmz = NaN*ones(3,1);
0268     phiT_pyk = NaN*ones(3,1);
0269     phiP_pyk = NaN*ones(3,1);
0270 %
0271     return
0272 end
0273 %
0274 Nperp = Nperp(it_disp);
0275 %
0276 %   Eigenvector
0277 %
0278 [Evec,Eval] = eig(D);
0279 %
0280 Eval0 = diag(Eval); %eigenvalues
0281 %
0282 iEval = find(abs(Eval0) == min(abs(Eval0)),1,'last'); %select eigenvalue lambda=0
0283 %
0284 Evec0 = Evec(:,iEval); %eigenvector
0285 %
0286 %   Polarization vectors
0287 %
0288 e_xyz = Evec0./norm(Evec0); %xyz
0289 %
0290 M_xypm = [1,i,0;1,-i,0;0,0,sqrt(2)]/sqrt(2);
0291 %
0292 e_pmz = M_xypm*e_xyz; %pmz
0293 %
0294 c = Npar/norm(N);
0295 s = Nperp/norm(N); 
0296 %
0297 M_xzpk = [c,0,-s;0,1,0;s,0,c];
0298 %
0299 e_pyk = M_xzpk*e_xyz; %pyk
0300 %
0301 if herm_mode == 1,
0302     %
0303     % linear absorption coefficient
0304     %
0305     XA = (X - X')/2;
0306     alphaphi_lin = -i*e_xyz'*XA*e_xyz;
0307     %
0308     % energy flow density
0309     %
0310     dYdNperp_sn = zeros(3,3,Ns,Nn);
0311     %
0312     Mlambdap_sn = 2*Mlambda_sn/Nperp;
0313     Gammapp_sn = ((besseli(n+2,slambda,1) + besseli(n-2,slambda,1))/4 + 3*besseli(n,slambda,1)/2 - besseli(n+1,slambda,1) - besseli(n-1,slambda,1)).*(exp(-i*imag(slambda))*ones(1,Nn));
0314     MGammapp_sn = reshape(Gammapp_sn,[1,1,Ns,Nn]);
0315     %
0316     dYdNperp_sn(1,1,:,:) = Mn_sn.^2.*Mlambdap_sn.*(MGammap_sn - MGamma_sn./Mlambda_sn).*MZ_sn./Mlambda_sn;
0317     dYdNperp_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*Mlambdap_sn.*MGammapp_sn.*MZ_sn;
0318     dYdNperp_sn(1,3,:,:) = -spar.*Mn_sn.*Mlambdap_sn.*(MGammap_sn - MGamma_sn/2./Mlambda_sn).*MZp_sn./sqrt(2*Mlambda_sn);
0319     dYdNperp_sn(2,1,:,:) = -dYdNperp_sn(1,2,:,:);
0320     dYdNperp_sn(2,2,:,:) = Mlambdap_sn.*((Mn_sn.^2./Mlambda_sn - 2).*MGammap_sn - Mn_sn.^2.*MGamma_sn./Mlambda_sn.^2 - 2*Mlambda_sn.*MGammapp_sn).*MZ_sn;
0321     dYdNperp_sn(2,3,:,:) = i*Mss_sn*spar.*Mlambdap_sn.*sqrt(Mlambda_sn/2).*(MGammap_sn/2./Mlambda_sn + MGammapp_sn).*MZp_sn;
0322     dYdNperp_sn(3,1,:,:) = dYdNperp_sn(1,3,:,:);
0323     dYdNperp_sn(3,2,:,:) = -dYdNperp_sn(2,3,:,:);
0324     dYdNperp_sn(3,3,:,:) = -Mlambdap_sn.*MGammap_sn.*Mzeta_sn.*MZp_sn;
0325     %
0326     dXdNperp_sn = Msa_sn.*Mzeta0_sn.*dYdNperp_sn;
0327     %
0328     dXdNperp_s = sum(dXdNperp_sn,4);
0329     %
0330     dXdNperp = sum(dXdNperp_s,3);
0331     %
0332     dXdNperp = (dXdNperp + dXdNperp')/2;%hermitian part only
0333     %
0334     dYdNpar_sn = zeros(3,3,Ns,Nn);
0335     %
0336     Mzetap_sn = -Mzeta_sn/Npar;
0337     Mzetap0_sn = -Mzeta0_sn/Npar;
0338     MZpp_sn = reshape(Zpp_sn,[1,1,Ns,Nn]);
0339     %
0340     dYdNpar_sn(1,1,:,:) = Mn_sn.^2.*MGamma_sn.*Mzetap_sn.*MZp_sn./Mlambda_sn;
0341     dYdNpar_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*MGammap_sn.*Mzetap_sn.*MZp_sn;
0342     dYdNpar_sn(1,3,:,:) = -spar.*Mn_sn.*MGamma_sn.*Mzetap_sn.*MZpp_sn./sqrt(2*Mlambda_sn);
0343     dYdNpar_sn(2,1,:,:) = -dYdNpar_sn(1,2,:,:);
0344     dYdNpar_sn(2,2,:,:) = (Mn_sn.^2.*MGamma_sn./Mlambda_sn - 2*Mlambda_sn.*MGammap_sn).*Mzetap_sn.*MZp_sn;
0345     dYdNpar_sn(2,3,:,:) = i*Mss_sn*spar.*sqrt(Mlambda_sn/2).*MGammap_sn.*Mzetap_sn.*MZpp_sn;
0346     dYdNpar_sn(3,1,:,:) = dYdNpar_sn(1,3,:,:);
0347     dYdNpar_sn(3,2,:,:) = -dYdNpar_sn(2,3,:,:);
0348     dYdNpar_sn(3,3,:,:) = -MGamma_sn.*Mzetap_sn.*(MZp_sn + Mzeta_sn.*MZpp_sn);
0349     %
0350     dXdNpar_sn = Msa_sn.*Mzetap0_sn.*Y_sn + Msa_sn.*Mzeta0_sn.*dYdNpar_sn;
0351     %
0352     dXdNpar_s = sum(dXdNpar_sn,4);
0353     %
0354     dXdNpar = sum(dXdNpar_s,3);
0355     %
0356     dXdNpar  = (dXdNpar + dXdNpar')/2;
0357     %
0358     phiT_xyz(1,1) = -e_xyz'*dXdNperp*e_xyz/2;
0359     phiT_xyz(2,1) = 0;
0360     phiT_xyz(3,1) = -e_xyz'*dXdNpar*e_xyz/2;
0361     phiP_xyz = N - real((e_xyz'*N)*e_xyz);
0362     %
0363     phiT_pmz = M_xypm*phiT_xyz; %pmz
0364     phiP_pmz = M_xypm*phiP_xyz; %pmz
0365     %
0366     phiT_pyk = M_xzpk*phiT_xyz; %pyk
0367     phiP_pyk = M_xzpk*phiP_xyz; %pyk
0368 else
0369     alphaphi_lin = NaN;
0370     phiT_xyz = NaN;
0371     phiP_xyz = NaN;
0372     phiT_pmz = NaN;
0373     phiP_pmz = NaN;
0374     phiT_pyk = NaN;
0375     phiP_pyk = NaN;
0376 end
0377 
0378 
0379 
0380 
0381 
0382 
0383 
0384 
0385
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.