S_anal_en

 This function calculates the avalanche source term with analytical
 formula for equilibrium of circular concentric flux surfaces

 INPUT
    - equilDKE
    - gridDKE: grid structure, XXmhu
    - Zmomcoef

 OUTPUT:
    - XXalpha_an: The term alpha for the bounce averaged source term for
    secondary runaways analytically for equilibrium of circular concentric flux surfaces

 By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

0001 function [XXalpha_an,XXmask_ctsa,thetastar_an,thetastar_an2] = S_anal_en(equilDKE,gridDKE,Zmomcoef,pnRA,dpnRA,npnRA,gammaRA)
0002 %
0003 % This function calculates the avalanche source term with analytical
0004 % formula for equilibrium of circular concentric flux surfaces
0005 %
0006 % INPUT
0007 %    - equilDKE
0008 %    - gridDKE: grid structure, XXmhu
0009 %    - Zmomcoef
0010 %
0011 % OUTPUT:
0012 %    - XXalpha_an: The term alpha for the bounce averaged source term for
0013 %    secondary runaways analytically for equilibrium of circular concentric flux surfaces
0014 %
0015 % By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0016 %
0017 %
0018 npsi = gridDKE.nr;
0019 %npn = gridDKE.npn;
0020 nmhu = gridDKE.nmhu;
0021 %
0022 epsilon=equilDKE.ap*equilDKE.xrho/equilDKE.Rp;
0023 costhetastar_an=NaN(npsi,npnRA,nmhu);
0024 %
0025 for ipsi=1:npsi
0026     for m=1:nmhu
0027         costhetastar_an(ipsi,:,m)=1/epsilon(ipsi) *(((1-gridDKE.mhu(m).^2).*(gammaRA+1).*(1+epsilon(ipsi)))./2-1);
0028     end
0029 end
0030 %
0031 XXmask_ctsa = costhetastar_an>-1 & costhetastar_an<1;
0032 XXcosthetastar_an=costhetastar_an.*XXmask_ctsa;
0033 XXcosthetastar_an(XXcosthetastar_an==0) = NaN;
0034 thetastar_an=acos(XXcosthetastar_an);
0035 thetastar_an2=2*pi-acos(XXcosthetastar_an);
0036 %
0037 [qe,me,mp,mn,e0,mu0] = pc_dke_yp;
0038 %
0039 Ip=1; %plasma current [MA]
0040 D = equilDKE.xrho.*mu0*Ip*1e6./(equilDKE.ap*2*pi); %Bp in cylindrical case, from Amperes law
0041 %
0042 for ipsi=1:npsi
0043     for m=1:nmhu
0044 dPsi1(ipsi,:,m)=epsilon(ipsi).*sin(thetastar_an(ipsi,:,m)).*(1+epsilon(ipsi))./(1+epsilon(ipsi).*cos(thetastar_an(ipsi,:,m))).^2;
0045 dPsi2(ipsi,:,m)=epsilon(ipsi).*sin(thetastar_an2(ipsi,:,m)).*(1+epsilon(ipsi))./(1+epsilon(ipsi).*cos(thetastar_an2(ipsi,:,m))).^2;
0046 Bp1(ipsi,:,m)=D(ipsi)./(1+epsilon(ipsi).*cos(thetastar_an(ipsi,:,m)));
0047 Bp2(ipsi,:,m)=D(ipsi)./(1+epsilon(ipsi).*cos(thetastar_an2(ipsi,:,m)));  
0048     end
0049 XXalpha_an(ipsi,:,:)=equilDKE.xB0(ipsi)*XXmask_ctsa(ipsi,:,:).*equilDKE.xrho(ipsi)./abs(dPsi1(ipsi,:,:))./Bp1(ipsi,:,:)...
0050     + equilDKE.xB0(ipsi)*XXmask_ctsa(ipsi,:,:).*equilDKE.xrho(ipsi)./abs(dPsi2(ipsi,:,:))./Bp2(ipsi,:,:); %should be normalized to ap instead ?, now its Rp
0051 end