colddisp_dke_jd

 This function solves the dispersion relation and the wave equation
 in the linear cold plasma model. It calculates the perpendicular index of
 refraction, the polarization and the power flow 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.

 The power flow is: phi = Re[e x (N x e*)] = Re[N - (N . e)e*] (electromagnetic approximation only) 

   INPUTS: 

       - Npar: parallel index of refraction [1,1]       
       - omega_rf: wave frequency (rad/s) [1,1] 
       - ne: electron density (m-3) [1,1]
       - zni: ion density (m-3) [1,i]
       - zZi: ion charges [1,i] 
       - zmi: ion masses (mp) [1,i]  
       - Bfield: magnetic field amplitude [1,1]
       - Nperp: Perpendicular refractive index (optional) [1,1]

   OUTPUTS:

       - Nperpp: perpendicular index of refraction (root +) [1,1]
       - Nperpm: perpendicular index of refraction (root -) [1,1] 
       - Kxyz: cold plasma dielectric tensor elements [3,3]
       - ep_xyz: polarization vector in the (x,y,z) coordinates (root +) [3,1]
       - em_xyz: polarization vector in the (x,y,z) coordinates (root -) [3,1] 
       - ep_pmz: polarization vector in the (p,m,z) coordinates (root +) [3,1] 
       - em_pmz: polarization vector in the (p,m,z) coordinates (root -) [3,1]  
       - ep_pyk: polarization vector in the (p,y,k) coordinates (root +) [3,1]  
       - em_pyk: polarization vector in the (p,y,k) coordinates (root -) [3,1] 
       - phip_xyz: normalized power flow vector in the (x,y,z) coordinates (root +) [3,1]
       - phim_xyz: normalized power flow vector in the (x,y,z) coordinates (root -) [3,1] 
       - phip_pmz: normalized power flow vector in the (p,m,z) coordinates (root +) [3,1] 
       - phim_pmz: normalized power flow vector in the (p,m,z) coordinates (root -) [3,1]  
       - phip_pyk: normalized power flow vector in the (p,y,k) coordinates (root +) [3,1]  
       - phim_pyk: normalized power flow vector in the (p,y,k) coordinates (root -) [3,1] 
       - sigmap: normalized energy density (root +) [1,1] 
       - sigmam: normalized energy density (root -) [1,1] 

 created 05/02/2003
 by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and  Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)

