LUKE - calculate the electron distribution function at the banana tip for magnetic ripple losses
Calculation the electron distribution function at the banana tip for magnetic ripple losses.
Upper or lower part of the plasma is chosen depending upon the vertical drift direction of the electrons.
INPUTS:
- ntheta: number of poloidal positions where the banan tips are calculated
- psin: normalized magnetic flux surfaces (given by LUKE) [1,npsi]
- pn: normalized momentum grid (given by LUKE)
- mhu: pitch-angle grid (given by LUKE)
- equil: toroidal MHD equilibrium structure (given by LUKE)
- XXf0 : 3-D electron distribution function solution of the Fokker-Planck equation (given by LUKE)
- mksa : MKSA data structure (given by LUKE)
OUTPUTS:
- ripple: output structure
ripple.psin_tip: normalized magnetic flux surface label where the density banana tips are calculated (0 at the plasma center, 1 at the separatrix)
ripple.rho_tip: normalized minor radius where the density banana tips are calculated (0 at the plasma center, 1 at the separatrix)
ripple.theta_tip: poloidal angle where the density banana tips are calculated
ripple.Ec_tip: kinetic energies at the banana tips [keV]
ripple.pn_tip: normalized momentum values at the banana tip (p/pth)
ripple.mhu_tip: pitch-angle at Bmin corresponding to the banana tip
ripple.XXf0_tip: Normalized electron distribution function of banana tips to ne_ref
ripple.R_tip: Major radius where the density banana tips are calculated
ripple.Z_tip: Vertical position where the density banana tips are calculated
ripple.ne_ref: reference density for XXf0_tip (m-3)
By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)0001 function ripple = bananatip_yp(ntheta,psin,dpsin,pn,mhu,equil,XXf0,pn2dpn,mksa) 0002 %LUKE - calculate the electron distribution function at the banana tip for magnetic ripple losses 0003 % 0004 % Calculation the electron distribution function at the banana tip for magnetic ripple losses. 0005 % Upper or lower part of the plasma is chosen depending upon the vertical drift direction of the electrons. 0006 % 0007 % INPUTS: 0008 % 0009 % - ntheta: number of poloidal positions where the banan tips are calculated 0010 % - psin: normalized magnetic flux surfaces (given by LUKE) [1,npsi] 0011 % - pn: normalized momentum grid (given by LUKE) 0012 % - mhu: pitch-angle grid (given by LUKE) 0013 % - equil: toroidal MHD equilibrium structure (given by LUKE) 0014 % - XXf0 : 3-D electron distribution function solution of the Fokker-Planck equation (given by LUKE) 0015 % - mksa : MKSA data structure (given by LUKE) 0016 % 0017 % OUTPUTS: 0018 % 0019 % - ripple: output structure 0020 % ripple.psin_tip: normalized magnetic flux surface label where the density banana tips are calculated (0 at the plasma center, 1 at the separatrix) 0021 % ripple.rho_tip: normalized minor radius where the density banana tips are calculated (0 at the plasma center, 1 at the separatrix) 0022 % ripple.theta_tip: poloidal angle where the density banana tips are calculated 0023 % ripple.Ec_tip: kinetic energies at the banana tips [keV] 0024 % ripple.pn_tip: normalized momentum values at the banana tip (p/pth) 0025 % ripple.mhu_tip: pitch-angle at Bmin corresponding to the banana tip 0026 % ripple.XXf0_tip: Normalized electron distribution function of banana tips to ne_ref 0027 % ripple.R_tip: Major radius where the density banana tips are calculated 0028 % ripple.Z_tip: Vertical position where the density banana tips are calculated 0029 % ripple.ne_ref: reference density for XXf0_tip (m-3) 0030 % 0031 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr) 0032 % 0033 if nargin < 7,error('Error: wrong number of input arguments in bananatip_yp.m');end 0034 % 0035 if isfield(equil,'equil_fit'), 0036 equil_fit = equil.equil_fit; 0037 else 0038 fitparam.mode_equil = 1;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y) 0039 fitparam.nharm = NaN;%Number of harmonics in the magnetic equilibrium interpolation (NaN, Inf or empty, nharm = ntheta-1) 0040 fitparam.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters 0041 fitparam.method = 'spline';%nearest,spline,pchip 0042 % 0043 equil_fit = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);%Build vectorized magnetic equilibrium structure 0044 % 0045 disp('----> Vectorization of the toroidal MHD equilibrium done.'); 0046 end 0047 % 0048 method = 'spline'; 0049 n_gauleg = 35;%Number of Legendre polynomials for projections of the electron velocity distribution 0050 n_pol = n_gauleg;%Necessary for accurate reconstruction at ksi = -1 0051 epsi_gauleg = 1e-14;%Accuracy of the Legendre projection 0052 % 0053 [x_leg,w_leg] = gauleg_dke_yp(n_gauleg,epsi_gauleg);%Gauss-Legendre abscissa and weigth calculations between -1 and + 1 0054 [Pm,cL] = leg_dke_yp(x_leg(:)',n_pol);%Legendre polynomial calculation 0055 % 0056 XXf0_c = zeros(size(XXf0,1),n_gauleg+1,size(XXf0,3)); 0057 % 0058 for ir = 1:size(XXf0,3),%radial position 0059 Xf0 = XXf0(:,:,ir);%Non-thermal distribution function at poloidal location 0060 Xf0_i = interp1(mhu,Xf0',x_leg(:),method)';%Interpolation of the distribution function on the Legendre polynomial pitch-angle grid 0061 Xf0_c = (Xf0_i.*(w_leg*ones(1,size(Xf0_i,1)))')*Pm';%Calculation of the Legendre coefficients for local values of the distribution function 0062 %Xf0_i1 = (Pm'*(Xf0_c.*(ones(size(Xf0_i,1),1)*cL'))')';%Distribution function on the Legendre polynomials pitch-angle grid back calculated 0063 %epsif0 = (Xf0_i-Xf0_i1)./Xf0_i;epsif0(isnan(epsif0)==1)=0; 0064 XXf0_c(:,:,ir) = Xf0_c; 0065 end 0066 % 0067 disp('----> Projection of the electron distribution function at Bmin on the Legendre polynomial basis done.'); 0068 % 0069 [~,~,~,~,~,~,~,mc2,~,~] = pc_dke_yp;%Universal physics constants 0070 % 0071 theta = linspace(0,2*pi,ntheta);%uniform poloidal angle grid grid 0072 dtheta = theta(2) - theta(1); 0073 % 0074 ripple.psin_tip = psin;%Fokker-Planck psi grid 0075 % 0076 ripple.rho_tip = psi2rho_jd(equil,ripple.psin_tip);%Fokker-Planck rho grid 0077 ripple.theta_tip = theta;%poloidal angle value at which the banana tip is located 0078 % 0079 ripple.Ec_tip = mc2*(sqrt(1 + pn.*pn*mksa.betath_ref^2) - 1);%Fast electron kinetic energy in keV at the banana tip; 0080 ripple.pn_tip = pn; 0081 % 0082 ripple.mhu_tip = zeros(size(XXf0,3),ntheta); 0083 % 0084 ripple.XXf0_tip = zeros(size(XXf0,1),ntheta,size(XXf0,3));%at the banana tip 0085 % 0086 Xpn2dpn = pn2dpn(:)*ones(1,ntheta); 0087 % 0088 for ir = 1:length(ripple.rho_tip), 0089 % 0090 x_a0 = ppval(equil_fit.x_fit.pp_a0,ripple.rho_tip(ir)); 0091 x_an = ppval(equil_fit.x_fit.pp_an,ripple.rho_tip(ir)); 0092 x_bn = ppval(equil_fit.x_fit.pp_bn,ripple.rho_tip(ir)); 0093 % 0094 y_a0 = ppval(equil_fit.y_fit.pp_a0,ripple.rho_tip(ir)); 0095 y_an = ppval(equil_fit.y_fit.pp_an,ripple.rho_tip(ir)); 0096 y_bn = ppval(equil_fit.y_fit.pp_bn,ripple.rho_tip(ir)); 0097 % 0098 ripple.R_tip(:,ir) = calcval_yp(equil_fit,ripple.theta_tip,x_a0,x_an,x_bn)' + equil.Rp; 0099 ripple.Z_tip(:,ir) = calcval_yp(equil_fit,ripple.theta_tip,y_a0,y_an,y_bn)' + equil.Zp; 0100 % 0101 B_a0 = ppval(equil_fit.B_fit.pp_a0,ripple.rho_tip(ir)); 0102 B_an = ppval(equil_fit.B_fit.pp_an,ripple.rho_tip(ir)); 0103 B_bn = ppval(equil_fit.B_fit.pp_bn,ripple.rho_tip(ir)); 0104 % 0105 BP_a0 = ppval(equil_fit.BP_fit.pp_a0,ripple.rho_tip(ir)); 0106 BP_an = ppval(equil_fit.BP_fit.pp_an,ripple.rho_tip(ir)); 0107 BP_bn = ppval(equil_fit.BP_fit.pp_bn,ripple.rho_tip(ir)); 0108 % 0109 r_a0 = ppval(equil_fit.r_fit.pp_a0,ripple.rho_tip(ir)); 0110 r_an = ppval(equil_fit.r_fit.pp_an,ripple.rho_tip(ir)); 0111 r_bn = ppval(equil_fit.r_fit.pp_bn,ripple.rho_tip(ir)); 0112 % 0113 calpha_a0 = ppval(equil_fit.calpha_fit.pp_a0,ripple.rho_tip(ir)); 0114 calpha_an = ppval(equil_fit.calpha_fit.pp_an,ripple.rho_tip(ir)); 0115 calpha_bn = ppval(equil_fit.calpha_fit.pp_bn,ripple.rho_tip(ir)); 0116 % 0117 B_tip = calcval_yp(equil_fit,ripple.theta_tip,B_a0,B_an,B_bn); 0118 Bmin = calcval_yp(equil_fit,0,B_a0,B_an,B_bn); 0119 BP_tip = calcval_yp(equil_fit,ripple.theta_tip,BP_a0,BP_an,BP_bn); 0120 r_tip = calcval_yp(equil_fit,ripple.theta_tip,r_a0,r_an,r_bn); 0121 calpha_tip = calcval_yp(equil_fit,ripple.theta_tip,calpha_a0,calpha_an,calpha_bn); 0122 % 0123 ripple.mhu_tip(ir,:) = sqrt(1 - Bmin./B_tip);%pitch-angle value at Bmin which corresponds to the poloidal angle value at which the banana tip is located 0124 % 0125 [Pm,cL] = leg_dke_yp(ripple.mhu_tip(ir,:),n_pol);%Value of the Legendre polynomial on the mhu_tip grid 0126 % 0127 ripple.XXf0_tip(:,:,ir) = (Pm'*(XXf0_c(:,:,ir).*(ones(size(XXf0,1),1)*cL'))')';%the distribution of electron banana tips on the grid 0128 % 0129 XJ_tip = ones(size(Xpn2dpn,1),1)*abs(r_tip./BP_tip./calpha_tip)'; 0130 % 0131 ripple.DNeDphi_tip(:,:,ir) = 4*pi*mksa.ne_ref*XJ_tip.*ripple.XXf0_tip(:,:,ir).*Xpn2dpn*dpsin(ir)*dtheta;%the incremental density of electron banana tips on the grid 0132 % 0133 end 0134 % 0135 ripple.ne_ref = mksa.ne_ref; 0136 % 0137