Description of banana_dke

0001 function [YTdt,YTmask,YTF,YTPSIB,YTmhu0,YTmhu,YTx,xymaskT,nt2] = banana_dke_jd(ymhu0,ap,Xx,Xy,XBx,XBy,XBphi);
0002 %
0003 % Calculation of the poloidal positions of the two banana tips for accurate numerical integration over the poloidal angle
0004 % close to the trapped/passing limit. All output quantities are given on this new grid
0005 %   (see the notice for explanations concerning notations)
0006 %
0007 %   INPUT:
0008 %
0009 %       - ymhu0: pitch-angle grid
0010 %       - Plasma magnetic equilibrium parameters
0011 %
0012 %   OUTPUT:
0013 %
0014 %       - YTdt: dtheta
0015 %       - YTmask: Mask on the integration over B
0016 %       - YTF: B*mhu0*[(R-Rp)^2+(Z-Zp)^2]/|BR*(Z-Zp)-BZ*(R-Rp)|/mhu/ap, integrant that appears in the numerical bounce integrals (see notice)
0017 %       - YTPSIB: Ratio of B-field to its minimum on flux surface, B/Bmin
0018 %       - YTmhu0: Pitch-angle at Bmin
0019 %       - YTmhu: Pitch-angle cosine of (psi,theta,mhu) at B
0020 %       - YTR: Major radius of the flux surface (psi,theta) at B
0021 %       - xymaskT: Trapped domain
0022 %       - nt2: number of points in the new theta grid
0023 %
0024 % by Joan Decker <jodecker@mit.edu> (MIT/RLE) and  Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0025 %
0026 [np,nt] = size(Xx);
0027 nmhu = length(ymhu0);
0028 %
0029 theta = linspace(0,2*pi,nt);%Poloidal angular grid
0030 %
0031 XBT = abs(XBphi);%Toroidal B-field
0032 XBP = sqrt(XBx.^2 + XBy.^2);%Poloidal B-field
0033 XB = sqrt(XBT.^2 + XBP.^2);%Total B-field
0034 %
0035 xBmin = min(XB.'); % Minimum of B on FS
0036 xBmax = max(XB.'); % Maximum of B on FS
0037 %
0038 y1mmhu02 = 1 - ymhu0.^2;
0039 y1mmhu02(y1mmhu02 == 0) = eps;%Avoid numerical singularities
0040 xyBb = xBmin.'*(1./y1mmhu02);%B-field at turning point
0041 %
0042 xmhu0T = sqrt(1 - xBmin./xBmax);%Pitch-angle at the trapped-passing boundary
0043 %
0044 xymhu0 = repmat(ymhu0,[np,1]);
0045 xymhu0T = repmat(xmhu0T.',[1,nmhu]);
0046 %
0047 xymaskT = abs(xymhu0) <= xymhu0T - 10*eps;%Select trapped particles (the eps is for access to very core plasma calculations)
0048 %
0049 XF = zeros(nt,np,5);
0050 XF(:,:,1) = XBx.';
0051 XF(:,:,2) = XBy.';
0052 XF(:,:,3) = XBphi.';
0053 XF(:,:,4) = Xx.';
0054 XF(:,:,5) = Xy.';
0055 %
0056 [xythetaT1,xythetaT2,xyF1,xyF2] = interp_theta_jd(theta,XBx.',XBy.',XBphi.',xyBb,XF);
0057 %
0058 xythetaT1(~xymaskT) = theta(1);
0059 xythetaT2(~xymaskT) = theta(1);
0060 %
0061 if any(isnan(xythetaT1(:))) | any(isnan(xythetaT2(:))),
0062     error('The theta interpolation failed.');
0063 end
0064 %
0065 Ytheta = repmat(theta.',[1,np,nmhu]);
0066 %
0067 YBx = repmat(XBx.',[1,1,nmhu]);    
0068 YBy = repmat(XBy.',[1,1,nmhu]);
0069 YBphi = repmat(XBphi.',[1,1,nmhu]);
0070 Yx = repmat(Xx.',[1,1,nmhu]);
0071 Yy = repmat(Xy.',[1,1,nmhu]);
0072 %
0073 xyBxT1 = xyF1(:,:,1);
0074 xyBxT2 = xyF2(:,:,1);
0075 %
0076 xyByT1 = xyF1(:,:,2);
0077 xyByT2 = xyF2(:,:,2);
0078 %
0079 xyBphiT1 = xyF1(:,:,3);
0080 xyBphiT2 = xyF2(:,:,3);
0081 %
0082 xyxT1 = xyF1(:,:,4);
0083 xyxT2 = xyF2(:,:,4);
0084 %
0085 xyyT1 = xyF1(:,:,5);
0086 xyyT2 = xyF2(:,:,5);
0087 %
0088 xyBxT1(~xymaskT) = YBx(1,~xymaskT);
0089 xyBxT2(~xymaskT) = YBx(1,~xymaskT);
0090 %
0091 xyByT1(~xymaskT) = YBy(1,~xymaskT);
0092 xyByT2(~xymaskT) = YBy(1,~xymaskT);
0093 %
0094 xyBphiT1(~xymaskT) = YBphi(1,~xymaskT);
0095 xyBphiT2(~xymaskT) = YBphi(1,~xymaskT);
0096 %
0097 xyxT1(~xymaskT) = Yx(1,~xymaskT);
0098 xyxT2(~xymaskT) = Yx(1,~xymaskT);
0099 %
0100 xyyT1(~xymaskT) = Yy(1,~xymaskT);
0101 xyyT2(~xymaskT) = Yy(1,~xymaskT);
0102 %
0103 YTtheta = [Ytheta;reshape(xythetaT1,[1,np,nmhu]);reshape(xythetaT2,[1,np,nmhu])];
0104 %
0105 YTBx = [YBx;reshape(xyBxT1,[1,np,nmhu]);reshape(xyBxT2,[1,np,nmhu])];
0106 YTBy = [YBy;reshape(xyByT1,[1,np,nmhu]);reshape(xyByT2,[1,np,nmhu])];
0107 YTBphi = [YBphi;reshape(xyBphiT1,[1,np,nmhu]);reshape(xyBphiT2,[1,np,nmhu])];
0108 YTx = [Yx;reshape(xyxT1,[1,np,nmhu]);reshape(xyxT2,[1,np,nmhu])];
0109 YTy = [Yy;reshape(xyyT1,[1,np,nmhu]);reshape(xyyT2,[1,np,nmhu])];
0110 %
0111 clear YBx YBy YBphi Yx Yy Ytheta
0112 %
0113 nt2 = nt + 2;%Two poloidal angle values are added corresponding to banana tips
0114 %
0115 % Sort according to theta
0116 %
0117 [YTtheta,YTiT] = sort(YTtheta);
0118 %
0119 csum1 = repmat(0:np-1,[nt2,1,nmhu]);
0120 csum2 = repmat(reshape(0:nmhu-1,[1,1,nmhu]),[nt2,np,1]);
0121 YTiT = YTiT + nt2*(csum1 + np*csum2);
0122 %
0123 YTBx = YTBx(YTiT);
0124 YTBy = YTBy(YTiT);
0125 YTBphi = YTBphi(YTiT);
0126 YTx = YTx(YTiT);
0127 YTy = YTy(YTiT);
0128 %
0129 % Half grid in theta for integration
0130 %
0131 YTBx = (YTBx(1:nt2-1,:,:) + YTBx(2:nt2,:,:))/2;
0132 YTBy = (YTBy(1:nt2-1,:,:) + YTBy(2:nt2,:,:))/2;
0133 YTBphi = (YTBphi(1:nt2-1,:,:) + YTBphi(2:nt2,:,:))/2;
0134 YTx = (YTx(1:nt2-1,:,:) + YTx(2:nt2,:,:))/2;
0135 YTy = (YTy(1:nt2-1,:,:) + YTy(2:nt2,:,:))/2;
0136 %
0137 YTB = sqrt(YTBx.^2 + YTBy.^2 + YTBphi.^2); 
0138 %
0139 nt2 = nt2 - 1;
0140 %
0141 YTBmin = repmat(xBmin,[nt2,1,nmhu]);
0142 YTmhu0 = repmat(reshape(ymhu0,[1,1,nmhu]),[nt2,np,1]);
0143 %
0144 YTPSIB = YTB./YTBmin;%Ratio of B-field to its minimum on flux surface
0145 YTBb =  repmat(reshape(xyBb,[1,np,nmhu]),[nt2,1,1]);%B-field at turning point
0146 %
0147 YTmask = YTB < YTBb;%Mask on the integration over B
0148 %
0149 clear YTBb
0150 %
0151 YTmhu = sign(YTmhu0).*sqrt(1 - YTPSIB.*(1 - YTmhu0.^2));%Pitch-angle cosine of (psi,theta,mhu)
0152 YTmhu(~YTmask) = 0;
0153 YTmhu(YTmhu == 0) = eps;%Avoid numerical singularities
0154 YTF = YTB.*(YTx.^2 + YTy.^2).*YTmhu0./abs(YTBx.*YTy - YTBy.*YTx)./YTmhu/ap;
0155 YTF(find(isinf(YTF))) = 1;%To avoid Inf and NaN
0156 YTdt = diff(YTtheta);
0157 %
banana_dke_jd

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
