0001 function fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul)
0002
0003 epsi_kulsrud = [0.04;0.06;0.08;0.10];
0004 Zi_kulsrud =[1,2,3,10];
0005 RR_kulsrud = 1e-4*[ 0.01914,NaN ,NaN ,NaN ;...
0006 0.5411 ,0.2611,NaN ,NaN ;...
0007 3.177 ,1.735 ,1.047,0.09 ;...
0008 10.04 ,5.839 ,3.757,0.449];
0009
0010 iepsi = find(abs((alpha*betath2 - epsi_kulsrud)./epsi_kulsrud) < 10*eps);
0011
0012 if opt.coll_mode == 0 && opt.bounce_mode == 0 && abs((betath2 - 1e-6)/1e-6) < 10*eps && ~isempty(iepsi) && any(Zi == Zi_kulsrud),
0013 RR_k = RR_kulsrud(iepsi,Zi == Zi_kulsrud);
0014 else
0015 RR_k = '';
0016 end
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 if grid.np_tail == 0,
0028 pnmax_S = grid.pnmax_S;
0029 elseif imag(grid.pnmax_S_tail) == 0,
0030 pnmax_S = grid.pnmax_S_tail;
0031 else
0032 pnmax_S = grid.pnmax_S*imag(grid.pnmax_S_tail);
0033 end
0034
0035 format
0036
0037 if p_opt > 0,
0038 filename = ['res_',id_simul,'.log'];
0039
0040 if exist(filename,'file')
0041 delete(filename);
0042 end
0043
0044 diary(filename);
0045 end
0046
0047 nr = length(RR);
0048 np = length(Ec);
0049
0050 disp('Simulation parameters : ')
0051 disp(['coll_mode = ',num2str(opt.coll_mode),' ; bounce_mode = ',num2str(opt.bounce_mode),' ; equil = ',id_equil])
0052 disp(['Te/mc2 = ',num2str(betath2),' ; E/Ec = ',num2str(alpha),' ; E/ED = ',num2str(alpha*betath2),' ; Zi = ',num2str(Zi)])
0053 if opt.synchro_mode == 0,
0054 disp(['No synchrotron force'])
0055 elseif imag(wpr) == 0,
0056 disp(['Synchrotron force on : wc2/wp2 = ',num2str(wpr)]);
0057 elseif imag(wpr) > 0,
0058 disp(['Synchrotron force on : Bhat = ',num2str(imag(wpr))]);
0059 else
0060 disp(['Synchrotron force on : 1/taur = ',num2str(-imag(wpr))]);
0061 end
0062 disp(['NR limit : pc/pT = ',num2str(1/sqrt(alpha*betath2)),' ; pmax/pc = ',num2str(pnmax_S*sqrt(alpha*betath2))])
0063 disp(['Relativ. : pc/pT = ',num2str(1/sqrt((alpha - 1)*betath2)),' ; pmax/pc = ',num2str(pnmax_S*sqrt((alpha - 1)*betath2))])
0064 disp(' ')
0065 if nr == 1,
0066 disp(['--> r/a = ',num2str(rho_S)])
0067 disp(['--> LUKE RR from bound calc. (normalized to ne and tauc) : ',num2str(RR)])
0068 disp(['--> LUKE RR from exact calc. (normalized to ne and tauc) : ',num2str(tRR.re(end) + tRR.ri(end))])
0069 disp(['--> LUKE ext. RE fraction : ',num2str(tnorm.re(end))])
0070 disp(['--> LUKE int. RE fraction : ',num2str(tnorm.ri(end))])
0071 disp(['--> LUKE tot. RE fraction : ',num2str(tnorm.re(end) + tnorm.ri(end))])
0072
0073 if opt.synchro_mode,
0074 fpar = f{end}(:,end).';
0075 ipparnan = find(fpar < 0,1);
0076 if ~isempty(ipparnan),
0077 fpar(ipparnan:end) = NaN;
0078 end
0079 ipparmin = find(diff(fpar) > 0,1);
0080 if ~isempty(ipparmin),
0081 disp(['--> minimum found in f(xi=1) : p/pT = ',num2str(fgrid.pn(ipparmin)),' ; Ec = ',num2str(Ec(ipparmin)),' MeV']);
0082
0083 ipparmax = find(diff(fpar) < 0 & 1:np-1 > ipparmin,1);
0084
0085 if ~isempty(ipparmax),
0086 disp(['--> maximum found in f(xi=1) : p/pT = ',num2str(fgrid.pn(ipparmax)),' ; Ec = ',num2str(Ec(ipparmax)),' MeV']);
0087 end
0088 end
0089 end
0090 else
0091 rho = (rho_S(1:end-1) + rho_S(2:end))/2;
0092 disp('--> r/a | bound. RR | exact RR')
0093 disp(num2str([rho.',RR.',tRR.re(:,end) + tRR.ri(:,end)]))
0094 end
0095 if ~isempty(RR_k),
0096 disp(['--> Kulsrud RR : ',num2str(RR_k)])
0097 end
0098 disp('----------------------')
0099
0100 if p_opt > 0,
0101 diary off
0102 end
0103
0104 if opt.timevol == 0,
0105
0106 figure(1),clf
0107
0108 fM = f{1}(:,end);
0109 f = f{2}(:,end);
0110
0111 xlab = 'E_c (MeV)';
0112 ylab = 'f(E_c,\xi = 1)';
0113 xlim = [0,max(Ec)];
0114 if isfield(graph,'ylim_f'),
0115 ylim = graph.ylim_f;
0116 else
0117 ylim = [eps,max(fM)];
0118 graph.ytick_f = 10.^(-18:2:2);
0119 end
0120 tit = '';
0121 siz = 20+14i;
0122
0123 leg = {'f_M','f'};
0124 graph1D_jd(Ec,fM,0,1,'','','',NaN,xlim,ylim,'-','none','b',0.5,siz);
0125 graph1D_jd(Ec,f,0,1,xlab,ylab,tit,leg,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0126
0127 if isfield(graph,'ytick_f'),
0128 set(gca,'ytick',graph.ytick_f)
0129 end
0130 set(gca,'YMinorGrid','off')
0131 set(gca,'YMinorTick','on')
0132
0133 print_jd(p_opt,['fig_',id_simul,'_f1'],'./figures',1)
0134
0135 else
0136
0137 ntn = length(tn);
0138
0139 if isfield(graph,'tnc') && graph.tnc,
0140 tn = tn*betath2^(3/2);
0141 RRfac = betath2^(-3/2);
0142 else
0143 RRfac = 1;
0144 end
0145
0146 if opt.tnmin == 0,
0147 xlim = [0,tn(end)];
0148 xtick = NaN;
0149 mask = 1:ntn;
0150 else
0151 xlim = [tn(2),tn(end)];
0152 xtick = 10.^(floor(log10(tn(2))):2:ceil(log10(tn(end))));
0153 mask = 2:ntn;
0154 end
0155
0156 siz = 20+14i;
0157 tit = '';
0158 if ~isfield(opt,'color') || opt.color,
0159 colors = {'b','r',[0,0.75,0],[0.75,0,0.75],[0,0.75,0.75],[0.75,0.75,0],'k'};
0160
0161 styles = {'-','-','-','-','--','-','-'};
0162 sizes = [2,0.5,0.5,0.5,2];
0163 colstr = '_col';
0164 else
0165 colors = {'k','k','k','k','k','k','k','k','k'};
0166 styles = {'-','--','-.','-','--','-.'};
0167 colstr = '';
0168 sizes = [2,0.5,0.5,0.5,2];
0169 end
0170
0171 figure(1),clf
0172
0173 xlab = 't/\tau_c';
0174 ylab = 'n/n_e';
0175 if isfield(graph,'ylim_norm'),
0176 ylim = graph.ylim_norm;
0177 else
0178 ylim = 10.^[-9,1];
0179 graph.ytick_norm = 10.^(-8:2:0);
0180 end
0181
0182 leg = {'ext. RE','int. RE','bulk','total'};
0183
0184 graph1D_jd(tn(mask),tnorm.re(1,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,styles{3},'none',colors{3},sizes(3),siz,gca);
0185 graph1D_jd(tn(mask),tnorm.ri(1,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,styles{2},'none',colors{2},sizes(2),siz,gca);
0186
0187 graph1D_jd(tn(mask),tnorm.b(1,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,styles{4},'none',colors{4},sizes(4),siz,gca);
0188 graph1D_jd(tn(mask),tnorm.tot(1,mask),opt.tnmin > 0,1,xlab,ylab,tit,leg,xlim,ylim,styles{5},'none',colors{1},sizes(1),siz,gca,0.9,0.7,0.7);
0189
0190 if nr > 1,
0191 graph1D_jd(tn(mask),tnorm.re(end,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','r',0.5,siz,gca);
0192 graph1D_jd(tn(mask),tnorm.ri(end,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','b',0.5,siz,gca);
0193
0194 graph1D_jd(tn(mask),tnorm.b(end,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','g',0.5,siz,gca);
0195 graph1D_jd(tn(mask),tnorm.tot(end,mask),opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','k',0.5,siz,gca);
0196 end
0197
0198 if ~isnan(xtick),
0199 set(gca,'xtick',xtick)
0200 end
0201
0202 if isfield(graph,'ytick_norm'),
0203 set(gca,'ytick',graph.ytick_norm)
0204 end
0205 set(gca,'YMinorGrid','off')
0206 set(gca,'YMinorTick','on')
0207 set(gca,'XMinorGrid','off')
0208 set(gca,'XMinorTick','on')
0209
0210 figure(2),clf
0211
0212 xlab = 't/\tau_c';
0213 ylab = '\Gamma_R\tau_c';
0214 if isfield(graph,'ylim_gamma'),
0215 ylim = graph.ylim_gamma;
0216 else
0217 ylim = 10.^[-12,-2];
0218 graph.ytick_gamma = 10.^(-12:2:-2);
0219 end
0220
0221 leg = {'ext. RR','int. RR','total RR'};
0222 if ~isnan(RR_k),
0223 leg = ['Kulsrud',leg];
0224 end
0225
0226 if ~isempty(RR_k),
0227
0228 graph1D_jd(xlim,[RR_k,RR_k]*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0229 end
0230 graph1D_jd(tn(2:ntn),tRR.re(1,:)*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,styles{3},'none',colors{3},sizes(3),siz,gca);
0231 graph1D_jd(tn(2:ntn),tRR.ri(1,:)*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,styles{4},'none',colors{4},sizes(4),siz,gca);
0232 graph1D_jd(tn(2:ntn),tRR.re(1,:)*RRfac + tRR.ri(1,:)*RRfac,opt.tnmin > 0,1,xlab,ylab,tit,leg,xlim,ylim,styles{5},'none',colors{1},sizes(1),siz,gca,0.9,0.7,0.7);
0233
0234 if nr > 1,
0235 graph1D_jd(tn(2:ntn),tRR.re(end,:)*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','r',0.5,siz,gca);
0236 graph1D_jd(tn(2:ntn),tRR.ri(end,:)*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','b',0.5,siz,gca);
0237 graph1D_jd(tn(2:ntn),tRR.re(end,:)*RRfac + tRR.ri(end,:)*RRfac,opt.tnmin > 0,1,'','','',NaN,xlim,ylim,'--','none','k',0.5,siz,gca);
0238 end
0239
0240 if ~isnan(xtick),
0241 set(gca,'xtick',xtick)
0242 end
0243
0244 if isfield(graph,'ytick_gamma'),
0245 set(gca,'ytick',graph.ytick_gamma)
0246 end
0247 set(gca,'YMinorGrid','off')
0248 set(gca,'YMinorTick','on')
0249 set(gca,'XMinorGrid','off')
0250 set(gca,'XMinorTick','on')
0251
0252 figure(3),clf
0253
0254 if ~isfield(graph,'xlog'),
0255 graph.xlog = 0;
0256 end
0257
0258 fM = reshape(f{1}(:,end,:),[np,nr]);
0259
0260 if ~isfield(graph,'xec') || graph.xec,
0261 xlab = '$\mathcal{E}_k$ (MeV)';
0262 ylab = '$f(\mathcal{E}_k,\xi = 1)$';
0263
0264 x = Ec;
0265 else
0266 xlab = '$p/(mc)$';
0267 ylab = '$f(p,\xi = 1)$';
0268
0269 x = fgrid.pn*sqrt(betath2);
0270 end
0271
0272 if ~isfield(graph,'xlim'),
0273 if graph.xlog,
0274 graph.xlim = [min(x),max(x)];
0275 else
0276 graph.xlim = [0,max(x)];
0277 end
0278 end
0279
0280 if isfield(graph,'ylim_f'),
0281 ylim = graph.ylim_f;
0282 else
0283 ylim = [eps,max(fM(:))];
0284 graph.ytick_f = 10.^(-18:2:2);
0285 end
0286
0287 ntdisp = length(graph.itdisp);
0288 leg = cell(1,ntdisp+1);
0289 leg{1} = ['t/\tau_c = 0'];
0290 graph1D_jd(x,fM(:,1),graph.xlog,1,'','','',NaN,graph.xlim,ylim,styles{1},'none',colors{1},sizes(1),siz,gca);
0291 for iit = 1:ntdisp-1,
0292 leg{1+iit} = ['t/\tau_c = ',num2str(tn(1+graph.itdisp(iit)))];
0293 graph1D_jd(x,f{1+graph.itdisp(iit)}(:,end,1),graph.xlog,1,'','','',NaN,graph.xlim,ylim,styles{1+iit},'none',colors{1+iit},sizes(1+iit),siz);
0294 end
0295 leg{1+ntdisp} = ['t/\tau_c = ',num2str(tn(1+graph.itdisp(ntdisp)))];
0296 graph1D_jd(x,f{1+graph.itdisp(ntdisp)}(:,end,1),graph.xlog,1,xlab,ylab,tit,leg,graph.xlim,ylim,styles{1+ntdisp},'none',colors{1+ntdisp},sizes(1+ntdisp),siz,gca,0.9,0.7,0.7);
0297
0298 if nr > 1,
0299 graph1D_jd(x,fM(:,end),graph.xlog,1,'','','',NaN,graph.xlim,ylim,'--','none',colors{1},0.5,siz,gca);
0300 for iit = 1:ntdisp,
0301 graph1D_jd(x,f{1+graph.itdisp(iit)}(:,end,end),graph.xlog,1,'','','',NaN,graph.xlim,ylim,'--','none',colors{1+iit},0.5,siz);
0302 end
0303 end
0304
0305 if isfield(graph,'xtick'),
0306 set(gca,'xtick',graph.xtick)
0307 end
0308 if isfield(graph,'ytick_f'),
0309 set(gca,'ytick',graph.ytick_f)
0310 end
0311 set(gca,'XMinorGrid','off')
0312 set(gca,'XMinorTick','on')
0313 set(gca,'YMinorGrid','off')
0314 set(gca,'YMinorTick','on')
0315 set(get(gca,'xlabel'),'interpreter','latex')
0316 set(get(gca,'ylabel'),'interpreter','latex')
0317
0318 print_jd(p_opt,['fig_',id_simul,'_nr',colstr],'./figures',1)
0319 print_jd(p_opt,['fig_',id_simul,'_rr',colstr],'./figures',2)
0320 print_jd(p_opt,['fig_',id_simul,'_f1',colstr],'./figures',3)
0321
0322 if nr > 1,
0323
0324 figure(4),clf
0325
0326 xlim = [0,1];
0327 xlab = 'r/a';
0328 ylab = 'n/n_e';
0329 if isfield(graph,'ylim_norm'),
0330 ylim = graph.ylim_norm;
0331 else
0332 ylim = 10.^[-9,1];
0333 graph.ytick_norm = 10.^(-8:2:0);
0334 end
0335
0336 leg = {'ext. RE','int. RE','bulk','total'};
0337
0338 graph1D_jd(rho,tnorm.re(:,end),0,1,'','','',NaN,xlim,ylim,'-','none','r',2,siz,gca);
0339 graph1D_jd(rho,tnorm.ri(:,end),0,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0340 graph1D_jd(rho,tnorm.b(:,end),0,1,'','','',NaN,xlim,ylim,'-','none',[0,0.5,0],2,siz,gca);
0341 graph1D_jd(rho,tnorm.tot(:,end),0,1,xlab,ylab,tit,leg,xlim,ylim,'-','none','k',2,siz,gca,0.9,0.7,0.7);
0342
0343 set(gca,'xtick',0:0.2:1)
0344
0345 if isfield(graph,'ytick_norm'),
0346 set(gca,'ytick',graph.ytick_norm)
0347 end
0348 set(gca,'YMinorGrid','off')
0349 set(gca,'YMinorTick','on')
0350
0351 figure(5),clf
0352
0353 ylim = NaN;
0354
0355 leg = cell(1,ntdisp);
0356 for iit = 1:ntdisp-1,
0357 leg{iit} = ['t/\tau_c = ',num2str(tn(1+graph.itdisp(iit)))];
0358 graph1D_jd(rho,tnorm.re(:,1+graph.itdisp(iit)) + tnorm.ri(:,1+graph.itdisp(iit)),0,0,'','','',NaN,xlim,ylim,'-','none',colors{iit},2,siz);
0359 end
0360 leg{ntdisp} = ['t/\tau_c = ',num2str(tn(1+graph.itdisp(ntdisp)))];
0361 graph1D_jd(rho,tnorm.re(:,1+graph.itdisp(ntdisp)) + tnorm.ri(:,1+graph.itdisp(ntdisp)),0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{ntdisp},2,siz,gca,0.9,0.7,0.7);
0362
0363 set(gca,'xtick',0:0.2:1)
0364
0365
0366 elseif isfield(graph,'dp_cyl') && graph.dp_cyl > 0,
0367
0368 if ~isfield(graph,'pnmax_disp'),
0369 graph.pnmax_disp = pnmax_S;
0370 end
0371
0372 f02D = f{end};
0373 fM2D = f{1};
0374
0375 calc_fcyl_jd(f02D,sqrt(betath2),fgrid.pn,fgrid.mhu,graph.dp_cyl,1,1,graph.pnmax_disp,tacc{end},pmhucrit,[6,0,0]);
0376 set(gca,'xlim',graph.xlim_disp);
0377 set(gca,'xtick',graph.xtick_disp);
0378 set(gca,'ylim',graph.ylim_disp);
0379 set(gca,'ytick',graph.ytick_disp);
0380
0381 figpos = get(gcf,'position');
0382
0383 newfigpos = [figpos(1),figpos(2),figpos(3),figpos(3)*1/2];
0384 set(gcf,'position',newfigpos)
0385
0386
0387
0388 dke_out.ZXXS = ZXXS;
0389 dke_out.xRR_fsav = RR;
0390 dke_out.betath_ref = sqrt(betath2);
0391
0392 gridDKE.xmhubounce2 = 0;
0393
0394 stream_dke_jd(dke_out,opt,fgrid,gridDKE,1,graph.dp_cyl,graph,tacc{end},pmhucrit,7);
0395
0396
0397
0398 set(gca,'xlim',graph.xlim_disp);
0399 set(gca,'xtick',graph.xtick_disp);
0400 set(gca,'ylim',graph.ylim_disp);
0401 set(gca,'ytick',graph.ytick_disp);
0402 set(gcf,'position',newfigpos)
0403
0404 print_jd(p_opt,['fig_',id_simul,'_f2',colstr],'./figures',6)
0405 print_jd(p_opt,['fig_',id_simul,'_stream',colstr],'./figures',7)
0406
0407 end
0408
0409 end
0410 