0001 function plot_fluct_jd(itn,xrho,pn,txJ,txP,XXfM,tXXf0,simul,Te0,opt)
0002
0003 if 0,
0004
0005 nt = length(itn);
0006 npn = length(pn);
0007
0008 ppar = [fliplr(-pn(2:end)),pn(1:end)];
0009 tfpar = NaN(nt,2*npn - 1);
0010
0011 ir_P = find(sum(txP,1) == max(sum(txP,1)));
0012
0013 wavestr = [simul.path,'wave_results/waves_',simul.id];
0014
0015 data = load([wavestr,'_1.mat'],'waves');
0016
0017 iy = 0;
0018 for iw = 1:length(data.waves),
0019 iiy = 0;
0020 for ib = 1:length(data.waves{iw}.rayinit.launch.bNpar0),
0021 for ir = 1:length(data.waves{iw}.rayinit.launch.rZ0),
0022
0023 iy = iy + 1;
0024 iiy = iiy + 1;
0025
0026 yiw(iy) = iw;
0027 yNpar0(iy) = data.waves{iw}.rayinit.launch.bNpar0(ib);
0028 ydNpar0(iy) = data.waves{iw}.rayinit.launch.bdNpar0(ib);
0029 ytheta0(iy) = data.waves{iw}.rayinit.ytheta0(iiy);
0030
0031 if data.waves{iw}.rayinit.launch.tail.mode == 3,
0032
0033 if data.waves{iw}.rayinit.launch.tail.bNparmax_tail == 0,
0034 Nparmax_tail = 6.5/sqrt(Te0);
0035 else
0036 Nparmax_tail = data.waves{iw}.rayinit.launch.tail.Nparmax_tail;
0037 end
0038
0039 [tail(iy).Npar0_tail,tail(iy).dNpar0_tail,tail(iy).P0_2piRp_tail] = calc_tail_jd(yNpar0(iy),ydNpar0(iy),1,...
0040 Nparmax_tail,data.waves{iw}.rayinit.launch.tail.bn_tail,data.waves{iw}.rayinit.launch.tail.bP_tail,data.waves{iw}.rayinit.launch.tail.bopt_tail,0);
0041
0042 end
0043
0044 end
0045 end
0046 end
0047
0048 ny = iy;
0049
0050 yNpar = NaN(nt,ny);
0051 ydNpar = NaN(nt,ny);
0052 yprop = false(nt,ny);
0053
0054 fparM = [flipud(XXfM(2:end,1,ir_P));XXfM(1:end,end,ir_P)].';
0055
0056 for it = 1:nt,
0057
0058 tfpar(it,:) = [flipud(tXXf0{it}(2:end,1,ir_P));tXXf0{it}(1:end,end,ir_P)].';
0059
0060 data = load([wavestr,'_',num2str(it),'.mat'],'waves');
0061
0062 iy = 0;
0063 for iw = 1:length(data.waves),
0064 for iiy = 1:length(data.waves{iw}.rayinit.yNpar0),
0065
0066 iy = iy + 1;
0067
0068 yNpar(it,iy) = data.waves{iw}.rayinit.yNpar0(iiy);
0069 ydNpar(it,iy) = data.waves{iw}.rayinit.ydNpar0(iiy);
0070
0071 end
0072
0073 for iray = 1:length(data.waves{iw}.rays),
0074
0075 iy = find(yiw == iw & yNpar(it,:) == data.waves{iw}.rays{iray}.rayinits{1}.yNpar0 & ytheta0 == data.waves{iw}.rays{iray}.rayinits{1}.ytheta0);
0076
0077 if isempty(iy) || length(iy) > 1,
0078 disp('warning : pb in ray identification.')
0079 keyboard
0080 else
0081 yprop(it,iy) = true;
0082 end
0083
0084 end
0085 end
0086
0087 disp(['time index ',num2str(it),'/',num2str(nt),' processed.'])
0088
0089 end
0090
0091 clear tXXf0 data
0092
0093 else
0094
0095 load toto2.mat
0096 keyboard
0097
0098 end
0099
0100 colors = repmat({'r';'b';'k';[0,0.5,0]},[1,4]);
0101 markers = repmat({'+','o','s','d'},[4,1]);
0102
0103 yNpar = abs(yNpar);
0104
0105
0106 it = 1;
0107
0108 nx = 101;
0109 nbins = 10;
0110 pmax = 15;
0111 dit = 10;
0112
0113 fac0 = 1/3;
0114 res = 600;
0115
0116 fig.position_slider = [0,0,1200,400];
0117 ax1.position_slider = [30,100,320,270];
0118 ax2.position_slider = [440,100,320,270];
0119 ax3.position_slider = [850,100,320,270];
0120
0121 fig.position_movie = [0,0,1200,360];
0122 ax1.position_movie = [30,60,320,270];
0123 ax2.position_movie = [440,60,320,270];
0124 ax3.position_movie = [850,60,320,270];
0125
0126 slider.position = [150,10,1000,30];
0127 txt.position = [30,10,40,30];
0128 itxt.position = [80,10,40,30];
0129
0130 slider.callback = @slider_callback;
0131 itxt.callback = @txt_callback;
0132
0133 if opt.mov == 0,
0134 fig.handle = figure('name','fluctuations','visible','on','menubar','figure','numbertitle','off','resize','off','position',fig.position_slider,'toolbar','none','units','pixels');
0135 else
0136 fig.handle = figure('name','fluctuations','visible','on','menubar','figure','numbertitle','off','resize','off','position',fig.position_movie,'toolbar','none','units','pixels');
0137 end
0138
0139 subplot(1,3,1);
0140
0141 xlim = [0,1];
0142 xlab = 'r/a';
0143 ylim = 1.1*[min([txJ(:);txP(:)]),max([txJ(:);txP(:)])];
0144 ylab = '';
0145 tit = '';
0146 leg_JP = {'J (MA/m^2)','P (MW/m^3)'};
0147
0148 graph1D_jd(xrho,txJ(it,:),0,0,'','','',NaN,xlim,ylim,'-','+','r',2,20+14i);
0149 graph1D_jd(xrho,txP(it,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','o','b',2,20+14i);
0150
0151 subplot(1,3,2);
0152
0153 xlim = [-pmax,pmax];
0154 xlab = 'p/p_T';
0155 ylim = [1e-8,1];
0156 ylab = 'f(p) [p_{\perp} = 0]';
0157 tit = '';
0158 leg_f = {'f_M','f_{LH}'};
0159
0160 graph1D_jd(ppar,fparM,0,1,'','','',NaN,xlim,ylim,'-','none','b',0.5,20+14i);
0161 graph1D_jd(ppar,tfpar(it,:),0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,20+14i);
0162
0163 subplot(1,3,3);
0164
0165 if ~isscalar(opt.iy) || ~any(opt.iy == 1:ny) || isempty(tail(opt.iy).Npar0_tail),
0166 opt.iy = 0;
0167 end
0168
0169 if opt.iy == 0,
0170
0171 xlim = [0,nt+1];
0172 xlab = 'time index';
0173 ylim = NaN;
0174 ylab = 'N_{||0}';
0175 tit = '';
0176
0177 yNpar_true = yNpar;
0178 yNpar_false = yNpar;
0179
0180 yNpar_true(~yprop) = NaN;
0181 yNpar_false(yprop) = NaN;
0182
0183 for iy = 1:ny,
0184 graph1D_jd(1:nt,yNpar_true(:,iy),0,0,'','','',NaN,xlim,ylim,'none',markers{iy},colors{iy},2,20+14i);
0185 graph1D_jd(1:nt,yNpar_false(:,iy),0,0,'','','',NaN,xlim,ylim,'none',markers{iy},colors{iy},0.5,20+14i);
0186 end
0187
0188 ylim = get(gca,'ylim');
0189
0190 graph1D_jd([it,it],ylim,0,0,xlab,ylab,tit,NaN,xlim,ylim,'--','none','k',2,20+14i);
0191
0192 else
0193
0194 [xN,xdPdN,xlim] = tail_spec(tail(opt.iy),nx);
0195 xdPdN0 = spec0(yNpar(it,opt.iy),ydNpar(it,opt.iy),xN);
0196
0197 leg_Npar = {'Hist.','Fluct.','Tail'};
0198
0199
0200
0201
0202 xdPdN = xdPdN*diff(xlim)/nbins;
0203 xdPdN0 = xdPdN0*diff(xlim)/nbins*fac0;
0204
0205 edges = [linspace(xlim(1),xlim(2),nbins+1)];
0206
0207
0208
0209 counts = histc(yNpar(1:it,opt.iy),edges)/it;
0210
0211 bar(edges,counts,'histc');
0212
0213 xlab = 'N_{||}';
0214 ylim = NaN;
0215 ylab = 'dPdN_{||} (norm.)';
0216 tit = '';
0217
0218 graph1D_jd(xN,xdPdN0,0,0,'','','',NaN,xlim,ylim,'-','none',[0,0.5,0],2,20+14i);
0219 graph1D_jd(xN,xdPdN,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,20+14i);
0220
0221 end
0222
0223 axes = get(gcf,'children');
0224
0225 ax1.handle = axes(3);
0226 ax2.handle = axes(2);
0227 ax3.handle = axes(1);
0228
0229 if opt.mov == 0,
0230
0231 set(ax1.handle,'units','pixels','position',ax1.position_slider)
0232 set(ax2.handle,'units','pixels','position',ax2.position_slider)
0233 set(ax3.handle,'units','pixels','position',ax3.position_slider)
0234
0235 legend(ax1.handle,leg_JP)
0236 legend(ax2.handle,leg_f)
0237
0238 if opt.iy,
0239 legend(ax3.handle,leg_Npar)
0240 Npar = yNpar(:,opt.iy);
0241 dNpar = ydNpar(:,opt.iy);
0242 else
0243 Npar = NaN;
0244 dNpar = NaN;
0245 end
0246
0247 slider.handle = uicontrol(fig.handle,'Style','slider','Position',slider.position,'min',1,'max',nt,'sliderstep',[1/(nt-1),10/(nt-1)],'Value',it);
0248 txt.handle = uicontrol(fig.handle,'Style','text','Position',txt.position,'string',['# = '],'FontSize',20);
0249 itxt.handle = uicontrol(fig.handle,'callback',itxt.callback,'Style','edit','Position',itxt.position,'string',num2str(it),'FontSize',20);
0250
0251 set(slider.handle,'callback',{slider.callback,itxt.handle,ax1.handle,txJ,txP,ax2.handle,tfpar,ax3.handle,Npar,dNpar,fac0});
0252 set(itxt.handle,'callback',{itxt.callback,slider.handle,ax1.handle,txJ,txP,ax2.handle,tfpar,ax3.handle,Npar,dNpar,fac0});
0253
0254 else
0255
0256 set(ax1.handle,'units','pixels','position',ax1.position_movie)
0257 set(ax2.handle,'units','pixels','position',ax2.position_movie)
0258 set(ax3.handle,'units','pixels','position',ax3.position_movie)
0259
0260 legend(ax1.handle,leg_JP)
0261 legend(ax2.handle,leg_f)
0262
0263 if opt.iy,
0264 legend(ax3.handle,leg_Npar)
0265 Npar = yNpar(:,opt.iy);
0266 dNpar = ydNpar(:,opt.iy);
0267 else
0268 Npar = NaN;
0269 dNpar = NaN;
0270 end
0271
0272 iit = 0;
0273
0274 for it = 1:dit:nt,
0275
0276 iit = iit + 1;
0277
0278 update_plot(it,ax1.handle,txJ,txP,ax2.handle,tfpar,ax3.handle,Npar,dNpar,fac0);
0279
0280
0281
0282 savename = ['Fig_fluct_frame_',repmat('0',[1,4-ceil(log10(iit+1))]),num2str(iit)];
0283 eval(['print -dpng -r',num2str(res),' ',savename,'.png']);
0284
0285 end
0286
0287 end
0288
0289 if opt.spec,
0290
0291 for iy = 1:ny,
0292
0293 figure(100+iy),clf
0294
0295 if ~isempty(tail(iy).Npar0_tail),
0296
0297 [xN,xdPdN,xlim] = tail_spec(tail(iy),nx);
0298
0299 leg = {'Fluct.','Tail'};
0300
0301 xdPdN = xdPdN*nt*diff(xlim)/nbins;
0302
0303 edges = [linspace(xlim(1),xlim(2),nbins+1)];
0304 counts = histc(yNpar(:,iy),edges);
0305 bar(edges,counts,'histc');
0306
0307 graph1D_jd(xN,xdPdN,0,0,'N_{||}','dPdN_{||} (norm.)','',leg,xlim,NaN,'-','none','r',2,20+14i,gca,0.9,0.7,0.7);
0308
0309 end
0310
0311 end
0312
0313 end
0314
0315 if isfield(opt,'itcomp') && length(opt.itcomp) > 1,
0316
0317 figure(1001),clf
0318
0319 locxlim = [0,1];
0320 xlab = 'r/a';
0321 ylim = 1.1*[min(min([txJ(opt.itcomp,:);txP(opt.itcomp,:)])),max(max([txJ(opt.itcomp,:);txP(opt.itcomp,:)]))];
0322 ylab = '';
0323 tit = '';
0324
0325 nitcomp = length(opt.itcomp);
0326 locleg = cell(1,2*nitcomp);
0327 for locit = 1:nitcomp-1,
0328 locleg{1+2*(locit - 1)} = ['it. # ',num2str(opt.itcomp(locit)),': J (MA/m^2)'];
0329 locleg{2+2*(locit - 1)} = ['it. # ',num2str(opt.itcomp(locit)),': P (MW/m^3)'];
0330 graph1D_jd(xrho,txJ(opt.itcomp(locit),:),0,0,'','','',NaN,locxlim,ylim,'-','+',colors{locit},2,20+14i);
0331 graph1D_jd(xrho,txP(opt.itcomp(locit),:),0,0,'','','',NaN,locxlim,ylim,'-','o',colors{locit},2,20+14i);
0332 end
0333 locleg{1+2*(nitcomp - 1)} = ['it. # ',num2str(opt.itcomp(nitcomp)),': J (MA/m^2)'];
0334 locleg{2+2*(nitcomp - 1)} = ['it. # ',num2str(opt.itcomp(nitcomp)),': P (MW/m^3)'];
0335 graph1D_jd(xrho,txJ(opt.itcomp(nitcomp),:),0,0,'','','',NaN,locxlim,ylim,'-','+',colors{nitcomp},2,20+14i);
0336 graph1D_jd(xrho,txP(opt.itcomp(nitcomp),:),0,0,xlab,ylab,tit,locleg,locxlim,ylim,'-','o',colors{nitcomp},2,20+14i,gca,0.9,0.7,0.7);
0337
0338
0339 if opt.iy > 0,
0340
0341 xdPdN0comp = NaN(nitcomp,nx);
0342 for locit = 1:nitcomp,
0343 xdPdN0comp(locit,:) = spec0(yNpar(opt.itcomp(locit),opt.iy),ydNpar(opt.itcomp(locit),opt.iy),xN);
0344 end
0345 xdPdN0comp = xdPdN0comp*diff(xlim)/nbins*fac0;
0346 counts = histc(yNpar(:,opt.iy),edges)/nt;
0347
0348 figure(1002),clf
0349
0350 xlab = 'N_{||}';
0351 ylim = NaN;
0352 ylab = 'dPdN_{||} (norm.)';
0353 tit = '';
0354
0355 locleg = cell(1,nitcomp);
0356 for locit = 1:nitcomp-1,
0357 locleg{locit} = ['it. # ',num2str(opt.itcomp(locit))];
0358 graph1D_jd(xN,xdPdN0comp(locit,:),0,0,'','','',NaN,xlim,ylim,'-','none',colors{locit},2,20+14i);
0359 end
0360 locleg{nitcomp} = ['it. # ',num2str(opt.itcomp(nitcomp))];
0361 graph1D_jd(xN,xdPdN0comp(nitcomp,:),0,0,xlab,ylab,tit,locleg,xlim,ylim,'-','none',colors{nitcomp},2,20+14i,gca,0.9,0.7,0.7);
0362
0363 hold on
0364 bar(edges,counts,0.6,'histc');
0365 hold off
0366
0367 figure(1003),clf
0368
0369 xlab = 'N_{||}';
0370 ylim = NaN;
0371 ylab = 'counts';
0372 tit = '';
0373
0374 graph1D_jd(NaN,NaN,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,20+14i,gca,0.9,0.7,0.7);
0375
0376 hold on
0377 bar(edges,counts*nt,0.6,'histc');
0378 hold off
0379
0380 end
0381
0382 end
0383
0384 end
0385
0386 function slider_callback(handle,eventdata,itxt_handle,ax1_handle,txJ,txP,ax2_handle,tfpar,ax3_handle,Npar,dNpar,fac0)
0387
0388 value = get(handle,'value');
0389
0390
0391
0392 it = round(value);
0393
0394 set(itxt_handle,'string',num2str(it));
0395
0396 update_plot(it,ax1_handle,txJ,txP,ax2_handle,tfpar,ax3_handle,Npar,dNpar);
0397
0398 end
0399
0400 function txt_callback(handle,eventdata,slider_handle,ax1_handle,txJ,txP,ax2_handle,tfpar,ax3_handle,Npar,dNpar,fac0)
0401
0402 txt = get(handle,'string');
0403
0404
0405
0406 nt = get(slider_handle,'max');
0407
0408 it = round(str2double(txt));
0409 it(it < 1) = 1;
0410 it(it > nt) = nt;
0411
0412 set(handle,'string',num2str(it));
0413 set(slider_handle,'Value',it);
0414
0415 update_plot(it,ax1_handle,txJ,txP,ax2_handle,tfpar,ax3_handle,Npar,dNpar);
0416
0417 end
0418
0419 function update_plot(it,ax1_handle,txJ,txP,ax2_handle,tfpar,ax3_handle,Npar,dNpar,fac0)
0420
0421 lines = get(ax1_handle,'children');
0422 set(lines(2),'ydata',txJ(it,:));
0423 set(lines(1),'ydata',txP(it,:));
0424
0425 lines = get(ax2_handle,'children');
0426 set(lines(1),'ydata',tfpar(it,:));
0427
0428 lines = get(ax3_handle,'children');
0429
0430 if length(get(lines(1),'xdata')) == 2,
0431 set(lines(1),'xdata',[it,it]);
0432 else
0433
0434 xlim = get(ax3_handle,'xlim');
0435
0436 nbins = size(get(lines(3),'xdata'),2) - 1;
0437
0438
0439 edges = [linspace(xlim(1),xlim(2),nbins+1)];
0440
0441
0442
0443 counts = histc(Npar(1:it),edges)/it;
0444
0445 delete(lines(3)),
0446 hold on
0447 bar(ax3_handle,edges,counts,'histc');
0448 hold off
0449
0450 xN = get(lines(2),'xdata');
0451
0452 xdPdN0 = spec0(Npar(it),dNpar(it),xN);
0453
0454
0455
0456 xdPdN0 = xdPdN0*diff(xlim)/nbins*fac0;
0457
0458 set(lines(2),'ydata',xdPdN0);
0459
0460
0461 lines = get(ax3_handle,'children');
0462 set(ax3_handle,'children',lines([2,3,1]));
0463
0464 end
0465
0466 end
0467
0468 function [xN,xdPdN,xlim] = tail_spec(tail,nx)
0469
0470 Npar0_tail = tail.Npar0_tail;
0471 dNpar0_tail = tail.dNpar0_tail;
0472 P0_2piRp_tail = tail.P0_2piRp_tail;
0473
0474 n_tail = length(Npar0_tail);
0475
0476 Npar0_tail = abs(Npar0_tail);
0477
0478
0479 xlim = [min(Npar0_tail) - 2*max(dNpar0_tail),max(Npar0_tail) + 2*max(dNpar0_tail)];
0480
0481 xN = linspace(xlim(1),xlim(2),nx);
0482
0483 xdPdN = zeros(1,nx);
0484 for in = 1:n_tail,
0485 xdPdN = xdPdN + P0_2piRp_tail(in)*exp(-(xN - Npar0_tail(in)).^2/dNpar0_tail(in).^2)/(sqrt(pi)*dNpar0_tail(in));
0486 end
0487
0488 end
0489
0490 function xdPdN = spec0(Npar,dNpar,xN)
0491
0492 xdPdN = exp(-(xN - Npar).^2/dNpar.^2)/(sqrt(pi)*dNpar);
0493
0494 end
0495
0496
0497
0498
0499
0500
0501