hottensor_jd

 This function calculates the kinetic plasma susceptibility tensor. 

 - (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: perpendicular index of refraction [1,1] 
       - herm_mode: mode for the calculation: (0) full susceptibility tensor 
                                              (1) hermitian part 
           >>> default: 0
       - nmax: maximum harmonic number
           >>> default: 10
       - NZ: number of grid points in Z function integration calculation 
           >>> default: 1000

   OUTPUTS:

       - X: kinetic plasma susceptibility tensor elements [3,3]
       - dXdNpar: derivative of X with respect to Npar
       - dXdNperp: derivative of X with respect to Nperp

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

0001 function [X,dXdNpar,dXdNperp] = hottensor_jd(sa,sy,sbeta,Npar,Nperp,herm_mode,nmax,NZ)
0002 %
0003 % This function calculates the kinetic plasma susceptibility tensor.
0004 %
0005 % - (x,y,z): We assume that the magnetic field is in the z direction and
0006 %         the perpendicular index of refraction in the x direction.
0007 % - (p,m,z): p and m refer to the left hand and right hand rotating fields
0008 %         respectively.
0009 % - (p,y,k): p and y are the transverse components of the polarization,
0010 %         with respect to the wave vector. k is the longitudinal component.
0011 %
0012 %   INPUTS:
0013 %
0014 %       - sa: ratio (wps/w)^2 for each species [Ns,1]
0015 %       - sy: ratio wcs/w for each species [Ns,1]
0016 %           >>> Warning: sy < 0 for negatively charged species
0017 %       - sbeta: ratio (Ts/mc2)^(1/2) for each species [Ns,1]
0018 %       - Npar: parallel index of refraction [1,1]
0019 %       - Nperp: perpendicular index of refraction [1,1]
0020 %       - herm_mode: mode for the calculation: (0) full susceptibility tensor
0021 %                                              (1) hermitian part
0022 %           >>> default: 0
0023 %       - nmax: maximum harmonic number
0024 %           >>> default: 10
0025 %       - NZ: number of grid points in Z function integration calculation
0026 %           >>> default: 1000
0027 %
0028 %   OUTPUTS:
0029 %
0030 %       - X: kinetic plasma susceptibility tensor elements [3,3]
0031 %       - dXdNpar: derivative of X with respect to Npar
0032 %       - dXdNperp: derivative of X with respect to Nperp
0033 %
0034 % created 05/02/2003
0035 % by Joan Decker <joan.decker@cea.fr> and
0036 %    Yves Peysson <yves.peysson@cea.fr>
0037 %
0038 if nargin < 8,
0039    NZ = 1000;
0040 end
0041 if nargin < 7,
0042    nmax = 10;
0043 end
0044 if nargin < 6,
0045    herm_mode = 0;
0046 end
0047 if nargin < 5,
0048    error('Not enough arguments');
0049 end
0050 %
0051 % Initializations
0052 %
0053 sa = sa(:);
0054 sy = sy(:);
0055 sbeta = sbeta(:);
0056 %
0057 Ns = length(sa);
0058 %
0059 if length(sy) ~= Ns || length(sbeta) ~= Ns,
0060     error('Argument size not compatible')
0061 end
0062 %
0063 sy(sy.^2 == 1) = sy(sy.^2 == 1) + eps;%to avoid singularity
0064 %
0065 ss = sign(sy);
0066 spar = sign(Npar);
0067 %
0068 n = (-nmax:nmax);
0069 Nn = 2*nmax+1;
0070 n0 = nmax+1;
0071 %
0072 % Calculation of the plasma dispersion function
0073 %
0074 zeta_sn = (1 - abs(sy)*n)./(sqrt(2)*abs(Npar)*(sbeta*ones(1,Nn)));
0075 %
0076 if exist(['erfcc.',mexext],'file') == 3,%fortran file
0077     Z_sn = zeros(Ns,Nn);
0078     for is = 1:Ns,
0079         %
0080         x = [zeta_sn(is,:),zeros(1,1024 - Nn)];
0081         %
0082         [yr,yi] = erfcc(Nn,real(x),imag(x));
0083         %
0084         Z_sn(is,:) = sqrt(pi)*(i*yr(1:Nn) - yi(1:Nn));
0085     end
0086     Zp_sn = -2*(1 + zeta_sn.*Z_sn);
0087     Zpp_sn = -2*(Z_sn - 2*zeta_sn.*(1 + zeta_sn.*Z_sn));
0088     %
0089     maskp_sn = abs(zeta_sn) > sqrt(1e-2/(2*eps));
0090     maskpp_sn = abs(zeta_sn) > sqrt(sqrt(1e-2/(4*eps)));
0091     %
0092     if any(maskp_sn),
0093         disp(['Warning : there are ',num2str(sum(maskp_sn(:))),'/',num2str(length(maskp_sn(:))),...
0094             ' points outside the plasma dispersion function first derivative calculation range'])
0095     end
0096     %
0097     if any(maskpp_sn),
0098         disp(['Warning : there are ',num2str(sum(maskpp_sn(:))),'/',num2str(length(maskpp_sn(:))),...
0099             ' points outside the plasma dispersion function second derivative calculation range'])
0100     end
0101     %
0102 else
0103     %
0104     error('Please compile the function erfcc.f as a mex file.');
0105     %
0106     % Warning: Zplasma must be checked first. Second derivative
0107     % calculations are wrong for imag. part and for small arguments.
0108     %
0109     %[Z_sn,Zp_sn,Zpp_sn] = Zplasma(zeta_sn,NZ);
0110     %
0111 end
0112 %
0113 %   Dielectric tensor elements
0114 %
0115 Y_sn = zeros(3,3,Ns,Nn);
0116 %
0117 Msa_sn = shiftdim(repmat(sa,[1,Nn,3,3]),2);
0118 Mss_sn = reshape(repmat(ss,[1,Nn]),[1,1,Ns,Nn]);
0119 Mn_sn = reshape(repmat(n,[Ns,1]),[1,1,Ns,Nn]);
0120 %
0121 Mzeta_sn = reshape(zeta_sn,[1,1,Ns,Nn]);
0122 Mzeta0_sn = repmat(Mzeta_sn(:,:,:,n0),[3,3,1,Nn]);
0123 MZ_sn = reshape(Z_sn,[1,1,Ns,Nn]);
0124 MZp_sn = reshape(Zp_sn,[1,1,Ns,Nn]);
0125 %
0126 slambda = (Nperp*sbeta./sy).^2;
0127 Gamma_sn = besseli(n,slambda,1).*(exp(-i*imag(slambda))*ones(1,Nn));
0128 Gammap_sn = ((besseli(n+1,slambda,1) + besseli(n-1,slambda,1))/2 - besseli(n,slambda,1)).*(exp(-i*imag(slambda))*ones(1,Nn));
0129 %
0130 MGamma_sn = reshape(Gamma_sn,[1,1,Ns,Nn]);
0131 MGammap_sn = reshape(Gammap_sn,[1,1,Ns,Nn]);
0132 %
0133 Mlambda_sn = reshape(repmat(slambda,[1,Nn]),[1,1,Ns,Nn]);
0134 %
0135 Y_sn(1,1,:,:) = Mn_sn.^2.*MGamma_sn.*MZ_sn./Mlambda_sn;%Yxx
0136 Y_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*MGammap_sn.*MZ_sn;%Yxy
0137 Y_sn(1,3,:,:) = -spar.*Mn_sn.*MGamma_sn.*MZp_sn./sqrt(2*Mlambda_sn);%Yxz
0138 Y_sn(2,1,:,:) = -Y_sn(1,2,:,:);%Yyx
0139 Y_sn(2,2,:,:) = (Mn_sn.^2.*MGamma_sn./Mlambda_sn - 2*Mlambda_sn.*MGammap_sn).*MZ_sn;%Yyy
0140 Y_sn(2,3,:,:) = i*Mss_sn*spar.*sqrt(Mlambda_sn/2).*MGammap_sn.*MZp_sn;%Yyz
0141 Y_sn(3,1,:,:) = Y_sn(1,3,:,:);%Yzx
0142 Y_sn(3,2,:,:) = -Y_sn(2,3,:,:);%Yzy
0143 Y_sn(3,3,:,:) = -MGamma_sn.*Mzeta_sn.*MZp_sn;%Yzz
0144 %
0145 X_sn = Msa_sn.*Mzeta0_sn.*Y_sn;%Xsn
0146 %
0147 X_s = sum(X_sn,4);%Xs
0148 %
0149 X = sum(X_s,3);
0150 %
0151 if herm_mode == 1,
0152     %
0153     X = (X + X')/2;
0154     %
0155     dYdNperp_sn = zeros(3,3,Ns,Nn);
0156     %
0157     Mlambdap_sn = 2*Mlambda_sn/Nperp;
0158     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));
0159     MGammapp_sn = reshape(Gammapp_sn,[1,1,Ns,Nn]);
0160     %
0161     dYdNperp_sn(1,1,:,:) = Mn_sn.^2.*Mlambdap_sn.*(MGammap_sn - MGamma_sn./Mlambda_sn).*MZ_sn./Mlambda_sn;
0162     dYdNperp_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*Mlambdap_sn.*MGammapp_sn.*MZ_sn;
0163     dYdNperp_sn(1,3,:,:) = -spar.*Mn_sn.*Mlambdap_sn.*(MGammap_sn - MGamma_sn/2./Mlambda_sn).*MZp_sn./sqrt(2*Mlambda_sn);
0164     dYdNperp_sn(2,1,:,:) = -dYdNperp_sn(1,2,:,:);
0165     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;
0166     dYdNperp_sn(2,3,:,:) = i*Mss_sn*spar.*Mlambdap_sn.*sqrt(Mlambda_sn/2).*(MGammap_sn/2./Mlambda_sn + MGammapp_sn).*MZp_sn;
0167     dYdNperp_sn(3,1,:,:) = dYdNperp_sn(1,3,:,:);
0168     dYdNperp_sn(3,2,:,:) = -dYdNperp_sn(2,3,:,:);
0169     dYdNperp_sn(3,3,:,:) = -Mlambdap_sn.*MGammap_sn.*Mzeta_sn.*MZp_sn;
0170     %
0171     dXdNperp_sn = Msa_sn.*Mzeta0_sn.*dYdNperp_sn;
0172     %
0173     dXdNperp_s = sum(dXdNperp_sn,4);
0174     %
0175     dXdNperp = sum(dXdNperp_s,3);
0176     %
0177     dXdNperp = (dXdNperp + dXdNperp')/2;%hermitian part only
0178     %
0179     dYdNpar_sn = zeros(3,3,Ns,Nn);
0180     %
0181     Mzetap_sn = -Mzeta_sn/Npar;
0182     Mzetap0_sn = -Mzeta0_sn/Npar;
0183     MZpp_sn = reshape(Zpp_sn,[1,1,Ns,Nn]);
0184     %
0185     dYdNpar_sn(1,1,:,:) = Mn_sn.^2.*MGamma_sn.*Mzetap_sn.*MZp_sn./Mlambda_sn;
0186     dYdNpar_sn(1,2,:,:) = i*Mss_sn.*Mn_sn.*MGammap_sn.*Mzetap_sn.*MZp_sn;
0187     dYdNpar_sn(1,3,:,:) = -spar.*Mn_sn.*MGamma_sn.*Mzetap_sn.*MZpp_sn./sqrt(2*Mlambda_sn);
0188     dYdNpar_sn(2,1,:,:) = -dYdNpar_sn(1,2,:,:);
0189     dYdNpar_sn(2,2,:,:) = (Mn_sn.^2.*MGamma_sn./Mlambda_sn - 2*Mlambda_sn.*MGammap_sn).*Mzetap_sn.*MZp_sn;
0190     dYdNpar_sn(2,3,:,:) = i*Mss_sn*spar.*sqrt(Mlambda_sn/2).*MGammap_sn.*Mzetap_sn.*MZpp_sn;
0191     dYdNpar_sn(3,1,:,:) = dYdNpar_sn(1,3,:,:);
0192     dYdNpar_sn(3,2,:,:) = -dYdNpar_sn(2,3,:,:);
0193     dYdNpar_sn(3,3,:,:) = -MGamma_sn.*Mzetap_sn.*(MZp_sn + Mzeta_sn.*MZpp_sn);
0194     %
0195     dXdNpar_sn = Msa_sn.*Mzetap0_sn.*Y_sn + Msa_sn.*Mzeta0_sn.*dYdNpar_sn;
0196     %
0197     dXdNpar_s = sum(dXdNpar_sn,4);
0198     %
0199     dXdNpar = sum(dXdNpar_s,3);
0200     %
0201     dXdNpar  = (dXdNpar + dXdNpar')/2;
0202     %
0203 end