This function calculates the avalanche source term
INPUT
- gridDKE: grid structure, XXmhu
- display: dkedisplay structure
- Zmomcoef
- Zbouncecoef
- equilDKE
- radialDKE
OUTPUT:
- The full avalanche operator S(p,mhu,l) (source term for runaways)
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 [XXSavalanches_norm,xnrem_init,xnrep_init,xnre_out_norm] = avalanche_en(dkepath,dkeparam,dkedisplay,mksa,gridDKE,Zmomcoef,Zbouncecoef,ohm,equilDKE,radialDKE) 0002 % 0003 % This function calculates the avalanche source term 0004 % 0005 % INPUT 0006 % - gridDKE: grid structure, XXmhu 0007 % - display: dkedisplay structure 0008 % - Zmomcoef 0009 % - Zbouncecoef 0010 % - equilDKE 0011 % - radialDKE 0012 % 0013 % OUTPUT: 0014 % - The full avalanche operator S(p,mhu,l) (source term for runaways) 0015 % 0016 % By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker 0017 % (CEA-DRFC, joan.decker@cea.fr) 0018 % 0019 xrho = equilDKE.xrho; 0020 xrhoT = equilDKE.xrhoT; 0021 xrhoP = equilDKE.xrhoP; 0022 xrhoV = equilDKE.xrhoV; 0023 % 0024 if isfield(ohm,'xnrem_init') && isfield(ohm,'xnrep_init'), 0025 if length(ohm.rho) == 1,%local calculation or uniform value 0026 xnrem_init = ohm.xnrem_init*ones(size(xrho)); 0027 xnrep_init = ohm.xnrep_init*ones(size(xrho)); 0028 elseif ~isfield(ohm,'rhotype') || strcmp(ohm.rhotype,'g'),%geometric rho 0029 xnrem_init = interp1(ohm.rho,ohm.xnrem_init,xrho,'spline') 0030 xnrep_init = interp1(ohm.rho,ohm.xnrep_init,xrho,'spline') 0031 elseif strcmp(ohm.rhotype,'p'), 0032 xnrem_init = interp1(ohm.rho,ohm.xnrem_init,xrhoP,'spline') 0033 xnrep_init = interp1(ohm.rho,ohm.xnrep_init,xrhoP,'spline') 0034 elseif strcmp(ohm.rhotype,'t'), 0035 xnrem_init = interp1(ohm.rho,ohm.xnrem_init,xrhoT,'spline') 0036 xnrep_init = interp1(ohm.rho,ohm.xnrep_init,xrhoT,'spline') 0037 elseif strcmp(ohm.rhotype,'v'), 0038 xnrem_init = interp1(ohm.rho,ohm.xnrem_init,xrhoV,'spline') 0039 xnrep_init = interp1(ohm.rho,ohm.xnrep_init,xrhoV,'spline') 0040 else 0041 error('unknown rho type in ohm structure') 0042 end 0043 else 0044 xnrem_init = zeros(size(xrho)); 0045 xnrep_init = zeros(size(xrho)); 0046 end 0047 % 0048 if dkeparam.avalanche_mode == 1, %if 1, avalanche mode on 0049 % 0050 npsi = gridDKE.nr; 0051 npn = gridDKE.npn; 0052 nmhu = gridDKE.nmhu; 0053 % 0054 pnmin2_KO = dkeparam.pnmin2_KO; % low threshold on secondary electrons for the knock-on operator 0055 pnmax2_KO = dkeparam.pnmax2_KO; % high threshold on secondary electrons for the knock-on operator 0056 % 0057 %extending the grid so that we can count the avalanche created beyond pnmax_S, if pnmax2_KO > pnmax_S 0058 % 0059 pnRA = [gridDKE.pn,gridDKE.pn(end) + gridDKE.dpn(end):gridDKE.dpn(end):pnmax2_KO]; 0060 npnRA = length(pnRA); 0061 dpnRA = [gridDKE.dpn,gridDKE.dpn(end)*ones(1,npnRA - npn)]; 0062 % 0063 pnRA2 = pnRA.^2; 0064 XXpnRA2 = permute(repmat(pnRA.^2,[npsi,1,nmhu]),[2,3,1]); 0065 XXpnRA = permute(repmat(pnRA,[nmhu 1]),[2,1]); 0066 gammaRA = sqrt(1+(pnRA*mksa.betath_ref).^2); 0067 % 0068 if dkeparam.bounce_mode > 0, 0069 %[XXalpha_an,XXmask_ctsa,thetastar_an,thetastar_an2]=S_anal_en(equilDKE,gridDKE,Zmomcoef, pnRA, dpnRA, npnRA, gammaRA); %calculates alpha and theta* with analytical formula in circ concentric case 0070 tic 0071 [XXthetastar1,XXthetastar2,XXalpha] = thetastar_calc_dke_en(equilDKE,gridDKE,Zmomcoef, pnRA, dpnRA, npnRA, gammaRA);%calc the critical poloidal angle value theta_star for runaway avalanches with interp1 0072 toc 0073 else 0074 XXalpha=ones(npsi,npnRA,nmhu); 0075 XXthetastar1=0; 0076 XXthetastar2=0; 0077 XXalpha_an=ones(npsi,npnRA,nmhu); 0078 end 0079 % 0080 Lambd = mksa.lnc_e_ref; 0081 ne_ref=mksa.ne_ref; 0082 beta=mksa.betath_ref; 0083 ap=equilDKE.ap; 0084 Rp=equilDKE.Rp; 0085 %nr=0.05*ne; 0086 % 0087 B0=equilDKE.xB0; 0088 Bmax = max((sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2 + equilDKE.XBphi.^2)),[],2)'; 0089 % 0090 mhu0T=sqrt(1-B0./Bmax); 0091 % 0092 [~,~,~,~,~,xbetath] = normcoef_dke_yp(mksa,equilDKE,gridDKE); 0093 for ipsi=1:npsi 0094 if ~isnan(pnmin2_KO),% S(p>=pnmin2_KO) 0095 if imag(pnmin2_KO) == 0, 0096 pthmask = pnRA >= pnmin2_KO; %pmin normalized to reference thermal momentum 0097 else 0098 pthmask = pnRA >= imag(pnmin2_KO)*xbetath(ipsi)/mksa.betath_ref; %pmin normalized to actual thermal momentum 0099 end 0100 else 0101 pthmask = pnRA >= 0; % mask set in main_dke_yp using XXRintmask 0102 end 0103 % 0104 if pnmax2_KO <= gridDKE.pn(end), 0105 pthmask = pthmask & pnRA <= pnmax2_KO; 0106 end 0107 % 0108 Srp(ipsi,:) = beta^3*beta^3*1/(4*pi)/Lambd./(beta*pnRA.*gammaRA.*(gammaRA-1).^2).*pthmask; % rosenbluth (S*), first beta3 is from the Jacobian 0109 for ip0=1:length(pnRA) 0110 p0mask = pnRA >= pnRA(ip0); 0111 I(ipsi,ip0) = 2*pi*integral_dke_jd(dpnRA,p0mask.*Srp(ipsi,:).*pnRA2,2); 0112 end 0113 end 0114 % 0115 % figure;graph1D_jd(pnRA,I/I(1),0,0,'p>p0','integral over S [p0,pmax]','',NaN,NaN,NaN,'-','.','b',2); 0116 % 0117 XXSrp=permute(repmat(Srp, [1,1,nmhu]),[1,2,3]); 0118 % 0119 delt = zeros(npnRA,nmhu); 0120 % 0121 %delta function in the cylindrical + analytical case 0122 % 0123 for ip=1:npnRA 0124 xistar(ip) = sqrt((gammaRA(ip)-1)/(gammaRA(ip)+1)); 0125 dxistar=gridDKE.dmhu(1); 0126 xigrid = linspace(xistar(ip)-1/2*dxistar,xistar(ip)+1/2*dxistar,2); 0127 fdel = 1*nmhu/2; 0128 mhu_S = linspace(gridDKE.mhu(1)-gridDKE.dmhu(1)/2, gridDKE.mhu(end)+gridDKE.dmhu(end)/2,nmhu+1); 0129 delt(ip,:) = redistrib_jd(xigrid,fdel,mhu_S,0);%integral_dke_jd(gridDKE.dmhu',delt(1,ip,:),3) 0130 end 0131 delt = permute(repmat(delt, [1,1,npsi]),[3,1,2]); 0132 %[dummy,I] = min(abs(gridDKE.mhu - sqrt((gammaRA(ip)-1)/(gammaRA(ip)+1)))) 0133 %delt(:,ip,I) = 1; 0134 %delt(:,ip,nmhu-I+1) = 1; 0135 % 0136 if dkeparam.bounce_mode > 0, 0137 C=NaN(npsi,npnRA,nmhu); 0138 % 0139 for ipsi=1:npsi 0140 for m=1:nmhu 0141 C(ipsi,:,m)=2/pi./(gammaRA+1)./(Zbouncecoef.Xxlambda(m,ipsi).*Zbouncecoef.xqtilde(ipsi)).*abs(gridDKE.mhu(m))./(1-gridDKE.mhu(m).^2).^2; 0142 % 0143 %XXS_an=XXSrp.*C.*XXalpha_an; 0144 ximask = zeros(npsi,npnRA,nmhu); 0145 ximask(:,:,nmhu/2+1:end) = 1; %positiv xi 0146 XXbouncef=C.*XXalpha.*ximask; %numerical 0147 end 0148 end 0149 % 0150 else 0151 C=1*delt; %cylindrical case 0152 XXbouncef=C.*XXalpha; 0153 end 0154 % 0155 XXSl = XXSrp.*XXbouncef; %<S()> numerical 0156 % 0157 XXSrp_xis=XXSrp.*delt; %3D analytical rosenbluth operator, non bounce averaged 0158 XXSl2=XXSrp_xis; 0159 XXSavalanches_normtot = permute(XXSl,[2,3,1]); 0160 % 0161 %check the operator, compare with theory RP NF 1997 0162 alpha = 1:0.5:20;%E/Ec 0163 Zi = -1; 0164 betath2 = mksa.betath_ref^2;% 0165 % 0166 for ia = 1:length(alpha) 0167 pc = sqrt(2)/sqrt((alpha(ia)-1)*betath2); 0168 %calculate secondary growth rate pmax<p<pext NORMALIZED to nr 0169 XXRintmask = XXpnRA >= pc; 0170 XXmhu0mask = gridDKE.mhup >= mhu0T; 0171 Srp_int_norm_nr(ia) = mksa.nhu_ref*2*pi*(Zbouncecoef.xqtilde(1,1)./Zbouncecoef.xqhat(1,1)).*integral_dke_jd(gridDKE.dmhu,XXmhu0mask'.*Zbouncecoef.Xxlambda.*integral_dke_jd(dpnRA,XXSavalanches_normtot(:,:,1).*XXRintmask.*XXpnRA2(:,:,1),1),1)'; % 0172 end 0173 % 0174 RR_an = RR_rosenbluth_en(alpha, Zi, mksa.ne_ref, mksa.lnc_e_ref, equilDKE.ap*equilDKE.xrho(1)/equilDKE.Rp); %analytic AVALANCHE growth rate Rosenbluth Putvinski NF 1997 0175 % 0176 % figure;graph1D_jd(alpha,Srp_int_norm_nr./RR_an,0,0,'E/E_c','growth rate (s^-1)','Numerical growth rate',NaN,NaN,NaN,'-','o','b',2);hold on 0177 % graph1D_jd(alpha,RR_an./RR_an,0,0,'','','',NaN, NaN,NaN,'-','none','r',3);hold off 0178 % figure;graph1D_jd(alpha,Srp_int_norm_nr,0,0,'E/E_c','growth rate (s^-1)','Numerical growth rate',NaN,NaN,NaN,'-','o','b',2);hold on 0179 % graph1D_jd(alpha,RR_an,0,0,'','','',NaN, NaN,NaN,'-','none','r',3);hold off 0180 % 0181 if npnRA > npn, 0182 xnre_out_norm(1,gridDKE.rdke) = 2*pi*(Zbouncecoef.xqtilde./Zbouncecoef.xqhat).*integral_dke_jd(gridDKE.dmhu,Zbouncecoef.Xxlambda.*integral_dke_jd(dpnRA((npn+1):end),XXSavalanches_normtot((npn+1):end,:,gridDKE.rdke).*XXpnRA2((npn+1):end,:,gridDKE.rdke),1),1)'; % runaway electrons created outisde pmax THIS IS A RATE! (per collision time) 0183 else 0184 xnre_out_norm(1,gridDKE.rdke) = 0; 0185 end 0186 % 0187 XXSavalanches_norm = XXSavalanches_normtot(1:npn,:,:); %the avalanche part with p < pmax 0188 % 0189 % 0190 %display section: 0191 % if dkedisplay.display_mode>0 0192 %ir=1; %nr of the flux surface to be displayed 0193 % if Zbouncecoef.XXlambda~=1 0194 % 0195 % ir=1 0196 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar1(ir,:,:)),'p','mhu','theta^*_1',NaN,NaN,0,NaN,30); colorbar; 0197 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar2(ir,:,:)),'p','mhu','theta^*_2',NaN,NaN,0,NaN,30); colorbar; 0198 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXS(ir,:,:)),'p','mhu','S',NaN,NaN,0,NaN,30); colorbar; 0199 % for ipsi=1:npsi 0200 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXbouncef(ipsi,:,:)),'p','mhu','Bouncefactor',NaN,NaN,0,NaN,100); colorbar; 0201 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXalpha(ipsi,:,:)),'p','mhu','XXalpha',NaN,NaN,0,NaN,100); colorbar; 0202 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(C(ipsi,:,:)),'p','mhu','C',NaN,NaN,0,NaN,100); colorbar; 0203 % end 0204 0205 % for ipsi=1:npsi 0206 % figure;graph2D_jd(gridDKE.xpn.*gridDKE.Xmhu,gridDKE.xpn.*sqrt(1-gridDKE.Xmhu.^2),squeeze(XXbouncef(ipsi,:,:))','p_{par}/p_{th}','p_{perp}/p_{th}','BF',NaN,NaN,0,NaN,30); colorbar; 0207 % hold on 0208 % h = plot(gridDKE.pn.*mhuT(ipsi),gridDKE.pn.*sqrt(1-mhuT(ipsi).^2)); 0209 % %h2 = plot(-gridDKE.pn.*mhuT(ipsi),gridDKE.pn.*sqrt(1-mhuT(ipsi).^2)); 0210 % set(h,'Color','k','LineWidth',2) 0211 % %set(h2,'Color','k','LineWidth',2) 0212 % end 0213 % 0214 % 0215 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(thetastar_an2(ir,:,:)),'p','mhu','theta^*_2 anal.',NaN,NaN,0,NaN,30); colorbar; 0216 % figure;graph2D_jd(pnRA,gridDKE.mhu,squeeze(XXalpha_an(ir,:,:)),'p','mhu','XXalpha analytical',NaN,NaN,0,NaN,100); colorbar; 0217 % figure;graph2D_jd(gridDKE.xpn.*gridDKE.XXmhu,gridDKE.xpn.*sqrt(1-gridDKE.XXmhu.^2),squeeze(thetastar_an(ir,:,:))','p_{par}/p_{th}','p_{perp}/p_{th}','theta^*_1 anal.',NaN,NaN,0,NaN,30); colorbar; 0218 % %figure;graph2D_jd(gridDKE.pn(:),gridDKE.XXmhu(1,:,l),XXthetastar2(:,:,l),'p','mhu','theta^*_2',NaN,NaN,0,NaN,30);colorbar; 0219 % %figure;graph2D_jd(gridDKE.pn(:),gridDKE.XXmhu(1,:,l),log(XXS(:,:,l)),'p','mhu','log(S)',NaN,NaN,1,NaN,30);colorbar; 0220 % %figure;graph2D_jd(Zmomcoef.xz.*gridDKE.XXmhu,Zmomcoef.xz.*sqrt(1-gridDKE.XXmhu.^2),log(XXS(:,:,l)),'p_{par}/p_{th}','p_{perp}/p_{th}','log(S)',NaN,NaN,1,NaN,30);colorbar; 0221 % %figure;graph2D_jd(gridDKE.pn(:),gridDKE.XXmhu(1,:,l),C(:,:,l),'p/p_{th}','mhu','C',NaN,NaN,0,NaN,30); 0222 % %figure;graph2D_jd(gridDKE.pn(:),gridDKE.XXmhu(1,:,l),XXalpha(:,:,l),'p/p_{th}','mhu','XXalpha',NaN,NaN,0,NaN,50);colorbar; 0223 % %check that the thetas appear in the same place in matrix: 0224 % figure;hold on 0225 % for n=1:npsi 0226 % graph1D_jd(1:nmhu,sum(~isnan(squeeze(XXthetastar1(n,:,:)))),0,0,'','','',NaN,NaN,NaN,'-','o','b') 0227 % graph1D_jd(1:nmhu,sum(~isnan(squeeze(XXthetastar2(n,:,:)))),0,0,'column in theta matrix','nr of non NaN per row','',NaN,NaN,NaN,'-','.','r') 0228 % end 0229 % % figure; graph2D_jd(gridDKE.pn(:).*gridDKE.XXmhu(1,:,l),gridDKE.pn(:).*sqrt(1-gridDKE.XXmhu(1,:,l).^2),XXthetastar1(:,:,l)'/pi,'p_{par}/p_{th}','p_{perp}/p_{th}','theta^*_1/pi',NaN,NaN,0,NaN,60);colorbar; 0230 % % figure; graph2D_jd(gridDKE.pn(:).*gridDKE.XXmhu(1,:,l),gridDKE.pn(:).*sqrt(1-gridDKE.XXmhu(1,:,l).^2),XXthetastar2(:,:,l)'/pi,'p_{par}/p_{th}','p_{perp}/p_{th}','theta^*_2/pi',NaN,NaN,0,NaN,60);colorbar; 0231 % end 0232 % % 0233 % figure 0234 % ylabel('XXalpha') 0235 % hold on 0236 % for h=1:nmhu 0237 % plot(XXalpha(1,:,h),'r*') 0238 % end 0239 % figure 0240 % ylabel('XXalpha anal.') 0241 % hold on 0242 % for h=1:nmhu 0243 % plot(XXalpha_an(1,:,h),'r*') 0244 % end 0245 % figure 0246 % for g=1:npsi 0247 % XS=squeeze(XXS(g,:,:)); 0248 % XS(isnan(XS)) = 0; 0249 % XSint = integral_dke_jd(gridDKE.dmhu,XS,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0250 % graph1D_jd(gridDKE.pn,XSint,0,1,'p','S_{bounce}','S on each flux surface',NaN,NaN,NaN,'-','.') 0251 % hold all, 0252 % end 0253 % 0254 % figure 0255 % for g=1:npsi 0256 % XS=squeeze(XXS(g,:,:)); 0257 % XS(isnan(XS)) = 0; 0258 % XSint = integral_dke_jd(gridDKE.pn,XS,1)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0259 % graph1D_jd(gridDKE.mhu,XSint,0,0,'mhu','S_{bounce}','S on each flux surface',NaN,NaN,NaN,'-','.') 0260 % hold all, 0261 % end 0262 % figure 0263 % for g=1:npsi 0264 % XSan=squeeze(XXS_an(g,:,:)); 0265 % XSan(isnan(XSan)) = 0; 0266 % XSan_int = integral_dke_jd(gridDKE.dmhu,XSan,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0267 % graph1D_jd(gridDKE.pn,XSan_int,0,1,'p','S_{bounce} analytical','S on each flux surface',NaN,NaN,NaN,'-','.') 0268 % end 0269 % %compare analytic alpha 0270 % figure 0271 % for g=1:npsi 0272 % XSan=squeeze(XXS_an(g,:,:)); 0273 % XSan(isnan(XSan)) = 0; 0274 % XSan_int = integral_dke_jd(gridDKE.pn,XSan,1)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0275 % graph1D_jd(gridDKE.mhu,XSan_int,0,0,'mhu','S_{bounce} anal.','S on each flux surface',NaN,NaN,NaN,'-','.') 0276 % hold all, 0277 % end 0278 % 0279 % figure 0280 % % 0281 % for g=1:npsi 0282 % XS = squeeze(XXSl(g,:,:)); 0283 % XS(isnan(XS)) = 0; 0284 % XSint = integral_dke_jd(gridDKE.dmhu,XS,2); 0285 % %semilogy(gridDKE.pn,XSint,'.-') 0286 % %legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(ir))); 0287 % graph1D_jd(pnRA,XSint,0,1,'p','S_{bounce} LUKE','S on each flux surface',NaN,NaN,NaN,'-','.') 0288 % hold all, 0289 % end 0290 % semilogy(pnRA, Srp(g,:),'rx');title('Rosenbluth theory (red x)'); xlabel('p'); ylabel('S') 0291 % % 0292 % figure 0293 % % 0294 % for g=1:npsi 0295 % XSrp = squeeze(XXSrp_xis(g,:,:)); %integrate the 3D analytical S_RP, should get same as plotting the 2D Srp 0296 % XSrpint = integral_dke_jd(gridDKE.dmhu,XSrp,2);%/(gridDKE.dmhu(1)); %integrate over mhu, pos n neg : /2 0297 % %semilogy(gridDKE.pn,XSint,'.-') 0298 % %legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(ir))); 0299 % graph1D_jd(pnRA,XSrpint,0,1,'p','S_{RP}','S on each flux surface',NaN,NaN,NaN,'-','.') 0300 % hold all, 0301 % end 0302 % semilogy(pnRA, Srp(g,:),'rx');title('S_RP compared to Rosenbluth theory (red x)'); xlabel('p'); ylabel('S') 0303 %as fcn of mhu, for each slize of p 0304 % figure 0305 % hold on 0306 % for ip=1:npn 0307 % XS = squeeze(XXS(1,:,:)); 0308 % graph1D_jd(gridDKE.mhu,XS(ip,:),0,0,'mhu','S','',NaN,NaN,NaN,'-','.','b') 0309 % end 0310 % % 0311 % figure 0312 % for g=1:npsi 0313 % Xbouncef = squeeze(XXbouncef(g,:,:)); 0314 % Xbouncef(isnan(Xbouncef)) = 0; 0315 % Xbouncefint = integral_dke_jd(ones(1,length(gridDKE.dmhu)),Xbouncef,2); %integrate over mhu, pos n neg : /2 0316 % graph1D_jd(gridDKE.pn,Xbouncefint,0,0,'p','Bounce factor','',NaN,NaN,NaN,'-','.','b') 0317 % legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(g))); 0318 % %graph1D_jd(Zmomcoef.xz,XS,0,1,'S_{bounce}','p','S on each flux surface',NaN,NaN,NaN,'-','.') 0319 % hold all, 0320 % end 0321 % figure 0322 % for g=1:npsi 0323 % Xalpha = squeeze(XXalpha(g,:,:)); 0324 % Xalpha(isnan(Xalpha)) = 0; 0325 % Xalphaint = integral_dke_jd(gridDKE.dmhu,Xalpha,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0326 % graph1D_jd(gridDKE.pn,Xalphaint,0,0,'p','Xalpha','',NaN,NaN,NaN,'-','.','b') 0327 % legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(g))); 0328 % %graph1D_jd(Zmomcoef.xz,XS,0,1,'S_{bounce}','p','S on each flux surface',NaN,NaN,NaN,'-','.') 0329 % hold all, 0330 % end 0331 % figure 0332 % for g=1:npsi 0333 % Xalpha_anal=squeeze(XXalpha_an(g,:,:)); 0334 % Xalpha_anal(isnan(Xalpha_anal)) = 0; 0335 % Xalpha_anal_int = integral_dke_jd(gridDKE.dmhu,Xalpha_anal,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0336 % graph1D_jd(gridDKE.pn,Xalpha_anal_int,0,0,'p','Xalpha analytical','alpha on each flux surface',NaN,NaN,NaN,'-','.') 0337 % end 0338 % figure 0339 % for g=1:ir 0340 % Xmask = squeeze(XXmask_ctsa(g,:,:)); 0341 % Xmaskint = integral_dke_jd(gridDKE.dmhu,Xmask,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0342 % graph1D_jd(gridDKE.pn,Xmaskint,0,0,'p','Xmask','',NaN,NaN,NaN,'-','.','b') 0343 % legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(ir))); 0344 % %graph1D_jd(Zmomcoef.xz,XS,0,1,'S_{bounce}','p','S on each flux surface',NaN,NaN,NaN,'-','.') 0345 % hold all, 0346 % end 0347 % figure 0348 % hold on 0349 % for ip=1:npn 0350 % Xalpha = squeeze(XXalpha(3,:,:)); 0351 % graph1D_jd(gridDKE.mhu,Xalpha(ip,:),0,0,'mhu','Xalpha','',NaN,NaN,NaN,'-','.','b') 0352 % end 0353 % figure 0354 % hold on 0355 % for ip=1:npn 0356 % Xalphaan = squeeze(XXalpha_an(1,:,:)); 0357 % graph1D_jd(gridDKE.mhu,Xalphaan(ip,:),0,0,'mhu','Xalpha anal.','',NaN,NaN,NaN,'-','.','b') 0358 % end 0359 % figure 0360 % hold on 0361 % for g=1:npsi 0362 % XC = squeeze(C(g,:,:)); 0363 % XC(isnan(XC)) = 0; 0364 % XCint = integral_dke_jd(gridDKE.dmhu,XC,2)/sum(gridDKE.dmhu); %integrate over mhu, pos n neg : /2 0365 % graph1D_jd(gridDKE.pn,XCint,0,0,'p','XC','',NaN,NaN,NaN,'-','.','b') 0366 % legend(sprintf( '%s: %d', 'rho', gridDKE.xrho(g))); 0367 % %graph1D_jd(Zmomcoef.xz,XS,0,1,'S_{bounce}','p','S on each flux surface',NaN,NaN,NaN,'-','.') 0368 % hold all, 0369 % end 0370 % figure 0371 % hold on 0372 % for ip=1:npn 0373 % XC = squeeze(C(3,:,:)); 0374 % graph1D_jd(gridDKE.mhu,XC(ip,:),0,0,'mhu','C','',NaN,NaN,NaN,'-','.','b') 0375 % end 0376 % 0377 % end 0378 else %if no avalanche mode is specified, zero contribution 0379 XXSavalanches_norm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr); 0380 xnre_out_norm = zeros(1,gridDKE.nr); 0381 end