f3tempmodel

R5-X2 - Calculation of the electron momentum distribution function using the 3 temperatures model

 Calculation of the electron momentum distribution function using the 3
 temperatures model. Simple modeling of the fast electron tail driven by
 the LH wave, useful for benchmarking the Fast ELectron Bremsstrahlung
 calculations.

INPUTS:

   - equilDKE: magnetic equilibrium from LUKE [structure]
   - momentumDKE: momentum space grid definitions from LUKE [structure]
   - lambda: fraction of fast electrons [1,1]
   - pmin : boundary between Maxwellian distribution and tail distribution [nr,1]
   - Tfor: Forward temperature (keV) [1,1] 
   - Tperp: Perpendicular temperature (keV) [1,1]
   - Tback: Backward temperature (keV) [1,1]
   - Ecut: Up limit of the distribution (keV) [1,1]
   - mode: test mode activated if 1 [1,1]

OUTPUT:

   - XXf: distribution function normalized to unity [npn,nmhu,nr]
   - xnorm: normalisation of the distribution function (must be unity) [nr,1]

by Y. Peysson (CEA/IRFM) (yves.peysson@cea.fr) and J. Decker (CEA/IRFM) (joan.decker@cea.fr)
 Modified by F.-X Duthoit (CEA/IRFM) (francois-xavier.duthoit@cea.fr)

