fproc_runaway_scan_2D

0001 function fproc_runaway_scan_2D(scanstr1,scanstr2,scanvar,scanval,scanlab,RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,iscandisp,xlim,xlog,xtick,p_opt,id_simul_scan,nan_enforce)
0002 %
0003 nscan1 = length(scanval.(scanstr1));
0004 nscan2 = length(scanval.(scanstr2));
0005 %
0006 snorm.re = NaN(nscan1,nscan2);
0007 snorm.ri = NaN(nscan1,nscan2);
0008 snorm.rtot = NaN(nscan1,nscan2);
0009 snorm.b = NaN(nscan1,nscan2);
0010 snorm.tot = NaN(nscan1,nscan2);
0011 %
0012 sRR.re = NaN(nscan1,nscan2);
0013 sRR.ri = NaN(nscan1,nscan2);
0014 sRR.rtot = NaN(nscan1,nscan2);
0015 %
0016 spnmin = NaN(nscan1,nscan2);
0017 spnmax = NaN(nscan1,nscan2);
0018 sEcmin = NaN(nscan1,nscan2);
0019 sEcmax = NaN(nscan1,nscan2);
0020 %
0021 salpha = NaN(nscan1,nscan2);
0022 sbetath2 = NaN(nscan1,nscan2);
0023 %
0024 spnmax_S = NaN(nscan1,nscan2);
0025 spnmax_S_tail = NaN(nscan1,nscan2);
0026 %
0027 sfpar = cell(nscan1,nscan2);
0028 %
0029 for iscan1 = 1:nscan1,
0030     for iscan2 = 1:nscan2,
0031         %
0032         snorm.re(iscan1,iscan2) = tnorm{iscan1,iscan2}.re(end);
0033         snorm.ri(iscan1,iscan2) = tnorm{iscan1,iscan2}.ri(end);
0034         snorm.rtot(iscan1,iscan2) = tnorm{iscan1,iscan2}.re(end) + tnorm{iscan1,iscan2}.ri(end);
0035         snorm.b(iscan1,iscan2) = tnorm{iscan1,iscan2}.b(end);
0036         snorm.tot(iscan1,iscan2) = tnorm{iscan1,iscan2}.tot(end);
0037         %
0038         sRR.re(iscan1,iscan2) = tRR{iscan1,iscan2}.re(end);
0039         sRR.ri(iscan1,iscan2) = tRR{iscan1,iscan2}.ri(end);
0040         sRR.rtot(iscan1,iscan2) = tRR{iscan1,iscan2}.re(end) + tRR{iscan1,iscan2}.ri(end);
0041         %
0042         if ~isempty(f{iscan1,iscan2}{end}),
0043             sfpar{iscan1,iscan2} = f{iscan1,iscan2}{end}(:,end).';
0044             np = length(sfpar{iscan1,iscan2});
0045             ipparnan = find(sfpar{iscan1,iscan2} < 0,1);
0046             if ~isempty(ipparnan),
0047                 sfpar{iscan1,iscan2}(ipparnan:end) = NaN;
0048             end
0049             ipparmin = find(diff(sfpar{iscan1,iscan2}) > 0,1);
0050             if ~isempty(ipparmin),
0051                 spnmin(iscan1,iscan2) = fgrid{iscan1,iscan2}.pn(ipparmin);
0052                 sEcmin(iscan1,iscan2) = Ec{iscan1,iscan2}(ipparmin);
0053                 %
0054                 ipparmax = find(diff(sfpar{iscan1,iscan2}) < 0 & 1:np-1 > ipparmin,1);
0055                 if ~isempty(ipparmax),
0056                     spnmax(iscan1,iscan2) = fgrid{iscan1,iscan2}.pn(ipparmax);
0057                     sEcmax(iscan1,iscan2) = Ec{iscan1,iscan2}(ipparmax);
0058                 end
0059             end
0060         end
0061         if strcmp(scanstr1,'alpha'),
0062             salpha(iscan1,iscan2) = scanval.(scanstr1)(iscan1);
0063         elseif strcmp(scanstr2,'alpha'),
0064             salpha(iscan1,iscan2) = scanval.(scanstr2)(iscan2);
0065         else
0066             salpha(iscan1,iscan2) = alpha;
0067         end
0068         if strcmp(scanstr1,'betath2'),
0069             sbetath2(iscan1,iscan2) = scanval.(scanstr1)(iscan1);
0070         elseif strcmp(scanstr2,'betath2'),
0071             sbetath2(iscan1,iscan2) = scanval.(scanstr2)(iscan2);
0072         else
0073             sbetath2(iscan1,iscan2) = betath2;
0074         end
0075         %
0076         if strcmp(scanstr1,'grid.pnmax_S'),
0077             spnmax_S(iscan1,iscan2) = scanval.(scanstr1)(iscan1);
0078         else
0079             spnmax_S(iscan1,iscan2) = grid.pnmax_S;
0080         end
0081         if strcmp(scanstr1,'grid.pnmax_S_tail'),
0082             spnmax_S_tail(iscan1,iscan2) = scanval.(scanstr1)(iscan1);
0083         else
0084             spnmax_S_tail(iscan1,iscan2) = grid.pnmax_S_tail;
0085         end
0086     end
0087 end
0088 %
0089 for ienforce = 1:size(nan_enforce,1),
0090     spnmin(nan_enforce(ienforce,1),nan_enforce(ienforce,2)) = NaN;
0091     spnmax(nan_enforce(ienforce,1),nan_enforce(ienforce,2)) = NaN;
0092 end
0093 %
0094 spc = 1./sqrt((salpha - 1).*sbetath2);
0095 sgammac = sqrt(salpha./(salpha - 1));
0096 %
0097 [~,~,~,~,~,~,~,mc2] = pc_dke_yp;%Universal physics constants
0098 sEcc = 1e-3*mc2*(sgammac - 1);
0099 %
0100 if grid.np_tail > 0,
0101     if imag(grid.pnmax_S_tail) == 0,
0102         spnmax_S = spnmax_S_tail;
0103     else
0104         spnmax_S = spnmax_S.*imag(spnmax_S_tail);
0105     end
0106 end
0107 %
0108 sgammamax_S = sqrt(1 + spnmax_S.^2.*sbetath2);
0109 sEcmax_S = 1e-3*mc2*(sgammamax_S - 1);
0110 %
0111 siz = 20+14i;
0112 colors = {'b','r','g','m','y','c','k'};
0113 %
0114 x1 = scanval.(scanstr1);
0115 x2 = scanval.(scanstr2);
0116 %
0117 xlim1 = NaN;
0118 xlim2 = NaN;
0119 %
0120 xlog1 = xlog.(scanstr1);
0121 xlog2 = xlog.(scanstr2);
0122 %
0123 xtick1 = xtick.(scanstr1);
0124 xtick2 = xtick.(scanstr2);
0125 %
0126 figure(1),clf
0127 %
0128 xlab1 = scanlab.(scanstr1);
0129 xlab2 = scanlab.(scanstr2);
0130 tit = 'n_r/n_e';
0131 %
0132 y = -log10(snorm.rtot);
0133 cont = -log10(graph.cont_norm);
0134 cont_lab = -log10(graph.cont_norm_lab);
0135 %
0136 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);%colorbar
0137 hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,0,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0138 %
0139 set(gca,'xtick',xtick1,'ytick',xtick2)
0140 %
0141 figure(2),clf
0142 %
0143 tit = '\Gamma_R\tau_c';
0144 %
0145 y = -log10(sRR.rtot);
0146 cont = -log10(graph.cont_gamma);
0147 cont_lab = -log10(graph.cont_gamma_lab);
0148 %
0149 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);%colorbar
0150 hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,0,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0151 %
0152 set(gca,'xtick',xtick1,'ytick',xtick2)
0153 %
0154 figure(3),clf
0155 %
0156 tit = 'p_c/p_T';
0157 %
0158 y = spc;
0159 cont = NaN;%-log10(graph.cont_gamma);
0160 % cont_lab = -log10(graph.cont_gamma_lab);
0161 %
0162 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0163 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0164 %
0165 set(gca,'xtick',xtick1,'ytick',xtick2)
0166 %
0167 %
0168 figure(4),clf
0169 %
0170 tit = 'p_{grid}/p_T';
0171 %
0172 y = spnmax_S;
0173 cont = NaN;%-log10(graph.cont_gamma);
0174 % cont_lab = -log10(graph.cont_gamma_lab);
0175 %
0176 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0177 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0178 %
0179 set(gca,'xtick',xtick1,'ytick',xtick2)
0180 %
0181 %
0182 figure(5),clf
0183 %
0184 tit = 'p_{min}/p_T';
0185 %
0186 y = spnmin;
0187 cont = NaN;%-log10(graph.cont_gamma);
0188 % cont_lab = -log10(graph.cont_gamma_lab);
0189 %
0190 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0191 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0192 %
0193 set(gca,'xtick',xtick1,'ytick',xtick2)
0194 %
0195 figure(6),clf
0196 %
0197 tit = 'p_{max}/p_T';
0198 %
0199 y = spnmax;
0200 cont = NaN;%-log10(graph.cont_gamma);
0201 % cont_lab = -log10(graph.cont_gamma_lab);
0202 %
0203 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0204 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0205 %
0206 set(gca,'xtick',xtick1,'ytick',xtick2)
0207 %
0208 %
0209 figure(7),clf
0210 %
0211 tit = 'E_c (keV)';
0212 %
0213 y = sEcc;
0214 cont = NaN;%-log10(graph.cont_gamma);
0215 % cont_lab = -log10(graph.cont_gamma_lab);
0216 %
0217 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0218 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0219 %
0220 set(gca,'xtick',xtick1,'ytick',xtick2)
0221 %
0222 %
0223 figure(8),clf
0224 %
0225 tit = 'E_{grid} (keV)';
0226 %
0227 y = sEcmax_S;
0228 cont = NaN;%-log10(graph.cont_gamma);
0229 % cont_lab = -log10(graph.cont_gamma_lab);
0230 %
0231 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0232 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0233 %
0234 set(gca,'xtick',xtick1,'ytick',xtick2)
0235 %
0236 %
0237 figure(9),clf
0238 %
0239 tit = 'E_{min} (keV)';
0240 %
0241 y = sEcmin;
0242 cont = NaN;%-log10(graph.cont_gamma);
0243 % cont_lab = -log10(graph.cont_gamma_lab);
0244 %
0245 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0246 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0247 %
0248 set(gca,'xtick',xtick1,'ytick',xtick2)
0249 %
0250 figure(10),clf
0251 %
0252 tit = 'E_{max} (keV)';
0253 %
0254 y = sEcmax;
0255 cont = NaN;%-log10(graph.cont_gamma);
0256 % cont_lab = -log10(graph.cont_gamma_lab);
0257 %
0258 ax = graph2D_jd(x1,x2,y,xlab1,xlab2,tit,xlim1,xlim2,1,NaN,cont,0,'--',NaN,0.5,siz,0.9);colorbar
0259 % hold on,graph2D_jd(x1,x2,y,'','','',xlim1,xlim2,0,NaN,cont_lab,1,'-','k',2,siz,1,NaN,ax);hold off%colorbar
0260 %
0261 set(gca,'xtick',xtick1,'ytick',xtick2)
0262 %
0263 figure(11),clf
0264 %
0265 tit = 'threshold';
0266 %
0267 y2min = NaN(1,nscan1);
0268 y2max = NaN(1,nscan1);
0269 %
0270 for iscan1 = 1:nscan1,
0271     testres = ~isnan(spnmax(iscan1,:));
0272     if any(testres),
0273         y2min(iscan1) = x2(find(testres,1,'first'));
0274         y2max(iscan1) = x2(find(testres,1,'last'));
0275     end
0276 end
0277 %
0278 leg = {'min','max'};
0279 %
0280 graph1D_jd(x1,y2min,xlog.(scanstr1),xlog.(scanstr2),'','','',NaN,xlim.(scanstr1),xlim.(scanstr2),'-','o','r',2,siz);
0281 graph1D_jd(x1,y2max,xlog.(scanstr1),xlog.(scanstr2),scanlab.(scanstr1),scanlab.(scanstr2),tit,leg,xlim.(scanstr1),xlim.(scanstr2),'-','o','r',2,siz,gca,0.9,0.7,0.7);
0282 %
0283 set(gca,'xtick',xtick1,'ytick',xtick2)
0284 %
0285 figure(12),clf
0286 %
0287 y1min = NaN(1,nscan2);
0288 y1max = NaN(1,nscan2);
0289 %
0290 for iscan2 = 1:nscan2,
0291     testres = ~isnan(spnmax(:,iscan2));
0292     if any(testres),
0293         y1min(iscan2) = x1(find(testres,1,'first'));
0294         y1max(iscan2) = x1(find(testres,1,'last'));
0295     end
0296 end
0297 %
0298 graph1D_jd(x2,y1min,xlog.(scanstr2),xlog.(scanstr1),'','','',NaN,xlim.(scanstr2),xlim.(scanstr1),'-','o','r',2,siz);
0299 graph1D_jd(x2,y1max,xlog.(scanstr2),xlog.(scanstr1),scanlab.(scanstr2),scanlab.(scanstr1),tit,leg,xlim.(scanstr2),xlim.(scanstr1),'-','o','r',2,siz,gca,0.9,0.7,0.7);
0300 %
0301 set(gca,'xtick',xtick2,'ytick',xtick1)
0302 %
0303 print_jd(p_opt,['fig_',id_simul_scan,'_nr'],['./figures',filesep,scanstr1,'_',scanstr2],1)
0304 print_jd(p_opt,['fig_',id_simul_scan,'_rr'],['./figures',filesep,scanstr1,'_',scanstr2],2)
0305 print_jd(p_opt,['fig_',id_simul_scan,'_pc'],['./figures',filesep,scanstr1,'_',scanstr2],3)
0306 print_jd(p_opt,['fig_',id_simul_scan,'_pgrid'],['./figures',filesep,scanstr1,'_',scanstr2],4)
0307 print_jd(p_opt,['fig_',id_simul_scan,'_pmin'],['./figures',filesep,scanstr1,'_',scanstr2],5)
0308 print_jd(p_opt,['fig_',id_simul_scan,'_pmax'],['./figures',filesep,scanstr1,'_',scanstr2],6)
0309 print_jd(p_opt,['fig_',id_simul_scan,'_Ec'],['./figures',filesep,scanstr1,'_',scanstr2],7)
0310 print_jd(p_opt,['fig_',id_simul_scan,'_Egrid'],['./figures',filesep,scanstr1,'_',scanstr2],8)
0311 print_jd(p_opt,['fig_',id_simul_scan,'_Emin'],['./figures',filesep,scanstr1,'_',scanstr2],9)
0312 print_jd(p_opt,['fig_',id_simul_scan,'_Emax'],['./figures',filesep,scanstr1,'_',scanstr2],10)
0313 print_jd(p_opt,['fig_',id_simul_scan,'_y2'],['./figures',filesep,scanstr1,'_',scanstr2],11)
0314 print_jd(p_opt,['fig_',id_simul_scan,'_y1'],['./figures',filesep,scanstr1,'_',scanstr2],12)
0315 %
0316 if isfield(iscandisp,scanstr1) && isfield(iscandisp,scanstr2),
0317     %
0318     xlab = 'E_c (MeV)';
0319     ylab = 'f(E_c,\xi = 1)';
0320     xlim = [0,max(sEcmax_S(:))];    
0321     if isfield(graph,'ylim_f'),
0322         ylim = graph.ylim_f;
0323     else
0324         ylim = [eps,2];
0325         graph.ytick_f = 10.^(-18:2:2);
0326     end    
0327     %
0328     nscandisp1 = length(iscandisp.(scanstr1));
0329     nscandisp2 = length(iscandisp.(scanstr2));
0330     %
0331     leg = cell(1,nscandisp2);
0332     for iiscan2 = 1:nscandisp2,
0333         iscan2 = iscandisp.(scanstr2)(iiscan2);
0334         leg{iiscan2} = [scanlab.(scanstr2),' = ',num2str(scanval.(scanstr2)(iscan2))];
0335     end
0336     %
0337     for iiscan1 = 1:nscandisp1,
0338         %
0339         figure(20+iiscan1),clf,
0340         iscan1 = iscandisp.(scanstr1)(iiscan1);
0341         %
0342         for iiscan2 = 1:nscandisp2-1,
0343             iscan2 = iscandisp.(scanstr2)(iiscan2);
0344             graph1D_jd(Ec{iscan1,iscan2},sfpar{iscan1,iscan2},0,1,'','','',NaN,xlim,ylim,'-','none',colors{iiscan2},2,siz);
0345         end
0346         iscan2 = iscandisp.(scanstr2)(nscandisp2);
0347         graph1D_jd(Ec{iscan1,iscan2},sfpar{iscan1,iscan2},0,1,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{nscandisp2},2,siz,gca,0.9,0.7,0.7);
0348         %
0349         if isfield(graph,'ytick_f'),
0350             set(gca,'ytick',graph.ytick_f)
0351         end
0352         set(gca,'YMinorGrid','off')
0353         set(gca,'YMinorTick','on')
0354         %
0355         print_jd(p_opt,['fig_',id_simul_scan,'_f_',num2str(iiscan1)],['./figures',filesep,scanstr1,'_',scanstr2],20+iiscan1)
0356         %
0357     end
0358     %
0359 end
0360 %