detrapmj3

 calculations of losscone angle and detrapping momentum

0001 function [thetabanane,thetaloss,pndetrap,Eseuil,delta] = detrapmj3(zeff,a,B0,R0,BeLi,Ip,ne,rho)
0002 
0003 % calculations of losscone angle and detrapping momentum
0004 
0005 %rho; Normalized magnetic surface (r/a)
0006 %a; Plasma minor radius (m)
0007 %R0; Major radius (m)
0008 %zeff
0009 %Ip; Plsama current (MA)
0010 %BeLi = 'beta + Li/2'
0011 %delta; Magnetic ripple depth
0012 %Eseuil; Threshold energy (keV)
0013 %thetabanane; Banana angle w.r.t. pperp direction (deg)
0014 %thetaloss; Magnetic ripple loss cone angle w.r.t. pperp direction (deg)
0015 %pndetrap; Normalized magnetic ripple detrapping momentum
0016 %ne; (/m3)
0017 %D0; Shaf. shift (m)
0018 D0    =  a^2*(BeLi-0.5)/(2*R0);   
0019 
0020 %magnetic ripple depth
0021    delta = deltacalmj3(R0 + a*rho,zeros(size(rho)));
0022 %
0023 % Calculation of loss cone angle 'thetaloss'
0024     thetaconemj3;
0025 
0026 % Calculation of threshold energy and normalised detrapping momentum 'pndetrap'
0027 
0028 E      = [0:1:350];
0029 le     = length(E);
0030 
0031 for k2=1:le;
0032     F(k2) = fracmircalcmj3(E(k2),delta,ne,zeff,a,R0,B0,NBOB);
0033 end 
0034 indF      = find(F<=1);
0035 Eseuil = sum(E(indF).*(1-F(indF)))/sum(1-F(indF));     
0036 
0037 
0038 %if thetaZBC~=0
0039    Es    = Eseuil;
0040    pndetrap  = sqrt((Es/511+1).^2-1);
0041 %end