0001  function [XXf,xnorm] = f3tempmodel(equilDKE,momentumDKE,lambda,pmin,Tfor,Tperp,Tback,Ecut,mode)
0002 %R5-X2 - Calculation of the electron momentum distribution function using the 3 temperatures model
0003 %
0004 % Calculation of the electron momentum distribution function using the 3
0005 % temperatures model. Simple modeling of the fast electron tail driven by
0006 % the LH wave, useful for benchmarking the Fast ELectron Bremsstrahlung
0007 % calculations.
0008 %
0009 %INPUTS:
0010 %
0011 %   - equilDKE: magnetic equilibrium from LUKE [structure]
0012 %   - momentumDKE: momentum space grid definitions from LUKE [structure]
0013 %   - lambda: fraction of fast electrons [1,1]
0014 %   - pmin : boundary between Maxwellian distribution and tail distribution [nr,1]
0015 %   - Tfor: Forward temperature (keV) [1,1]
0016 %   - Tperp: Perpendicular temperature (keV) [1,1]
0017 %   - Tback: Backward temperature (keV) [1,1]
0018 %   - Ecut: Up limit of the distribution (keV) [1,1]
0019 %   - mode: test mode activated if 1 [1,1]
0020 %
0021 %OUTPUT:
0022 %
0023 %   - XXf: distribution function normalized to unity [npn,nmhu,nr]
0024 %   - xnorm: normalisation of the distribution function (must be unity) [nr,1]
0025 %
0026 %by Y. Peysson (CEA/IRFM) (yves.peysson@cea.fr) and J. Decker (CEA/IRFM) (joan.decker@cea.fr)
0027 % Modified by F.-X Duthoit (CEA/IRFM) (francois-xavier.duthoit@cea.fr)
0028 %
0029 if nargin < 7
0030     error('Wrong number of input arguments in f3tempmodel');
0031 end
0032 %
0033 if nargin == 7,
0034     xnorm = [];
0035     mode = 0;
0036 end
0037 %
0038 Xpn = momentumDKE.pn'*ones(1,length(momentumDKE.mhu));
0039 Xmhu = ones(length(momentumDKE.pn),1)*momentumDKE.mhu;
0040 Xpperp2 = Xpn.*Xpn.*(1 - Xmhu.*Xmhu);
0041 Xppar = Xpn.*Xmhu;
0042 Xppar2 = Xppar.*Xppar;
0043 %
0044 Xmask = Xppar >= 0;
0045 %
0046 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;%Universal physics constants
0047 %
0048 % Normalisation reference values
0049 %
0050 T_ref = max(equilDKE.xTe);
0051 beta_ref = sqrt(T_ref/mc2);
0052 %
0053 % Temperature and energy normalisation
0054 %
0055 Ten = equilDKE.xTe/T_ref;
0056 Tforn = repmat(Tfor,length(equilDKE.xTe),1)/T_ref;
0057 Tperpn = repmat(Tperp,length(equilDKE.xTe),1)/T_ref;
0058 Tbackn = repmat(Tback,length(equilDKE.xTe),1)/T_ref;
0059 Ecutn = repmat(Ecut,length(equilDKE.xTe),1)/mc2;
0060 %
0061 pmax = 1/beta_ref*sqrt((Ecutn+1).^2-1);
0062 %
0063 ZTen = shiftdim(repmat(Ten',[1,size(Xppar,1),size(Xppar,2)]),1);
0064 ZTforn = shiftdim(repmat(Tforn,[1,size(Xppar,1),size(Xppar,2)]),1);
0065 ZTperpn = shiftdim(repmat(Tperpn,[1,size(Xppar,1),size(Xppar,2),1]),1);
0066 ZTbackn = shiftdim(repmat(Tbackn,[1,size(Xppar,1),size(Xppar,2),1]),1);
0067 ZEcutn = shiftdim(repmat(Ecutn,[1,size(Xppar,1),size(Xppar,2)]),1);
0068 %
0069 Zpn = repmat(Xpn,[1,1,length(equilDKE.xTe)]);
0070 Zpperp2 = repmat(Xpperp2,[1,1,length(equilDKE.xTe)]);
0071 Zppar2 = repmat(Xppar2,[1,1,length(equilDKE.xTe)]);
0072 %
0073 Zpmax = shiftdim(repmat(pmax,[1,size(Xppar,1),size(Xppar,2)]),1);
0074 Zpmin = shiftdim(repmat(pmin,[1,size(Xppar,1),size(Xppar,2)]),1);
0075 %
0076 fMmask = (Zpn - Zpmin < 0).*(Zpn >= 0);
0077 f3Tmask = (Zpn - Zpmax <= 0).*(Zpn - Zpmin >= 0);
0078 Zmask = repmat(Xmask,[1,1,length(equilDKE.xTe)]);
0079 %
0080 Zgamma = sqrt(Zpn.*Zpn*beta_ref^2 + 1);
0081 %
0082 % Maxwellian distribution
0083 %
0084 Zm = exp(-Zpn.*Zpn./((Zgamma + 1).*ZTen));
0085 %
0086 % 3T distribution
0087 %
0088 Zf_for = Zmask.*exp(-Zppar2./ZTforn/2).*exp(-Zpperp2./ZTperpn/2);
0089 Zf_back = (1-Zmask).*exp(-Zppar2./ZTbackn/2).*exp(-Zpperp2./ZTperpn/2);
0090 Zf = Zf_for + Zf_back;
0091 %
0092 % Normalisation coefficients
0093 %
0094 alphaM = beta_ref^3*(4*pi*exp((Ten*beta_ref^2).^(-1))*4/3*gamma(3/2+1)/sqrt(pi).*besselk(2,(Ten*beta_ref^2).^(-1)).*(Ten*beta_ref^2)).^(-1);
0095 alpha3T = (sqrt(2*pi^3)*Tperpn.*(sqrt(Tforn)+sqrt(Tbackn))).^(-1);
0096 ZalphaM = shiftdim(repmat(alphaM',[1,size(Xppar,1),size(Xppar,2)]),1);
0097 Zalpha3T = shiftdim(repmat(alpha3T,[1,size(Xppar,1),size(Xppar,2)]),1);
0098 %
0099 % Step vectors
0100 %
0101 dx = repmat(momentumDKE.dpn',[1,length(equilDKE.xTe)]);
0102 dxmhu = momentumDKE.dmhu;
0103 %
0104 % Norm calculations
0105 %
0106 fMnorm = alphaM'.*fnorm(fMmask.*Zm);
0107 f3Tnorm = alpha3T.*fnorm(f3Tmask.*Zf);
0108 %
0109 % Distribution weight coefficients
0110 %
0111 lambdaM = (1 - lambda)*fMnorm.^(-1);
0112 lambda3T = lambda*f3Tnorm.^(-1);
0113 ZlambdaM = shiftdim(repmat(lambdaM,[1,size(Xppar,1),size(Xppar,2)]),1);
0114 Zlambda3T = shiftdim(repmat(lambda3T,[1,size(Xppar,1),size(Xppar,2)]),1);
0115 %
0116 % Distribution function
0117 %
0118 XXf = ZlambdaM.*ZalphaM.*Zm.*fMmask + Zlambda3T.*Zalpha3T.*Zf.*f3Tmask;
0119 xnorm = fnorm(XXf);
0120 %
0121 % Norm function
0122 %
0123      function y = fnorm(x)
0124          y = 2*pi*integral_dke_jd(ones(size(Xppar,1),1),dx.*integral_dke_jd(dxmhu,Zpn.*Zpn.*x,2),1);
0125      end
0126  end