fproc_varscan

0001 function fproc_varscan(varargin)
0002 %
0003 if nargin == 3,
0004     [simulist,var,opt] = varargin{:};
0005 elseif nargin == 9,
0006     [basestr,file_str,shotnum,shotime,id_simul_o,locid_simul,locpath_simul,var,opt] = varargin{:};
0007 else
0008     error('fproc_varscan takes either 3 or 9 arguments')
0009 end
0010 %
0011 if nargin == 9 && ischar(locpath_simul),
0012     path_simul = [basestr,filesep,shotnum,filesep,locpath_simul];
0013 elseif nargin == 9 && iscell(locpath_simul),
0014     path_simul = locpath_simul{1};
0015 else
0016     path_simul = '';
0017 end
0018 %
0019 % select shot times if time scan
0020 %
0021 if nargin == 3,%simulation files directly given
0022     %
0023     flag = false(1,length(simulist));
0024     %
0025     for ivar = 1:length(simulist),
0026         if exist(simulist{ivar},'file'),
0027             flag(ivar) = true;
0028         end
0029     end
0030     %
0031     simulist = simulist(flag);
0032     nvar = length(simulist);
0033     %
0034     ivar_disp = [];
0035     for ivar = 1:nvar,
0036         var.liststr{ivar} = num2str(var.list(ivar));
0037         if any(abs(var.list(ivar) - var.disp) < 10*eps),
0038             ivar_disp = [ivar_disp,ivar];
0039         end
0040     end
0041     %
0042     clear basestr
0043     %
0044     figstr = 'fig';
0045     figdir = './';
0046     %
0047 elseif iscell(shotime),
0048     equilbase = file_str;
0049     id_simul_o = basestr;
0050     var.str = shotnum;
0051     %
0052     var.liststr = select_shotimes_jd([equilbase,basestr,'_',shotnum,'_'],[locid_simul,'.mat'],path_simul,shotime);
0053     nvar = length(var.liststr);
0054     %
0055     ivar_disp = [];
0056     for ivar = 1:nvar,
0057         ttimes(ivar) = str2num(var.liststr{ivar});
0058         if any(ttimes(ivar) == var.disp),
0059             ivar_disp = [ivar_disp,ivar];
0060         end
0061     end
0062     %
0063     if opt.xdisp == 0 || ~strcmp(basestr,'TS'),
0064         var.list = ttimes;
0065     else
0066         [nvol,tn]=tsbase(str2double(shotnum),'sprofnmoy');
0067         if isempty(nvol),
0068             [nvol,tn]=tsbase(str2double(shotnum),'snmoy');
0069             nvol = nvol/1e19;
0070         end
0071         var.list = interp1(tn,nvol,ttimes);
0072         var.lab = '<n>';
0073         var.unit = '\times 10^{19} m^{-3}';
0074     end
0075     %
0076     figstr = [equilbase,id_simul_o,'_',var.str,'_scan',locid_simul];
0077     figdir = path_simul;
0078     %
0079 else
0080     equilbase = [file_str,basestr,'_',shotnum,'_',shotime];
0081     nvar = length(var.list);
0082     ttimes = repmat(str2num(shotime),[1,nvar]);
0083     ivar_disp = [];
0084     for ivar = 1:nvar,
0085         var.liststr{ivar} = num2str(var.list(ivar));
0086         if any(abs(var.list(ivar) - var.disp) < 10*eps),
0087             ivar_disp = [ivar_disp,ivar];
0088         end
0089     end
0090     %
0091     figstr = [equilbase,id_simul_o,'_',var.str,'_scan',locid_simul];
0092     figdir = path_simul;
0093     %
0094 end
0095 %
0096 tol_P0 = 1;%0.04;%tolerance on absorbed power verification
0097 %
0098 %************************************************************************************************************************************
0099 %
0100 % Load structures
0101 %
0102 %[dkepath] = load_structures_yp('dkepath',id_dkepath,path_dkepath);
0103 %
0104 dRp = NaN(1,nvar);
0105 dI_rf = NaN(1,nvar);
0106 dI_exp = NaN(1,nvar);
0107 dP_rf = NaN(1,nvar);
0108 dP_lin_rf = NaN(1,nvar);
0109 dP0_rf = NaN(1,nvar);
0110 drho_J = NaN(1,nvar);
0111 ddrho_J = NaN(1,nvar);
0112 drho_P = NaN(1,nvar);
0113 dTe = NaN(1,nvar);
0114 dne = NaN(1,nvar);
0115 %
0116 ds_c3po = NaN(1,nvar);
0117 ds_luke = NaN(1,nvar);
0118 %
0119 dNpar0 = NaN(1,nvar);
0120 dNpar_c3po = NaN(1,nvar);
0121 dNpar_luke = NaN(1,nvar);
0122 %
0123 dxrho = cell(1,nvar);
0124 dxJ = cell(1,nvar);
0125 dxP = cell(1,nvar);
0126 dxdpdt_rf = cell(1,nvar);
0127 %
0128 drho_P_lin = NaN(1,nvar);
0129 dTe_lin = NaN(1,nvar);
0130 dne_lin = NaN(1,nvar);
0131 %
0132 if ~isfield(opt,'enforce'),
0133     opt.enforce = 0;
0134 end
0135 %
0136 for ivar = 1:nvar,
0137     %
0138     if nargin == 3,
0139         [datapath,id_simul] = fileparts(simulist{ivar});
0140         datapath = [datapath,filesep];
0141     else
0142         datapath = path_simul;        
0143         id_simul = [equilbase,id_simul_o,'_',var.str,'_',var.liststr{ivar},locid_simul];%Simulation ID
0144     end
0145     %
0146     data_proc = fproc_data(datapath,id_simul,opt);
0147     %
0148     if isempty(data_proc),
0149         continue
0150     end
0151     %
0152     dP_rf(ivar) = data_proc.scalar.p_rf_2piRp(1)*2*pi*data_proc.scalar.Rp;
0153     dP_lin_rf(ivar) = data_proc.scalar.p_rf_2piRp_lin(1)*2*pi*data_proc.scalar.Rp;
0154     dP0_rf(ivar) = data_proc.scalar.p0_rf_2piRp(1)*2*pi*data_proc.scalar.Rp;
0155     %
0156     if ~opt.enforce && data_proc.simul.residue_rf > 1e-2
0157         disp(['Case ',var.lab,' = ',var.liststr{ivar},' : LUKE run did not converge or numerical instability developed.'])
0158     elseif ~opt.enforce && abs(dP_rf(ivar) - dP0_rf(ivar))/abs(dP0_rf(ivar)) > tol_P0,
0159         disp(['Case ',var.lab,' = ',var.liststr{ivar},' : verify wave power absorption'])
0160         keyboard
0161     elseif ~opt.enforce && isfield(opt,'P0_min') && dP0_rf(ivar) < opt.P0_min,
0162         disp(['Case ',var.lab,' = ',var.liststr{ivar},' : some rays did not enter the plasma.'])
0163     else
0164         %
0165         dRp(ivar) = data_proc.scalar.Rp;
0166         dI_rf(ivar) = data_proc.scalar.I_tot;
0167         %
0168         if isfield(data_proc.scalar,'I_exp'),
0169             dI_exp(ivar) = data_proc.scalar.I_exp;
0170         end
0171         %
0172         if opt.rho == 'g',
0173             drho_J(ivar) = data_proc.scalar.jrhoG(1);
0174             ddrho_J(ivar) = data_proc.scalar.jdrhoG(1);
0175             if isfield(data_proc.wave,'wrhoG') && ~isempty(data_proc.wave.wrhoG),
0176                 drho_P(ivar) = data_proc.wave.wrhoG(1);
0177             end
0178             if isfield(data_proc.wave,'wrhoG_lin') && ~isempty(data_proc.wave.wrhoG_lin),
0179                 drho_P_lin(ivar) = data_proc.wave.wrhoG_lin(1);
0180             end
0181             dxrho{ivar} = data_proc.radial.xrhoG;
0182         elseif opt.rho == 'p',
0183             drho_J(ivar) = data_proc.scalar.jrhoP(1);
0184             ddrho_J(ivar) = data_proc.scalar.jdrhoP(1);
0185             if isfield(data_proc.wave,'wrhoP') && ~isempty(data_proc.wave.wrhoP),
0186                 drho_P(ivar) = data_proc.wave.wrhoP(1);
0187             end
0188             if isfield(data_proc.wave,'wrhoP_lin') && ~isempty(data_proc.wave.wrhoP_lin),
0189                 drho_P_lin(ivar) = data_proc.wave.wrhoP_lin(1);
0190             end
0191             dxrho{ivar} = data_proc.radial.xrhoP;
0192         elseif opt.rho == 't',
0193             drho_J(ivar) = data_proc.scalar.jrhoT(1);
0194             ddrho_J(ivar) = data_proc.scalar.jdrhoT(1);
0195             if isfield(data_proc.wave,'wrhoT') && ~isempty(data_proc.wave.wrhoT),
0196                 drho_P(ivar) = data_proc.wave.wrhoT(1);
0197             end
0198             if isfield(data_proc.wave,'wrhoT_lin') && ~isempty(data_proc.wave.wrhoT_lin),
0199                 drho_P_lin(ivar) = data_proc.wave.wrhoT_lin(1);
0200             end
0201             dxrho{ivar} = data_proc.radial.xrhoT;
0202         elseif opt.rho == 'v',
0203             drho_J(ivar) = data_proc.scalar.jrhoV(1);
0204             ddrho_J(ivar) = data_proc.scalar.wrhoV(1);
0205             if isfield(data_proc.wave,'wrhoV') && ~isempty(data_proc.wave.wrhoV),
0206                 drho_P(ivar) = data_proc.wave.wrhoG(1);
0207             end
0208             if isfield(data_proc.wave,'wrhoV_lin') && ~isempty(data_proc.wave.wrhoV_lin),
0209                 drho_P_lin(ivar) = data_proc.wave.wrhoV_lin(1);
0210             end
0211             dxrho{ivar} = data_proc.radial.xrhoV;
0212         else
0213             error('rho type not recognized')
0214         end
0215         %
0216         if isfield(data_proc.wave,'wTe') && ~isempty(data_proc.wave.wTe),
0217             dTe(ivar) = data_proc.wave.wTe(1);
0218         end
0219         if isfield(data_proc.wave,'wne') && ~isempty(data_proc.wave.wne),
0220             dne(ivar) = data_proc.wave.wne(1);
0221         end
0222         %
0223         if isfield(data_proc.wave,'wTe_lin') && ~isempty(data_proc.wave.wTe_lin),
0224             dTe_lin(ivar) = data_proc.wave.wTe_lin(1);
0225         end
0226         if isfield(data_proc.wave,'wne_lin') && ~isempty(data_proc.wave.wne_lin),
0227             dne_lin(ivar) = data_proc.wave.wne_lin(1);
0228         end
0229         %
0230         dxJ{ivar} = data_proc.radial.J_tot;
0231         dxP{ivar} = data_proc.radial.P_rf;
0232         %
0233         if isfield(data_proc.radial,'xdpdt_rf'),
0234             dxdpdt_rf{ivar} = data_proc.radial.xdpdt_rf;
0235         end
0236         %
0237         ds_c3po(ivar) = 0;
0238         ds_luke(ivar) = 0;
0239         %
0240         if strcmp(var.str,'bNpar0'),
0241             dNpar0(ivar) = var.list(ivar);
0242         elseif isfield(opt,'Npar0'),
0243             dNpar0(ivar) = opt.Npar0;
0244         end
0245         %
0246         dNpar_c3po(ivar) = 0;
0247         dNpar_luke(ivar) = 0;
0248         %
0249         if isfield(data_proc,'wave') && isfield(data_proc.wave,'wrays'),
0250             xs_c3po = [];
0251             xNpar_c3po = [];
0252             xP0_2piRp_c3po = [];
0253             xs_luke = [];
0254             xNpar_luke = [];
0255             xP0_2piRp_luke = [];
0256             for iw = 1:length(data_proc.wave.wrays),
0257                 for iray = 1:length(data_proc.wave.wrays{iw}),
0258                     xs_c3po = [xs_c3po,data_proc.wave.wrays{iw}{iray}.c3po.ss];
0259                     xP0_2piRp_c3po = [xP0_2piRp_c3po,data_proc.wave.wrays{iw}{iray}.c3po.P0_2piRp];
0260                     xNpar_c3po = [xNpar_c3po,abs(real(data_proc.wave.wrays{iw}{iray}.c3po.sNpar))];
0261                     if isfield(data_proc.wave.wrays{iw}{iray},'luke'),
0262                         xs_luke = [xs_luke,data_proc.wave.wrays{iw}{iray}.luke{end}.ss];
0263                         xNpar_luke = [xNpar_luke,abs(real(data_proc.wave.wrays{iw}{iray}.luke{end}.sNpar))];
0264                         xP0_2piRp_luke = [xP0_2piRp_luke,data_proc.wave.wrays{iw}{iray}.luke{end}.P0_2piRp];
0265                     end
0266                 end
0267             end
0268             %
0269             ds_c3po(ivar) = sum(xs_c3po.*xP0_2piRp_c3po)/sum(xP0_2piRp_c3po);
0270             dNpar_c3po(ivar) = sum(xNpar_c3po.*xP0_2piRp_c3po)/sum(xP0_2piRp_c3po);
0271             %
0272             if ~isempty(xs_luke),
0273                 ds_luke(ivar) = sum(xs_luke.*xP0_2piRp_luke)/sum(xP0_2piRp_luke);
0274                 dNpar_luke(ivar) = sum(xNpar_luke.*xP0_2piRp_luke)/sum(xP0_2piRp_luke);
0275             end
0276         end
0277         %
0278     end
0279 end
0280 %
0281 deta = (dne/1e19).*dI_rf.*dRp./dP_rf;
0282 dksirf = 3.27*deta./dTe;
0283 %
0284 % display results
0285 %
0286 if isfield(opt,'mom') && opt.mom == 'j',
0287     drho = drho_J;
0288     ddrho = ddrho_J;
0289 else
0290     drho = drho_P;
0291     ddrho = NaN(1,nvar);
0292 end
0293 %
0294 style0 = 'none';
0295 style = '-';
0296 style2 = '--';
0297 marker = 'none';
0298 marker1 = '+';
0299 marker2 = 'o';
0300 marker3 = 's';
0301 marker4 = '.';
0302 color1 = 'k';
0303 color2 = 'r';
0304 color3 = 'b';
0305 color4 = 'g';
0306 color5 = 'm';
0307 color6 = 'c';
0308 %
0309 colors = {color1,color2,color3,color4};
0310 markers = {marker1,marker2,marker3,marker4};
0311 %
0312 width = 0.5;
0313 width2 = 2;
0314 siz = 20+14i;
0315 %
0316 red = 0.9;
0317 lspace = 0.5;
0318 lspace2 = 0.7;
0319 bspace = 0.5;
0320 bspace2 = 0.6;
0321 %
0322 if ~isempty(var.unit),
0323     xlab = [var.lab,' (',var.unit,')'];
0324 else
0325     xlab = var.lab;
0326 end
0327 %
0328 if isfield(opt,'xlim'),
0329     xlim = opt.xlim;
0330 else
0331     xlim = [min(var.list),max(var.list)];
0332 end
0333 %
0334 if isfield(opt,'xtick'),
0335     xtick = opt.xtick;
0336 else
0337     xtick = var.list;
0338 end
0339 %
0340 tit0 = '';%[basestr,' #',shotnum];
0341 %
0342 figure(1),clf
0343 %
0344 ylab = '';
0345 %
0346 if isfield(opt,'ylim_Irho'),
0347     ylim = opt.ylim_Irho;
0348 else
0349     ylim = [0,1.2];
0350 end
0351 %
0352 if isfield(opt,'Ifac') && opt.Ifac > 1,
0353     Ifacstr = [' \times',num2str(opt.Ifac)];
0354 else
0355     opt.Ifac = 1;
0356     Ifacstr = '';
0357 end
0358 %
0359 leg = {['I_{LH}',Ifacstr,' (MA)'],'\rho_{dep}','Location','NorthWest'};
0360 tit = '';
0361 %
0362 graph1D_jd(var.list,dI_rf/opt.Ifac,0,0,'','','',NaN,xlim,ylim,style,marker1,color2,width2,siz);
0363 graph1D_jd(var.list,drho,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker2,color3,width2,siz,gca,red,lspace2,bspace2);
0364 set(gca,'xtick',xtick);
0365 if isfield(opt,'ytick_Irho') && ~any(isnan(opt.ytick_Irho)),
0366     set(gca,'ytick',opt.ytick_Irho);
0367 end
0368 %
0369 figure(101),clf
0370 %
0371 if isfield(opt,'ylim_I'),
0372     ylim = opt.ylim_I;
0373 else
0374     ylim = 1e3*[0,max(dI_rf)]*1.2;
0375 end
0376 %
0377 ylab = 'I_{LH} (kA)';
0378 tit = '';
0379 %
0380 if isfield(opt,'I_exp') && opt.I_exp == 1,
0381     graph1D_jd(var.list,1e3*dI_exp/opt.Ifac,0,0,'','','',NaN,xlim,ylim,style,marker2,color3,width2,siz);
0382     leg = {'EXP','LUKE'};
0383 else
0384     leg = NaN;
0385 end
0386 graph1D_jd(var.list,1e3*dI_rf/opt.Ifac,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0387 set(gca,'xtick',xtick);
0388 if isfield(opt,'ytick_I') && ~any(isnan(opt.ytick_I)),
0389     set(gca,'ytick',opt.ytick_I);
0390 end
0391 %
0392 figure(102),clf
0393 %
0394 if isfield(opt,'ylim_rho'),
0395     ylim = opt.ylim_rho;
0396 else
0397     ylim = [0,1.2];
0398 end
0399 %
0400 ylab = '\rho_{dep}';
0401 tit = '';
0402 %
0403 graph1D_jd(var.list,drho,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker2,color3,width2,siz,gca,red,lspace2,bspace2,ddrho/2);
0404 set(gca,'xtick',xtick);
0405 if isfield(opt,'ytick_rho') && ~any(isnan(opt.ytick_rho)),
0406     set(gca,'ytick',opt.ytick_rho);
0407 end
0408 %
0409 figure(2),clf
0410 %
0411 if isfield(opt,'ylim_eta'),
0412     ylim = opt.ylim_eta;
0413 else
0414     ylim = [0,max(deta)]*1.2;
0415 end
0416 %
0417 ylab = '\eta_{LH} (A/W/m^2) = n_{19}IR/P';
0418 tit = '';
0419 %
0420 graph1D_jd(var.list,deta,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0421 set(gca,'xtick',xtick);
0422 if isfield(opt,'ytick_eta') && ~any(isnan(opt.ytick_eta)),
0423     set(gca,'ytick',opt.ytick_eta);
0424 end
0425 %
0426 figure(3),clf
0427 %
0428 if isfield(opt,'ylim_xi'),
0429     ylim = opt.ylim_xi;
0430 else
0431     ylim = [0,max(dksirf)]*1.2;
0432 end
0433 %
0434 ylab = '\xi_{LH} = 32.7\times n_{20}IR/(PT)';
0435 tit = '';
0436 %
0437 graph1D_jd(var.list,dksirf,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0438 set(gca,'xtick',xtick);
0439 if isfield(opt,'ytick_xi') && ~any(isnan(opt.ytick_xi)),
0440     set(gca,'ytick',opt.ytick_xi);
0441 end
0442 %
0443 figure(4),clf
0444 %
0445 if isfield(opt,'ylim_P'),
0446     ylim = opt.ylim_P;
0447 else
0448     ylim = [0,max(dP0_rf)]*1.2;
0449 end
0450 %
0451 ylab = 'P (MW)';
0452 leg = {'Coupled power','Absorbed power (C3PO)','Absorbed power (LUKE)','Location','NorthWest'};
0453 %
0454 graph1D_jd(var.list,dP0_rf,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width2,siz);
0455 graph1D_jd(var.list,dP_lin_rf,0,0,'','','',NaN,xlim,ylim,style,marker2,color3,width2,siz);
0456 graph1D_jd(var.list,dP_rf,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0457 set(gca,'xtick',xtick);
0458 if isfield(opt,'ytick_P') && ~any(isnan(opt.ytick_P)),
0459     set(gca,'ytick',opt.ytick_P);
0460 end
0461 %
0462 % --------------------- Weighed ray length -----------------------
0463 %
0464 figure(5),clf
0465 %
0466 if isfield(opt,'ylim_s'),
0467     ylim = opt.ylim_s;
0468 else
0469     ylim = [0,max([ds_c3po,ds_luke])]*1.2;
0470 end
0471 %
0472 ylab = 'Weighed ray length (m)';
0473 leg = {'C3PO','LUKE'};
0474 %
0475 graph1D_jd(var.list,ds_c3po,0,0,'','','',NaN,xlim,ylim,style,marker2,color3,width2,siz);
0476 graph1D_jd(var.list,ds_luke,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0477 set(gca,'xtick',xtick);
0478 if isfield(opt,'ytick_s') && ~any(isnan(opt.ytick_s)),
0479     set(gca,'ytick',opt.ytick_s);
0480 end
0481 %
0482 % --------------------- N // -----------------------
0483 %
0484 figure(6),clf
0485 %
0486 if isfield(opt,'ylim_Npar'),
0487     ylim = opt.ylim_Npar;
0488 else
0489     ylim = [0,max([dNpar0,dNpar_c3po,dNpar_luke])]*1.2;
0490 end
0491 %
0492 ylab = '';
0493 leg = {'Initial <N//>','C3PO absorbed <N//>','LUKE absorbed <N//>'};
0494 %
0495 graph1D_jd(var.list,dNpar0,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width2,siz);
0496 graph1D_jd(var.list,dNpar_c3po,0,0,'','','',NaN,xlim,ylim,style,marker2,color3,width2,siz);
0497 graph1D_jd(var.list,dNpar_luke,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker1,color2,width2,siz,gca,red,lspace2,bspace2);
0498 set(gca,'xtick',xtick);
0499 if isfield(opt,'ytick_Npar') && ~any(isnan(opt.ytick_Npar)),
0500     set(gca,'ytick',opt.ytick_Npar);
0501 end
0502 %
0503 mask_disp = ~isnan(dI_rf(ivar_disp));
0504 ivar_disp = ivar_disp(mask_disp);
0505 var_disp = var.liststr(ivar_disp);
0506 ndisp = length(ivar_disp);
0507 %
0508 if ndisp >= 1,
0509     %
0510     figure(7),clf
0511     %
0512     if isfield(opt,'ylim_xJ'),
0513         ylim = opt.ylim_xJ;
0514     else
0515         ylim = NaN;
0516     end
0517     %
0518     xlab = {'r/a'};
0519     xlim = [0,1];
0520     ylab = 'J (MA/m^2)';
0521     tit = '';
0522     leg = {};
0523     %
0524     for idisp = 1:ndisp-1,
0525         graph1D_jd(dxrho{ivar_disp(idisp)},dxJ{ivar_disp(idisp)},0,0,'','','',NaN,xlim,ylim,style,markers{idisp},colors{idisp},width2,siz);
0526         leg = [leg,[var.lab,' = ',var_disp{idisp}]];
0527     end
0528     %
0529     leg = [leg,[var.lab,' = ',var_disp{ndisp}]];
0530     graph1D_jd(dxrho{ivar_disp(ndisp)},dxJ{ivar_disp(ndisp)},0,0,xlab,ylab,tit,leg,xlim,ylim,style,markers{ndisp},colors{ndisp},width2,siz,gca,red,lspace2,bspace2);
0531     %
0532     if isfield(opt,'xtick_xJ') && ~any(isnan(opt.xtick_xJ)),
0533         set(gca,'xtick',opt.xtick_xJ);
0534     end
0535     if isfield(opt,'ytick_xJ') && ~any(isnan(opt.ytick_xJ)),
0536         set(gca,'ytick',opt.ytick_xJ);
0537     end
0538     %
0539     figure(8),clf
0540     %
0541     if isfield(opt,'ylim_xP'),
0542         ylim = opt.ylim_xP;
0543     else
0544         ylim = NaN;
0545     end
0546     %
0547     ylab = 'P (MW/m^3)';
0548     %
0549     for idisp = 1:ndisp-1,
0550         graph1D_jd(dxrho{ivar_disp(idisp)},dxP{ivar_disp(idisp)},0,0,'','','',NaN,xlim,ylim,style,markers{idisp},colors{idisp},width2,siz);
0551     end
0552     graph1D_jd(dxrho{ivar_disp(ndisp)},dxP{ivar_disp(ndisp)},0,0,xlab,ylab,tit,leg,xlim,ylim,style,markers{ndisp},colors{ndisp},width2,siz,gca,red,lspace2,bspace2);
0553     if isfield(opt,'xtick_xP') && ~any(isnan(opt.xtick_xP)),
0554         set(gca,'xtick',opt.xtick_xP);
0555     end
0556     if isfield(opt,'ytick_xP') && ~any(isnan(opt.ytick_xP)),
0557         set(gca,'ytick',opt.ytick_xP);
0558     end
0559     %
0560     figure(9),clf
0561     %
0562     if isfield(opt,'ylim_xP'),
0563         ylim = opt.ylim_xP;
0564     else
0565         ylim = NaN;
0566     end
0567     %
0568     ylab = 'dp_{RF}/dt (N/m^3)';
0569     %
0570     for idisp = 1:ndisp-1,
0571         graph1D_jd(dxrho{ivar_disp(idisp)},dxdpdt_rf{ivar_disp(idisp)},0,0,'','','',NaN,xlim,ylim,style,markers{idisp},colors{idisp},width2,siz);
0572     end
0573     graph1D_jd(dxrho{ivar_disp(ndisp)},dxdpdt_rf{ivar_disp(ndisp)},0,0,xlab,ylab,tit,leg,xlim,ylim,style,markers{ndisp},colors{ndisp},width2,siz,gca,red,lspace2,bspace2);
0574     if isfield(opt,'xtick_xdpdt_rf') && ~any(isnan(opt.xtick_xdpdt_rf)),
0575         set(gca,'xtick',opt.xtick_xdpdt_rf);
0576     end
0577     if isfield(opt,'ytick_xdpdt_rf') && ~any(isnan(opt.ytick_xdpdt_rf)),
0578         set(gca,'ytick',opt.ytick_xdpdt_rf);
0579     end
0580     %
0581 end
0582 %
0583 p_opt = input_dke_yp('Printing option (-1->2)',-1,-1:2);%option for printing
0584 %
0585 print_jd(p_opt,[figstr,'_I_rho'],[datapath,filesep,'figures',filesep],1);
0586 print_jd(p_opt,[figstr,'_I'],[figdir,filesep,'figures',filesep],101);
0587 print_jd(p_opt,[figstr,'_rho'],[figdir,filesep,'figures',filesep],102);
0588 print_jd(p_opt,[figstr,'_eta'],[figdir,filesep,'figures',filesep],2);
0589 print_jd(p_opt,[figstr,'_ksi'],[figdir,filesep,'figures',filesep],3);
0590 print_jd(p_opt,[figstr,'_P'],[figdir,filesep,'figures',filesep],4);
0591 print_jd(p_opt,[figstr,'_s'],[figdir,filesep,'figures',filesep],5);
0592 print_jd(p_opt,[figstr,'_Npar'],[figdir,filesep,'figures',filesep],6);
0593 %
0594 if ndisp >= 1,
0595     %
0596     print_jd(p_opt,[figstr,'_xJ'],[figdir,filesep,'figures',filesep],7);
0597     print_jd(p_opt,[figstr,'_xP'],[figdir,filesep,'figures',filesep],8);
0598     %
0599 end
0600 %