jseries_dke_jd

 Calculates series J2m for bounce integrals
 
   INPUTS:
       
       - mhutrap: cosine of the trapping pitch angle [1,1]
       - mhu: grid in pitch angle cosines [1,nmhu]
       - m: m-1 is the maximun order for the series [1,1]

   OUTPUTS:
       
       - XJ: matrix of the series [m+1,nmhu]

 Joan Decker (jodecker@mit.edu) and Yves Peysson (yves.peysson@cea.fr)

0001 function XJ = jseries_dke_jd(mhutrap,mhu,m)
0002 %
0003 % Calculates series J2m for bounce integrals
0004 %
0005 %   INPUTS:
0006 %
0007 %       - mhutrap: cosine of the trapping pitch angle [1,1]
0008 %       - mhu: grid in pitch angle cosines [1,nmhu]
0009 %       - m: m-1 is the maximun order for the series [1,1]
0010 %
0011 %   OUTPUTS:
0012 %
0013 %       - XJ: matrix of the series [m+1,nmhu]
0014 %
0015 % Joan Decker (jodecker@mit.edu) and Yves Peysson (yves.peysson@cea.fr)
0016 %
0017 nmhu = length(mhu);
0018 XJ = zeros(m+1,nmhu);
0019 k = zeros(1,nmhu);
0020 e = zeros(1,nmhu);
0021 r = zeros(1,nmhu);
0022 r2 = zeros(1,nmhu);
0023 mhu2 = mhu.*mhu;
0024 mhutrap2 = mhutrap^2;
0025 %
0026 if mhutrap < 0 or mhutrap > 1,
0027     error('mhutrap must be between 0 and 1')
0028 end
0029 %
0030 if prod(mhu-mhutrap) == 0,
0031     imhu = find(mhu2 == mhutrap2);
0032     dmhu = diff(mhu);
0033      mhutrap = mhutrap - eps;%In order to avoid numerical singularity at mhutrap
0034      mhutrap2 = mhutrap^2;
0035 end
0036 %
0037 maskm = mhu2<mhutrap2;
0038 maskp = mhu2>mhutrap2;
0039 r = maskm.*mhu/mhutrap + maskp.*mhutrap./mhu;
0040 r2 = r.^2;
0041 %
0042 [k,e] = ellipke(r2);
0043 %
0044 XJ(1,:) = maskm.*abs(r).*k + maskp.*k;
0045 if m>0,
0046     XJ(2,:) = maskm.*abs(r).*(k-e) + maskp./r2.*(k-e);
0047     if m>1,
0048         for i=2:m,
0049             XJ(i+1,:) = 1/(2*i - 1)*((2*i - 2)*(1 + mhu2/mhutrap2).*XJ(i,:) - (2*i - 3)*mhu2/mhutrap2.*XJ(i-1,:));
0050         end
0051     end
0052 end
0053 if sum(sum(isnan(XJ)))~=0,error('NaN results in the J2m functions for bounce integrals');end
0054