0001 function [Nperpp,Nperpm,Kxyz,...
0002     ep_xyz,em_xyz,ep_pmz,em_pmz,ep_pyk,em_pyk,...
0003     phip_xyz,phim_xyz,phip_pmz,phim_pmz,phip_pyk,phim_pyk,...
0004     sigmap,sigmam,Dp,Dm,ps,ys]...
0005     = colddisp_dke_jd(Npar,omega_rf,ne,zni,zZi,zmi,Bfield,Nperp)
0006 %
0007 % This function solves the dispersion relation and the wave equation
0008 % in the linear cold plasma model. It calculates the perpendicular index of
0009 % refraction, the polarization and the power flow as a function of the wave
0010 % frequency and the parallel index of refraction.
0011 %
0012 % - (x,y,z): We assume that the magnetic field is in the z direction and the
0013 %           perpendicular index of refraction in the x direction.
0014 % - (p,m,z): p and m refer to the left hand and right hand rotating fields
0015 %           respectively.
0016 % - (p,y,k): p and y are the transverse components of the polarization, with
0017 %           respect to the wave vector. k is the longitudinal component.
0018 %
0019 % The power flow is: phi = Re[e x (N x e*)] = Re[N - (N . e)e*] (electromagnetic approximation only)
0020 %
0021 %   INPUTS:
0022 %
0023 %       - Npar: parallel index of refraction [1,1]
0024 %       - omega_rf: wave frequency (rad/s) [1,1]
0025 %       - ne: electron density (m-3) [1,1]
0026 %       - zni: ion density (m-3) [1,i]
0027 %       - zZi: ion charges [1,i]
0028 %       - zmi: ion masses (mp) [1,i]
0029 %       - Bfield: magnetic field amplitude [1,1]
0030 %       - Nperp: Perpendicular refractive index (optional) [1,1]
0031 %
0032 %   OUTPUTS:
0033 %
0034 %       - Nperpp: perpendicular index of refraction (root +) [1,1]
0035 %       - Nperpm: perpendicular index of refraction (root -) [1,1]
0036 %       - Kxyz: cold plasma dielectric tensor elements [3,3]
0037 %       - ep_xyz: polarization vector in the (x,y,z) coordinates (root +) [3,1]
0038 %       - em_xyz: polarization vector in the (x,y,z) coordinates (root -) [3,1]
0039 %       - ep_pmz: polarization vector in the (p,m,z) coordinates (root +) [3,1]
0040 %       - em_pmz: polarization vector in the (p,m,z) coordinates (root -) [3,1]
0041 %       - ep_pyk: polarization vector in the (p,y,k) coordinates (root +) [3,1]
0042 %       - em_pyk: polarization vector in the (p,y,k) coordinates (root -) [3,1]
0043 %       - phip_xyz: normalized power flow vector in the (x,y,z) coordinates (root +) [3,1]
0044 %       - phim_xyz: normalized power flow vector in the (x,y,z) coordinates (root -) [3,1]
0045 %       - phip_pmz: normalized power flow vector in the (p,m,z) coordinates (root +) [3,1]
0046 %       - phim_pmz: normalized power flow vector in the (p,m,z) coordinates (root -) [3,1]
0047 %       - phip_pyk: normalized power flow vector in the (p,y,k) coordinates (root +) [3,1]
0048 %       - phim_pyk: normalized power flow vector in the (p,y,k) coordinates (root -) [3,1]
0049 %       - sigmap: normalized energy density (root +) [1,1]
0050 %       - sigmam: normalized energy density (root -) [1,1]
0051 %
0052 % created 05/02/2003
0053 % by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and  Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0054 %
0055 [qe,me,mp,mn,e0] = pc_dke_yp;
0056 %
0057 if isnan(ne),%vacuum limit
0058     %
0059     ne = 1e-5*me*e0*omega_rf^2/qe^2;%so that the polarizations are separated
0060     %
0061 end
0062 %
0063 ns = [ne;zni(:)];
0064 qs = qe*[-1;zZi(:)];    
0065 if zZi(1) == 0,
0066    qs = qs*0;%vacuum case, no charge
0067 end
0068 %
0069 ms = [me;mp*zmi(:)];
0070 %
0071 wps = sqrt(ns.*qs.^2/e0./ms);
0072 wcs = qs*Bfield./ms;
0073 %
0074 ps = wps/omega_rf;
0075 ys = wcs/omega_rf;
0076 %
0077 ys(ys.^2 == 1) = ys(ys.^2 == 1) + eps;
0078 %
0079 %   Dielectric tensor elements
0080 %
0081 Kperp = 1 - sum(ps.^2./(1 - ys.^2));
0082 Kpar = 1 - sum(ps.^2);
0083 Kxy = sum(ps.^2.*ys./(1 - ys.^2));
0084 %
0085 Kxyz = [Kperp,-i*Kxy,0;i*Kxy,Kperp,0;0,0,Kpar];
0086 %
0087 A = Kperp;
0088 B = (Npar^2 - Kperp)*(Kperp + Kpar) + Kxy^2;
0089 C = Kpar*((Npar^2 - Kperp)^2 - Kxy^2);
0090 %
0091 %   Accessibility condition
0092 %
0093 if isnan(real(Npar)),
0094     Npar = sqrt(Kperp + (sqrt(4*Kperp*Kpar*Kxy^2*(Kxy^2 - (Kperp - Kpar)^2)) - Kxy^2*(Kperp + Kpar))/(Kperp - Kpar)^2);%warning: this solution is valid in the LH frequency range only
0095     Nperpp = Npar;%return accessibility condition
0096     return
0097 elseif isnan(imag(Npar)),
0098 %
0099 %   Cutoff condition
0100 %
0101     Npar = sqrt(Kperp + Kxy);
0102     %
0103     Nperpp = Npar;%return accessibility condition
0104     return    
0105 end
0106 %
0107 %   Dispersion relation
0108 %
0109 if nargin == 8,%Nperp is an input variable
0110     Nxyzp = [Nperp;0;Npar];
0111     Nxyzm = [Nperp;0;Npar];
0112     Nperpp = Nperp;
0113     Nperpm = Nperp;
0114 elseif nargin == 7,%Nperp is calculated from the cold dispersion relation
0115     Nperpp = sqrt((-B + sqrt(B^2-4*A*C))/(2*A));
0116     Nperpm = sqrt((-B - sqrt(B^2-4*A*C))/(2*A));
0117     %
0118     % Note: root + corresponds to the slow LH wave when Kperp > 0
0119     %
0120     Nxyzp = [Nperpp;0;Npar];         
0121     Nxyzm = [Nperpm;0;Npar]; 
0122 else
0123    error('Wrong number of input arguments in colddisp_dke_jd');
0124 end
0125 %
0126 %   Dispersion tensor for + root
0127 %
0128 N2p = Nperpp^2+Npar^2;
0129 N2m = Nperpm^2+Npar^2;
0130 %
0131 Dxyzp = Nxyzp*Nxyzp.'/N2p - sum(Nxyzp.^2/N2p)*diag([1,1,1]) + Kxyz/N2p;
0132 %
0133 %   Dispersion tensor for - root
0134 %
0135 Dxyzm = Nxyzm*Nxyzm.'/N2m - sum(Nxyzm.^2/N2m)*diag([1,1,1]) + Kxyz/N2m;
0136 %
0137 % verification (useful if nargin == 8)
0138 %
0139 Dp = det(Dxyzp);
0140 Dm = det(Dxyzm);
0141 %
0142 %   Eigenvectors
0143 %
0144 [Evecp,Evalp] = eig(Dxyzp);
0145 [Evecm,Evalm] = eig(Dxyzm);
0146 %
0147 Eval0p = diag(Evalp); %eigenvalues
0148 Eval0m = diag(Evalm); %eigenvalues
0149 %
0150 jp = find(abs(Eval0p) == min(abs(Eval0p))); %select eigenvalue lambda=0
0151 jm = find(abs(Eval0m) == min(abs(Eval0m))); %select eigenvalue lambda=0
0152 %
0153 Evec0p = Evecp(:,max(jp)); %eigenvector
0154 Evec0m = Evecm(:,max(jm)); %eigenvector
0155 %
0156 %   Polarization vectors
0157 %
0158 ep_xyz = Evec0p./norm(Evec0p); %xyz
0159 em_xyz = Evec0m./norm(Evec0m);
0160 %
0161 M_xypm = [1,i,0;1,-i,0;0,0,sqrt(2)]/sqrt(2);
0162 %
0163 ep_pmz = M_xypm*ep_xyz; %pmz
0164 em_pmz = M_xypm*em_xyz;
0165 %
0166 cp = Npar/norm(Nxyzp);
0167 sp = Nperpp/norm(Nxyzp); 
0168 cm = Npar/norm(Nxyzm);
0169 sm = Nperpm/norm(Nxyzm); 
0170 %
0171 Mp_xzpk = [cp,0,-sp;0,1,0;sp,0,cp];
0172 Mm_xzpk = [cm,0,-sm;0,1,0;sm,0,cm];
0173 %
0174 ep_pyk = Mp_xzpk*ep_xyz; %pyk
0175 em_pyk = Mm_xzpk*em_xyz;
0176 %
0177 %   Power flow phi = Re[e x (N x e*)] = Re[N - (N . e)e*]
0178 %
0179 Npmzp = M_xypm*Nxyzp;
0180 Npmzm = M_xypm*Nxyzm;
0181 Npykp = Mp_xzpk*Nxyzp;
0182 Npykm = Mm_xzpk*Nxyzm;
0183 %
0184 phip_xyz = real(Nxyzp - (Nxyzp.'*ep_xyz)*conj(ep_xyz));
0185 phim_xyz = real(Nxyzm - (Nxyzm.'*em_xyz)*conj(em_xyz)); 
0186 phip_pmz = real(Npmzp - (Npmzp.'*ep_pmz)*conj(ep_pmz));
0187 phim_pmz = real(Npmzm - (Npmzm.'*em_pmz)*conj(em_pmz)); 
0188 phip_pyk = real(Npykp - (Npykp.'*ep_pyk)*conj(ep_pyk)); 
0189 phim_pyk = real(Npykm - (Npykm.'*em_pyk)*conj(em_pyk));
0190 %
0191 %   Energy density
0192 %
0193 wdKperpdw = 2*sum(ps.^2./(1 - ys.^2).^2);
0194 wdKpardw = 2*sum(ps.^2);
0195 wdKxydw = -sum((3 - ys.^2).*ps.^2.*ys./(1 - ys.^2).^2);
0196 %
0197 wdKxyzdw = [wdKperpdw,-i*wdKxydw,0;i*wdKxydw,wdKperpdw,0;0,0,wdKpardw];
0198 %
0199 sigmap = ep_xyz'*(wdKxyzdw - 2*real(Nxyzp*Nxyzp.' - sum(Nxyzp.^2)*diag([1,1,1])))*ep_xyz/2;
0200 sigmam = em_xyz'*(wdKxyzdw - 2*real(Nxyzm*Nxyzm.' - sum(Nxyzm.^2)*diag([1,1,1])))*em_xyz/2;
0201 %