bounce_dke_jd

 calculates bounce averaging coefficients needed for solving the electron drift-kinetic-equation. 

   INPUTS:

       - ymhu0: pitch angle grid [1,nmhu]
       - plasma magnetic equilibrium

   OUTPUTS:
   

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

0001 function [xymaskT,xylambda_p0p0,xylambda_p2m1,xylambda_p0p3,xylambda_p2p2,xylambda_p2m2p2,xylambda_b_p1p1,xylambda_b_p1m1p2,xylambda_b_p1m1,xylambda_b_p1m2,xylambda_b_p2m2p2,xylambda_b_p3m2,xylambda_b_p2,xylambda_b_p1m3p4,xylambda_b_p1m2p2]...
0002 = bounce_dke_jd(equil_mode,ymhu0,equilDKE,xqtilde,xymaskT_in,sign)
0003 %
0004 % calculates bounce averaging coefficients needed for solving the electron drift-kinetic-equation.
0005 %
0006 %   INPUTS:
0007 %
0008 %       - ymhu0: pitch angle grid [1,nmhu]
0009 %       - plasma magnetic equilibrium
0010 %
0011 %   OUTPUTS:
0012 %
0013 %
0014 % by Joan Decker <jodecker@mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0015 %
0016 if nargin < 5,
0017     sign = 0;
0018     xymaskT_in = NaN;
0019 end
0020 if nargin == 5,
0021     error('Inconsistent number of arguments');
0022 end
0023 if nargin < 4,
0024     error('Not enough arguments');
0025 end
0026 %
0027 % Case of a numerical equilibrium
0028 %
0029 nmhu = length(ymhu0);
0030 nr = size(equilDKE.Xx,1);
0031 %
0032 if equil_mode == 1,
0033     [YMdt,YMmask,YMF,YMPSIB,YMmhu0,YMmhu,YMx,xymaskT,nm2] = banana_dke_jd(ymhu0,equilDKE.ap,equilDKE.Xx,equilDKE.Xy,equilDKE.XBx,equilDKE.XBy,equilDKE.XBphi);    
0034 elseif equil_mode == 0,
0035     for ir = 1:nr   
0036         xymaskT(ir,:) = ymhu0.^2 < equilDKE.xmhu0T2(ir);        
0037     end
0038 end            
0039 %
0040 if ~isnan(xymaskT_in),%for correct grid mask position
0041     xymaskT = xymaskT_in;
0042     for ir = 1:nr,
0043         xjmhut_min(ir) = min(find(xymaskT_in(ir,:) == 1));
0044         xjmhut_max(ir) = max(find(xymaskT_in(ir,:) == 1));
0045     end
0046     if sign == 1, xymaskT(:,max([xjmhut_min-1;ones(1,nr)])) = xymaskT_in(:,xjmhut_min);end
0047     if sign == -1,xymaskT(:,min([xjmhut_max+1;nmhu*ones(1,nr)])) = xymaskT_in(:,xjmhut_max);end
0048 end
0049 %
0050 if equil_mode == 1
0051 %
0052     xyqtilde = repmat(xqtilde',[1,nmhu]);
0053 %
0054     YMx0 = repmat(equilDKE.xx0,[nm2,1,nmhu]);
0055 %
0056     YMYM01 = YMmhu./YMmhu0;
0057     YMYM02 = YMYM01.*YMYM01;
0058     YMYM03 = YMYM01.*YMYM02;
0059 %
0060     YMPSIBM1 = 1.0./YMPSIB;
0061     YMPSIBM2 = YMPSIBM1.*YMPSIBM1;
0062     YMPSIBM3 = YMPSIBM1.*YMPSIBM2;
0063     YMPSIBP2 = YMPSIB.*YMPSIB;
0064     YMPSIBP3 = YMPSIBP2.*YMPSIB;
0065 %
0066     YMR0YMR1 = (1 + YMx0/equilDKE.Rp)./(1 + YMx/equilDKE.Rp);
0067     YMR0YMR2 = YMR0YMR1.*YMR0YMR1;
0068     YMR0YMR4 = YMR0YMR2.*YMR0YMR2;
0069 %
0070     YMmaskYMdtYMF = YMmask.*YMdt.*YMF;
0071     YMmaskYMdtYMF(YMmask == 0) = 0;%to deal with cases where YMF is NaN
0072     not_xymaskT = ~xymaskT;
0073 %
0074     Yweight_p0p0 = ones(size(YMYM01));%(YMmhu./YMmhu0).^(0).*YMPSIB.^(0);
0075     xylambda_p0p0 = reshape(sum(YMmaskYMdtYMF.*Yweight_p0p0)/(2*pi),[nr,nmhu])./xyqtilde;
0076 %
0077 %For flux surface averaged Ohmic electric field
0078 %
0079 
0080 %
0081     Yweight_p1m3p4 = YMYM01.*YMPSIBM3.*YMR0YMR4;
0082     xylambda_p1m3p4 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1m3p4)/(2*pi),[nr,nmhu])./xyqtilde;
0083     xylambda_b_p1m3p4 = not_xymaskT.*xylambda_p1m3p4;
0084 %
0085     Yweight_p1m2p2 = YMYM01.*YMPSIBM2.*YMR0YMR2;
0086     xylambda_p1m2p2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1m2p2)/(2*pi),[nr,nmhu])./xyqtilde;
0087     xylambda_b_p1m2p2 = not_xymaskT.*xylambda_p1m2p2;
0088 %
0089 %Fokker-Planck weights
0090 %
0091     Yweight_p2m1 = YMYM02.*YMPSIBM1;
0092     xylambda_p2m1 = reshape(sum(YMmaskYMdtYMF.*Yweight_p2m1)/(2*pi),[nr,nmhu])./xyqtilde;
0093 %
0094     Yweight_p1p1 = YMYM01.*YMPSIB;
0095     xylambda_p1p1 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1p1)/(2*pi),[nr,nmhu])./xyqtilde;
0096     xylambda_b_p1p1 = not_xymaskT.*xylambda_p1p1;
0097 %
0098     Yweight_p1m1p2 = YMYM01.*YMPSIBM1.*YMR0YMR2;
0099     xylambda_p1m1p2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1m1p2)/(2*pi),[nr,nmhu])./xyqtilde;
0100     xylambda_b_p1m1p2 = not_xymaskT.*xylambda_p1m1p2;%Equivalent form: xq.*xRmin.*xB0./xqtilde/R0./xBTmin
0101 %
0102     Yweight_p0p3 = YMPSIBP3;
0103     xylambda_p0p3 = reshape(sum(YMmaskYMdtYMF.*Yweight_p0p3)/(2*pi),[nr,nmhu])./xyqtilde;
0104 %
0105     Yweight_p2p2 = YMYM02.*YMPSIBP2;
0106     xylambda_p2p2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p2p2)/(2*pi),[nr,nmhu])./xyqtilde;
0107 %
0108 %Drift kinetic weights
0109 %
0110     Yweight_p1m1 = YMYM01.*YMPSIBM1;
0111     xylambda_p1m1 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1m1)/(2*pi),[nr,nmhu])./xyqtilde;
0112     xylambda_b_p1m1 = not_xymaskT.*xylambda_p1m1;
0113 %
0114     Yweight_p1m2 = YMYM01.*YMPSIBM2;
0115     xylambda_p1m2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p1m2)/(2*pi),[nr,nmhu])./xyqtilde;
0116     xylambda_b_p1m2 = not_xymaskT.*xylambda_p1m2;%In old notation: 1-delta/mhu0^2
0117 %
0118     Yweight_p2m2p2 = YMYM02.*YMPSIBM2.*YMR0YMR2;    
0119     xylambda_p2m2p2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p2m2p2)/(2*pi),[nr,nmhu])./xyqtilde;%************probl?me sur le courant en mode num�rique (grille uniforme ou non)
0120     xylambda_b_p2m2p2 = not_xymaskT.*xylambda_p2m2p2;
0121 %
0122     Yweight_p3m2 = YMYM03.*YMPSIBM2;
0123     xylambda_p3m2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p3m2)/(2*pi),[nr,nmhu])./xyqtilde;
0124     xylambda_b_p3m2 = not_xymaskT.*xylambda_p3m2;
0125 %
0126     Yweight_p2 = YMYM02;
0127     xylambda_p2 = reshape(sum(YMmaskYMdtYMF.*Yweight_p2)/(2*pi),[nr,nmhu])./xyqtilde;%In old notation: s*_tilde
0128     xylambda_b_p2 = not_xymaskT.*xylambda_p2;
0129 %
0130     ymask = find(ymhu0 == 0);
0131     if ~isempty(ymask),
0132         xylambda_p0p0(:,ymask) = eps;
0133         xylambda_p2m1(:,ymask) = eps;
0134         xylambda_p0p3(:,ymask) = eps;
0135         xylambda_p2p2(:,ymask) = eps;
0136         xylambda_b_p1p1(:,ymask) = eps;
0137         xylambda_b_p1m1p2(:,ymask) = eps;
0138         xylambda_b_p1m1(:,ymask) = eps;
0139         xylambda_b_p1m2(:,ymask) = eps;
0140         xylambda_b_p2m2p2(:,ymask) = eps;
0141         xylambda_b_p3m2(:,ymask) = eps;
0142         xylambda_b_p2(:,ymask) = eps;
0143     end
0144 else,%Ideal equilibrium with circular concentric flux-surfaces.
0145     m = 30;
0146     xrho = equilDKE.xx0/equilDKE.ap; %note: for the sake of bounce-averaged calculations in the
0147                           %      analytical limit of circ. conc.
0148                           %      flux-surfaces, we need theta0 = 0, thus
0149                           %      this definition of xrho. Of course, when
0150                           %      the configuration is circ. conc. ro start
0151                           %      with, the two definitions coincide.
0152     xia = xrho*equilDKE.ap/equilDKE.Rp;
0153 %
0154     ymhu0(find(ymhu0==0)) = eps;%Avoid numerical singularities
0155     ymhu2 = ymhu0.^2;
0156 %
0157     for ir = 1:nr
0158         if equilDKE.xmhu0T2(ir) == 0,
0159             xylambda_p0p0(ir,:) = zeros(size(ymhu0));
0160             xylambda_p2m1(ir,:) = zeros(size(ymhu0));
0161             xylambda_p0p3(ir,:) = zeros(size(ymhu0));
0162             xylambda_p2p2(ir,:) = zeros(size(ymhu0));
0163             xylambda_p2m2p2(ir,:) = zeros(size(ymhu0));
0164             xylambda_b_p1p1(ir,:) = zeros(size(ymhu0));
0165             xylambda_b_p1m1(ir,:) = zeros(size(ymhu0));
0166             xylambda_b_p1m2(ir,:) = zeros(size(ymhu0));
0167             xylambda_b_p3m2(ir,:) = zeros(size(ymhu0));     
0168         else
0169             YJ = jseries_dke_jd(sqrt(equilDKE.xmhu0T2(ir)),ymhu0,m);
0170 %
0171             chim = zeros(m+1,1);
0172             chim_tp = zeros(m+1,1);
0173             mhutrap2m = zeros(m+1,1);
0174             chim(1) = 1;
0175             chim_tp(1) = 1;
0176             mhutrap2m(1) = 1;
0177             for i = 1:m,
0178                 chim(i+1) = (2*i - 3)/(2*i)*chim(i);
0179                 chim_tp(i+1) = (2*i - 1)/(2*i)*chim_tp(i);
0180                 mhutrap2m(i+1) = equilDKE.xmhu0T2(ir)*mhutrap2m(i);
0181             end
0182 %
0183             Ymchim = chim*ones(1,nmhu);
0184             Ymchim_tp = chim_tp*ones(1,nmhu);
0185             Ymmhutrap2m = mhutrap2m*ones(1,nmhu);
0186             Ymmhu2 = ones(m-1,1)*ymhu2;
0187             Ymhu2 = ones(m,1)*ymhu2;
0188             Ymhu2(Ymhu2==0) = eps;
0189             Ymmhu2(Ymmhu2==0) = eps;
0190 %
0191             ysstar = 1/sqrt(1- equilDKE.xmhu0T2(ir)).*ones(1,nmhu);
0192             ylambda = 2/pi*sum(Ymchim(1:m,:).*Ymmhutrap2m(1:m,:).*YJ(1:m,:));
0193             ylambda(ylambda==0) = eps;
0194             ydeltalambda = 2/pi*sum(Ymchim(1:m,:).*Ymmhutrap2m(2:m+1,:).*YJ(2:m+1,:));
0195             ysstar_tp = 2/pi*sum((Ymchim(1:m,:) - Ymchim_tp(1:m,:).*(1-Ymhu2)).*Ymmhutrap2m(1:m,:)./Ymhu2.*YJ(1:m,:));
0196 %
0197             if m > 2,
0198                 ylambda_tp = 2/pi*(YJ(1,:) + sum((Ymchim(2:m,:) - Ymchim(1:m-1,:)./Ymmhu2).*Ymmhutrap2m(2:m,:).*YJ(2:m,:)) - Ymchim(m,:)./ymhu2.*Ymmhutrap2m(m+1,:).*YJ(m+1,:));
0199             elseif m == 2,
0200                 ylambda_tp = 2/pi*(YJ(1,:) + (Ymchim(2,:) - Ymchim(1,:)./Ymmhu2).*Ymmhutrap2m(2,:).*YJ(2,:) - Ymchim(2,:)./ymhu2.*Ymmhutrap2m(3,:).*YJ(3,:));
0201             end
0202 %
0203             xylambda_p0p0(ir,:) = ylambda;
0204             xylambda_p2m1(ir,:) = ylambda - ydeltalambda./ymhu2;
0205 %             xylambda_p0p3(ir,:) = zeros(size(ymhu0));%% TODO
0206 %             xylambda_p2p2(ir,:) = zeros(size(ymhu0));
0207             xylambda_p2m2p2(ir,:) = ysstar_tp;
0208             xylambda_b_p1p1(ir,:) = ~xymaskT(ir,:).*ysstar;
0209             xylambda_b_p1m1(ir,:) = ~xymaskT(ir,:)/(1 + xia(ir));
0210             xylambda_b_p1m2(ir,:) = ~xymaskT(ir,:)*(1 + xia(ir)^2/2)/(1 + xia(ir)).^2;
0211             xylambda_b_p3m2(ir,:) = ~xymaskT(ir,:).*(ymhu2 - xia(ir)*(1 - xia(ir)/2)/(1 + xia(ir)))/(1 + xia(ir))./ymhu2; 
0212         end
0213     end
0214     xylambda_b_p2m2p2 = ~xymaskT.*xylambda_p2m2p2;
0215     xylambda_b_p2 = xylambda_b_p2m2p2;
0216     xylambda_b_p1m1p2 = xylambda_b_p1p1;
0217     xylambda_b_p1m3p4 = xylambda_b_p1p1;
0218     xylambda_b_p1m2p2 = ones(size(xylambda_b_p1m3p4));
0219  end      
0220