0001 function [data_proc,ir_display] = proc_luke_jd(data,ir_display,p_opt,dp_cyl,pnmax,subtitle,opt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 7,
0013
0014
0015
0016 opt.peak = '';
0017 opt.wave = '';
0018 opt.jpeak = 0;
0019 opt.peakmode = 0;
0020 opt.blist = '';
0021
0022
0023
0024 opt.rho = '';
0025 opt.spec = 1;
0026 opt.its_rf = true;
0027 opt.diaryname = [];
0028 end
0029 if nargin < 6,
0030 subtitle = '';
0031 end
0032 if nargin < 5,
0033 pnmax = NaN;
0034 end
0035 if nargin < 4,
0036 dp_cyl = 0.1;
0037 end
0038 if nargin < 3,
0039 p_opt = -1;
0040 end
0041 if nargin < 2,
0042 ir_display = struct;
0043 end
0044
0045 if isstruct(ir_display),
0046
0047 opt = ir_display;
0048
0049 if isfield(opt,'ir_display'),
0050 ir_display = opt.ir_display;
0051 else
0052 ir_display = -1;
0053 end
0054
0055 if isfield(opt,'p_opt'),
0056 p_opt = opt.p_opt;
0057 end
0058 if isfield(opt,'pnmax'),
0059 pnmax = opt.pnmax;
0060 end
0061 if ~isfield(opt,'peak'),
0062 opt.peak = 'y';
0063 end
0064 if ~isfield(opt,'wave'),
0065 opt.wave = 'y';
0066 end
0067 if ~isfield(opt,'peakmode'),
0068 opt.peakmode = 1;
0069 end
0070 if ~isfield(opt,'rho'),
0071 opt.rho = 't';
0072 end
0073 if ~isfield(opt,'spec'),
0074 opt.spec = 0;
0075 end
0076 if ~isfield(opt,'its_rf'),
0077 opt.its_rf = [];
0078 end
0079 if ~isfield(opt,'blist'),
0080 opt.blist = [];
0081 end
0082 if ~isfield(opt,'diaryname'),
0083 opt.diaryname = [];
0084 end
0085 end
0086
0087 if ~isfield(opt,'peakmode'),
0088 opt.peakmode = 0;
0089 end
0090
0091 if ~isfield(opt,'jpeak'),
0092 opt.jpeak = 0;
0093 end
0094
0095 if ~isfield(opt,'rho') || isempty(opt.rho),
0096 opt.rho = input_dke_yp('Option for rho definition : ''g'' (geometric), ''p'' (poloidal flux), ''t'' (toroidal flux), ''v'' (volumic rho)','g',{'g','p','t','v'},'',[1,1]);
0097 end
0098
0099 if ~isfield(opt,'spec') || isempty(opt.spec)
0100 opt.spec = 1;
0101 end
0102
0103 if isfield(data,'lukestructs'),
0104 data = data.lukestructs;
0105 end
0106 if isfield(data,'output'),
0107 data = data.output;
0108 end
0109
0110 if ~isfield(opt,'gui'),
0111 opt.gui = false;
0112 end
0113
0114 if ~isfield(opt,'style'),
0115 opt.style = struct;
0116 end
0117 if ~isfield(opt.style,'distselectstyle'),
0118 opt.style.distselectstyle = struct;
0119 end
0120 if ~isfield(opt.style,'editselectstyle'),
0121 opt.style.editselectstyle = struct;
0122 end
0123
0124 style0 = 'none';
0125 style = '-';
0126 style2 = '--';
0127 marker = 'none';
0128 marker1 = '+';
0129 marker2 = 'o';
0130 marker3 = 's';
0131 marker4 = '.';
0132 color1 = 'k';
0133 color2 = 'r';
0134 color3 = 'b';
0135 color4 = 'g';
0136 color5 = 'm';
0137 color6 = 'c';
0138 color7 = 'y';
0139 colors = {color2,color3,color4,color5,color6,color7};
0140
0141 width = 0.5;
0142 width2 = 2;
0143
0144 if ~isfield(opt,'font'),
0145 opt.font = 20+14i;
0146 end
0147
0148 if isfield(data,'data_proc'),
0149
0150 data_proc = data.data_proc;
0151
0152 opt = conc_struct_jd(opt,selectfields_jd(data_proc.opt,{'wave','peak','jpeak','peakmode'}));
0153
0154 nr = length(data_proc.radial.xrhoG);
0155 nn_rf = length(data_proc.scalar.n_rf_list);
0156 if opt.wave == 'y',
0157 nw = length(data_proc.wave.wp0_rf_2piRp);
0158 end
0159 nit_rf = length(data_proc.simul.residue_rf);
0160
0161 if ~isfield(data_proc.scalar,'P_rf_rhoG'),
0162
0163 data_proc.scalar.P_rf_rhoG = NaN;
0164 data_proc.scalar.P_rf_rhoP = NaN;
0165 data_proc.scalar.P_rf_rhoT = NaN;
0166 data_proc.scalar.P_rf_rhoV = NaN;
0167
0168 data_proc.scalar.P_rf_drhoG = NaN;
0169 data_proc.scalar.P_rf_drhoP = NaN;
0170 data_proc.scalar.P_rf_drhoT = NaN;
0171 data_proc.scalar.P_rf_drhoV = NaN;
0172
0173 xpsin_S = zeros(1,nr+1);
0174 for ir = 1:nr,
0175 xpsin_S(ir+1) = 2*data_proc.radial.xrhoP(ir)^2 - xpsin_S(ir);
0176 end
0177 xrhoP_S = sqrt(xpsin_S);
0178 xrhoG_S = [0,interp1(data_proc.radial.xrhoP,data_proc.radial.xrhoG,xrhoP_S(2:end-1),'spline'),1];
0179 xrhoT_S = [0,interp1(data_proc.radial.xrhoP,data_proc.radial.xrhoT,xrhoP_S(2:end-1),'spline'),1];
0180 xrhoV_S = [0,interp1(data_proc.radial.xrhoP,data_proc.radial.xrhoV,xrhoP_S(2:end-1),'spline'),1];
0181
0182 data_proc.radial.xdrhoG = diff(xrhoG_S);
0183 data_proc.radial.xdrhoP = diff(xrhoP_S);
0184 data_proc.radial.xdrhoT = diff(xrhoT_S);
0185 data_proc.radial.xdrhoV = diff(xrhoV_S);
0186
0187 if opt.wave == 'y',
0188 data_proc.wave.wids = cell(1,nw);
0189
0190 data_proc.wave.wP_rf_rhoG = data_proc.wave.wrhoG;
0191 data_proc.wave.wP_rf_rhoP = data_proc.wave.wrhoP;
0192 data_proc.wave.wP_rf_rhoT = data_proc.wave.wrhoT;
0193 data_proc.wave.wP_rf_rhoV = data_proc.wave.wrhoV;
0194 end
0195
0196 end
0197
0198 else
0199
0200 Znorm = data.Znorm;
0201 Zcurr = data.Zcurr;
0202 ZP0 = data.ZP0;
0203 if isfield(data.dke_out,'Zmom'),
0204 Zmom = data.dke_out.Zmom;
0205 else
0206 Zmom = struct;
0207 end
0208 dke_out = data.dke_out;
0209 radialDKE = data.radialDKE;
0210 equilDKE = data.equilDKE;
0211 momentumDKE = data.momentumDKE;
0212 gridDKE = data.gridDKE;
0213 Zmomcoef = data.Zmomcoef;
0214 Zbouncecoef = data.Zbouncecoef;
0215 Zmripple = data.Zmripple;
0216 mksa = data.mksa;
0217
0218 if isfield(dke_out,'dkeparam')
0219 dkeparam = dke_out.dkeparam;
0220 elseif isfield(data,'dkeparam')
0221 dkeparam = data.dkeparam;
0222 else
0223 dkeparam = struct;
0224 end
0225
0226 if ~isfield(dke_out,'XXf0'),
0227 dke_out.XXf0 = legendre2f_yp(dkeparam,momentumDKE,dke_out.XXf0_interp);
0228 end
0229
0230 if ~isfield(dke_out,'XXfM'),
0231 dke_out.XXfM = legendre2f_yp(dkeparam,momentumDKE,dke_out.XXfM_interp);
0232 end
0233
0234 if ~isfield(dke_out,'XXfinit'),
0235 dke_out.XXfinit = legendre2f_yp(dkeparam,momentumDKE,dke_out.XXfinit_interp);
0236 end
0237
0238 if isfield(dke_out,'waves'),
0239 waves = proc_wave_alphaphi_luke_jd(data,dke_out.waves,1);
0240 else
0241 waves = {};
0242 end
0243
0244 xrhoG = equilDKE.xrho;
0245 xrhoP = equilDKE.xrhoP;
0246 xrhoT = equilDKE.xrhoT;
0247 xrhoV = equilDKE.xrhoV;
0248
0249
0250
0251 xrhoG_S = radialDKE.xrho_S_dke;
0252 xrhoP_S = sqrt(radialDKE.xpsin_S_dke);
0253 xrhoT_S = [0,interp1(xrhoG,xrhoT,radialDKE.xrho_S_dke(2:end-1),'spline'),1];
0254 xrhoV_S = [0,interp1(xrhoG,xrhoV,radialDKE.xrho_S_dke(2:end-1),'spline'),1];
0255
0256 xdrhoG = diff(xrhoG_S);
0257 xdrhoP = diff(xrhoP_S);
0258 xdrhoT = diff(xrhoT_S);
0259 xdrhoV = diff(xrhoV_S);
0260
0261 xTe = equilDKE.xTe;
0262 xne = equilDKE.xne;
0263
0264 nr = length(xrhoG);
0265 nw = length(waves);
0266
0267 yb = dke_out.yb;
0268
0269 if isfield(dke_out,'waveparam') && strcmp(dke_out.waveparam.model,'RT'),
0270 nfrags = length(dke_out.waveparam.biy_rf);
0271 else
0272 nfrags = NaN;
0273 end
0274
0275 if isfield(dkeparam,'timevol'),
0276 timevol = dkeparam.timevol;
0277 else
0278 timevol = 0;
0279 end
0280
0281 scocurr = sign(equilDKE.psia_apRp)*sign(equilDKE.XBphi(1,1));
0282
0283
0284
0285 if ~isempty(waves),
0286
0287 wxP_rf_lin = zeros(nw,nr);
0288 wxP_rf_lin_coll = zeros(nw,nr);
0289 wxP_rf_coll_luke = zeros(nw,nr);
0290
0291 wP_rf_rhoG_lin = zeros(1,nw);
0292 wP_rf_rhoP_lin = zeros(1,nw);
0293 wP_rf_rhoT_lin = zeros(1,nw);
0294 wP_rf_rhoV_lin = zeros(1,nw);
0295
0296 wP_rf_rhoG_lin_coll = zeros(1,nw);
0297 wP_rf_rhoP_lin_coll = zeros(1,nw);
0298 wP_rf_rhoT_lin_coll = zeros(1,nw);
0299 wP_rf_rhoV_lin_coll = zeros(1,nw);
0300
0301 wP_rf_rhoG_coll_luke = zeros(1,nw);
0302 wP_rf_rhoP_coll_luke = zeros(1,nw);
0303 wP_rf_rhoT_coll_luke = zeros(1,nw);
0304 wP_rf_rhoV_coll_luke = zeros(1,nw);
0305
0306 for iw = 1:nw,
0307
0308 if ~isfield(waves{iw},'model') || strcmp(waves{iw}.model,'RT'),
0309
0310 nb(iw) = length(waves{iw}.rays);
0311
0312 for iray = 1:nb(iw),
0313
0314 ray = waves{iw}.rays{iray};
0315
0316 if isfield(ray,'sP_2piRp_lin') && isfield(ray,'srho'),
0317 wxP_rf_lin(iw,:) = wxP_rf_lin(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin,xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0318 else
0319 wxP_rf_lin(iw,:) = NaN;
0320 end
0321
0322 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
0323
0324 if isfield(ray,'sP_2piRp_lin_coll') && isfield(ray,'srho'),
0325 wxP_rf_lin_coll(iw,:) = wxP_rf_lin_coll(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin_coll,xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0326 else
0327 wxP_rf_lin_coll(iw,:) = NaN;
0328 end
0329
0330 if isfield(ray,'sP_2piRp_coll_luke') && isfield(ray,'srho'),
0331 wxP_rf_coll_luke(iw,:) = wxP_rf_coll_luke(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_coll_luke{1},xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0332 else
0333 wxP_rf_coll_luke(iw,:) = NaN;
0334 end
0335
0336 end
0337
0338 end
0339
0340 elseif strcmp(waves{iw}.model,'FW'),
0341
0342 nr_fw = length(waves{iw}.rho);
0343 xwmask = zeros(nr,nr_fw);
0344
0345 for ir_fw = 1:nr_fw,
0346 if waves{iw}.rhotype == 'g',
0347 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoG_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoG_S(1:end-1)]))]).'./xdrhoG.';
0348 elseif waves{iw}.rhotype == 'p',
0349 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoP_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoP_S(1:end-1)]))]).'./xdrhoP.';
0350 elseif waves{iw}.rhotype == 't',
0351 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoT_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoT_S(1:end-1)]))]).'./xdrhoT.';
0352 elseif waves{iw}.rhotype == 'v',
0353 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoV_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoV_S(1:end-1)]))]).'./xdrhoV.';
0354 else
0355 error('rho type not recognized')
0356 end
0357 end
0358
0359 wpe_eve = 1e-6*(waves{iw}.pe_eve(1:end - 1) + waves{iw}.pe_eve(2:end))/2;
0360
0361 wxP_rf_lin(iw,:) = (xwmask*wpe_eve).';
0362
0363 else
0364 error('Wave type not recognized.')
0365 end
0366
0367 [wP_rf_rhoG_lin(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0368 [wP_rf_rhoP_lin(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0369 [wP_rf_rhoT_lin(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0370 [wP_rf_rhoV_lin(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0371
0372 [wP_rf_rhoG_lin_coll(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0373 [wP_rf_rhoP_lin_coll(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0374 [wP_rf_rhoT_lin_coll(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0375 [wP_rf_rhoV_lin_coll(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0376
0377 [wP_rf_rhoG_coll_luke(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0378 [wP_rf_rhoP_coll_luke(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0379 [wP_rf_rhoT_coll_luke(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0380 [wP_rf_rhoV_coll_luke(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0381
0382 end
0383
0384 wxp_rf_2piRp_lin = wxP_rf_lin.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0385 wp_rf_2piRp_lin = sum(wxp_rf_2piRp_lin,2);
0386 xp_rf_2piRp_lin = sum(wxp_rf_2piRp_lin,1);
0387 p_rf_2piRp_lin = sum(wp_rf_2piRp_lin);
0388
0389 wxp_rf_2piRp_lin_coll = wxP_rf_lin_coll.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0390 wp_rf_2piRp_lin_coll = sum(wxp_rf_2piRp_lin_coll,2);
0391 xp_rf_2piRp_lin_coll = sum(wxp_rf_2piRp_lin_coll,1);
0392 p_rf_2piRp_lin_coll = sum(wp_rf_2piRp_lin_coll);
0393
0394 wxp_rf_2piRp_coll_luke = wxP_rf_coll_luke.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0395 wp_rf_2piRp_coll_luke = sum(wxp_rf_2piRp_coll_luke,2);
0396 xp_rf_2piRp_coll_luke = sum(wxp_rf_2piRp_coll_luke,1);
0397 p_rf_2piRp_coll_luke = sum(wp_rf_2piRp_coll_luke);
0398
0399 if ~all(isnan(wP_rf_rhoG_lin)),
0400 wTe_lin = interp1(xrhoG,xTe,wP_rf_rhoG_lin);
0401 wne_lin = interp1(xrhoG,xne,wP_rf_rhoG_lin);
0402 else
0403 wTe_lin = NaN(1,nw);
0404 wne_lin = NaN(1,nw);
0405 end
0406
0407 if ~all(isnan(wP_rf_rhoG_lin_coll)),
0408 wTe_lin_coll = interp1(xrhoG,xTe,wP_rf_rhoG_lin_coll);
0409 wne_lin_coll = interp1(xrhoG,xne,wP_rf_rhoG_lin_coll);
0410 else
0411 wTe_lin_coll = NaN(1,nw);
0412 wne_lin_coll = NaN(1,nw);
0413 end
0414
0415 if ~all(isnan(wP_rf_rhoG_coll_luke)),
0416 wTe_coll_luke = interp1(xrhoG,xTe,wP_rf_rhoG_coll_luke);
0417 wne_coll_luke = interp1(xrhoG,xne,wP_rf_rhoG_coll_luke);
0418 else
0419 wTe_coll_luke = NaN(1,nw);
0420 wne_coll_luke = NaN(1,nw);
0421 end
0422
0423 [P_rf_rhoG_lin,P_rf_drhoG_lin] = calc_peak_jd(xrhoG,xdrhoG,xp_rf_2piRp_lin./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0424 [P_rf_rhoP_lin,P_rf_drhoP_lin] = calc_peak_jd(xrhoP,xdrhoP,xp_rf_2piRp_lin./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0425 [P_rf_rhoT_lin,P_rf_drhoT_lin] = calc_peak_jd(xrhoT,xdrhoT,xp_rf_2piRp_lin./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0426 [P_rf_rhoV_lin,P_rf_drhoV_lin] = calc_peak_jd(xrhoV,xdrhoV,xp_rf_2piRp_lin./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0427
0428 [P_rf_rhoG_lin_coll,P_rf_drhoG_lin_coll] = calc_peak_jd(xrhoG,xdrhoG,xp_rf_2piRp_lin_coll./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0429 [P_rf_rhoP_lin_coll,P_rf_drhoP_lin_coll] = calc_peak_jd(xrhoP,xdrhoP,xp_rf_2piRp_lin_coll./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0430 [P_rf_rhoT_lin_coll,P_rf_drhoT_lin_coll] = calc_peak_jd(xrhoT,xdrhoT,xp_rf_2piRp_lin_coll./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0431 [P_rf_rhoV_lin_coll,P_rf_drhoV_lin_coll] = calc_peak_jd(xrhoV,xdrhoV,xp_rf_2piRp_lin_coll./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0432
0433 [P_rf_rhoG_coll_luke,P_rf_drhoG_coll_luke] = calc_peak_jd(xrhoG,xdrhoG,xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0434 [P_rf_rhoP_coll_luke,P_rf_drhoP_coll_luke] = calc_peak_jd(xrhoP,xdrhoP,xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0435 [P_rf_rhoT_coll_luke,P_rf_drhoT_coll_luke] = calc_peak_jd(xrhoT,xdrhoT,xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0436 [P_rf_rhoV_coll_luke,P_rf_drhoV_coll_luke] = calc_peak_jd(xrhoV,xdrhoV,xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0437
0438 if ~isnan(P_rf_rhoG_lin)
0439 P_rf_Te_lin = interp1(xrhoG,xTe,P_rf_rhoG_lin);
0440 P_rf_ne_lin = interp1(xrhoG,xne,P_rf_rhoG_lin);
0441 else
0442 P_rf_Te_lin = NaN;
0443 P_rf_ne_lin = NaN;
0444 end
0445
0446 if ~isnan(P_rf_rhoG_lin_coll)
0447 P_rf_Te_lin_coll = interp1(xrhoG,xTe,P_rf_rhoG_lin_coll);
0448 P_rf_ne_lin_coll = interp1(xrhoG,xne,P_rf_rhoG_lin_coll);
0449 else
0450 P_rf_Te_lin_coll = NaN;
0451 P_rf_ne_lin_coll = NaN;
0452 end
0453
0454 if ~isnan(P_rf_rhoG_coll_luke)
0455 P_rf_Te_coll_luke = interp1(xrhoG,xTe,P_rf_rhoG_coll_luke);
0456 P_rf_ne_coll_luke = interp1(xrhoG,xne,P_rf_rhoG_coll_luke);
0457 else
0458 P_rf_Te_coll_luke = NaN;
0459 P_rf_ne_coll_luke = NaN;
0460 end
0461
0462 else
0463
0464 xp_rf_2piRp_lin = NaN;
0465 p_rf_2piRp_lin = NaN;
0466
0467 xp_rf_2piRp_lin_coll = NaN;
0468 p_rf_2piRp_lin_coll = NaN;
0469
0470 xp_rf_2piRp_coll_luke = NaN;
0471 p_rf_2piRp_coll_luke = NaN;
0472
0473 P_rf_rhoG_lin = NaN;P_rf_drhoG_lin = NaN;
0474 P_rf_rhoP_lin = NaN;P_rf_drhoP_lin = NaN;
0475 P_rf_rhoT_lin = NaN;P_rf_drhoT_lin = NaN;
0476 P_rf_rhoV_lin = NaN;P_rf_drhoV_lin = NaN;
0477
0478 P_rf_rhoG_lin_coll = NaN;P_rf_drhoG_lin_coll = NaN;
0479 P_rf_rhoP_lin_coll = NaN;P_rf_drhoP_lin_coll = NaN;
0480 P_rf_rhoT_lin_coll = NaN;P_rf_drhoT_lin_coll = NaN;
0481 P_rf_rhoV_lin_coll = NaN;P_rf_drhoV_lin_coll = NaN;
0482
0483 P_rf_rhoG_coll_luke = NaN;P_rf_drhoG_coll_luke = NaN;
0484 P_rf_rhoP_coll_luke = NaN;P_rf_drhoP_coll_luke = NaN;
0485 P_rf_rhoT_coll_luke = NaN;P_rf_drhoT_coll_luke = NaN;
0486 P_rf_rhoV_coll_luke = NaN;P_rf_drhoV_coll_luke = NaN;
0487
0488 P_rf_Te_lin = NaN;P_rf_ne_lin = NaN;
0489 P_rf_Te_lin_coll = NaN;P_rf_ne_lin_coll = NaN;
0490 P_rf_Te_coll_luke = NaN;P_rf_ne_coll_luke = NaN;
0491
0492 end
0493
0494
0495
0496 J_tot = Zcurr(end).x_0_fsav*mksa.j_ref*scocurr;
0497
0498
0499
0500 P_rf = ZP0(end).x_rf_fsav.'*mksa.P_ref;
0501 yP_rf = ZP0(end).xy_rf_fsav.'*mksa.P_ref;
0502 ynP_rf = ZP0(end).xyn_rf_fsav*mksa.P_ref;
0503 P_e = ZP0(end).x_e_fsav*mksa.P_ref;
0504
0505 if ~isempty(waves)
0506 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
0507 P_tot = xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke + P_e;
0508 else
0509 P_tot = P_rf + P_e;
0510 end
0511 else
0512 P_tot = P_e;
0513 end
0514
0515
0516
0517 if isfield(Zmom,'x_rf'),
0518 xdpdt_rf = Zmom(end).x_rf*mksa.dpdt_ref;
0519 else
0520 xdpdt_rf = NaN(1,nr);
0521 end
0522
0523 if isfield(dke_out,'residue_rf'),
0524 if iscell(dke_out.residue_rf),
0525 residue_rf = dke_out.residue_rf{end};
0526 else
0527 residue_rf = dke_out.residue_rf;
0528 end
0529 else
0530 residue_rf = NaN;
0531 end
0532
0533 nit_rf = length(residue_rf);
0534
0535 if isfield(dke_out,'zS'),
0536 zS = dke_out.zS;
0537
0538 if size(dke_out.zP_0_2piRp_mod,1) >= nit_rf,
0539 zP_0_2piRp_mod = dke_out.zP_0_2piRp_mod(1:nit_rf,end);
0540 if isfield(dke_out,'zP_0_2piRp_mod_coll')
0541 zP_0_2piRp_mod_coll = dke_out.zP_0_2piRp_mod_coll(1:nit_rf,end);
0542 end
0543 else
0544 zP_0_2piRp_mod = dke_out.zP_0_2piRp_mod(1,end);
0545 end
0546 else
0547 zS = '';
0548 zP_0_2piRp_mod = '';
0549 end
0550
0551 ny = size(yP_rf,1);
0552 nn_rf = size(ynP_rf,3);
0553
0554 if isfield(ZP0(end),'n_rf_list'),
0555 n_rf_list = ZP0(end).n_rf_list;
0556 elseif isfield(dke_out,'waveparam'),
0557 n_rf_list = dke_out.waveparam.n_rf_list;
0558 elseif isfield(dkeparam,'n_rf_list'),
0559 n_rf_list = dkeparam.n_rf_list;
0560 else
0561 n_rf_list = NaN(1,nn_rf);
0562 disp('WARNING : n_rf_list not determined.')
0563 end
0564
0565 xI = J_tot.*equilDKE.xdA_dke;
0566 I_tot = sum(xI);
0567
0568 [jrhoG,jdrhoG] = calc_peak_jd(xrhoG,xdrhoG,J_tot,equilDKE.xdA_dke,opt.jpeak,opt.peakmode);
0569 [jrhoP,jdrhoP] = calc_peak_jd(xrhoP,xdrhoP,J_tot,equilDKE.xdA_dke,opt.jpeak,opt.peakmode);
0570 [jrhoT,jdrhoT] = calc_peak_jd(xrhoT,xdrhoT,J_tot,equilDKE.xdA_dke,opt.jpeak,opt.peakmode);
0571 [jrhoV,jdrhoV] = calc_peak_jd(xrhoV,xdrhoV,J_tot,equilDKE.xdA_dke,opt.jpeak,opt.peakmode);
0572
0573 if sum(P_rf) ~= 0
0574 [P_rf_rhoG,P_rf_drhoG] = calc_peak_jd(xrhoG,xdrhoG,P_rf,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0575 [P_rf_rhoP,P_rf_drhoP] = calc_peak_jd(xrhoP,xdrhoP,P_rf,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0576 [P_rf_rhoT,P_rf_drhoT] = calc_peak_jd(xrhoT,xdrhoT,P_rf,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0577 [P_rf_rhoV,P_rf_drhoV] = calc_peak_jd(xrhoV,xdrhoV,P_rf,equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0578
0579 else
0580 P_rf_rhoG = NaN;P_rf_drhoG = NaN;
0581 P_rf_rhoP = NaN;P_rf_drhoP = NaN;
0582 P_rf_rhoT = NaN;P_rf_drhoT = NaN;
0583 P_rf_rhoV = NaN;P_rf_drhoV = NaN;
0584 end
0585
0586 if length(xrhoG) > 1,
0587 jTe = interp1(xrhoG,xTe,jrhoG);
0588 jne = interp1(xrhoG,xne,jrhoG);
0589
0590 if ~isnan(P_rf_rhoG)
0591 P_rf_Te = interp1(xrhoG,xTe,P_rf_rhoG);
0592 P_rf_ne = interp1(xrhoG,xne,P_rf_rhoG);
0593 else
0594 P_rf_Te = NaN;
0595 P_rf_ne = NaN;
0596 end
0597 else
0598 jTe = jrhoG;
0599 jne = jrhoG;
0600
0601 if ~isnan(P_rf_rhoG)
0602 P_rf_Te = P_rf_rhoG;
0603 P_rf_ne = P_rf_rhoG;
0604 else
0605 P_rf_Te = NaN;
0606 P_rf_ne = NaN;
0607 end
0608 end
0609
0610 if isempty(opt.peak) && sum(P_rf) ~= 0,
0611 opt.peak = input_dke_yp('Do you want to display specific information for each absorption peak? (y/n)','n',{'n','y'},'',[1,1]);
0612 end
0613
0614 if isempty(opt.wave) && sum(P_rf) ~= 0,
0615 opt.wave = input_dke_yp('Do you want to display specific information for each wave? (y/n)','n',{'n','y'},'',[1,1]);
0616 end
0617
0618 if nw == 1 && (~isfield(waves{1},'model') || strcmp(waves{1}.model,'RT')),
0619 if ~isfield(opt,'blist') || (isscalar(opt.blist) && isnan(opt.blist)),
0620 nb = length(waves{1}.rays);
0621 opt_blist = input(['You can group the ',num2str(nb),' rays by vectors in cell (leave empty otherwise) :']);
0622 else
0623 opt_blist = opt.blist;
0624 end
0625 else
0626 opt_blist = [];
0627 end
0628
0629 xp_rf_2piRp = P_rf.*equilDKE.xdV_2piRp_dke;
0630 yxp_rf_2piRp = yP_rf.*repmat(equilDKE.xdV_2piRp_dke,[ny,1]);
0631 yxnp_rf_2piRp = ynP_rf.*repmat(equilDKE.xdV_2piRp_dke.',[1,ny,nn_rf]);
0632
0633 p_rf_2piRp = sum(xp_rf_2piRp);
0634 np_rf_2piRp = sum(reshape(sum(yxnp_rf_2piRp,2),[nr,nn_rf]),1);
0635 xnp_rf_2piRp = reshape(sum(yxnp_rf_2piRp,2),[nr,nn_rf]);
0636
0637 p_e_2piRp = sum(P_e.*equilDKE.xdV_2piRp_dke);
0638 p_tot_2piRp = sum(P_tot.*equilDKE.xdV_2piRp_dke);
0639 EtaRp = I_tot/(p_tot_2piRp*2*pi);
0640 EtaRFRp = jne*1e-19*I_tot/(p_rf_2piRp*2*pi);
0641 ne_bar0 = sum(gradient(equilDKE.xrho).*equilDKE.xne)*1e-19;
0642 EtaRFRp_nebar0 = ne_bar0*I_tot/(p_rf_2piRp*2*pi);
0643 ksirf = 32.7*I_tot.*P_rf_ne*1e-20./(P_rf_Te*2*pi.*p_rf_2piRp);
0644
0645 eta = NaN(1,nr);
0646 eta(P_tot ~= 0) = J_tot(P_tot ~= 0)./P_tot(P_tot ~= 0);
0647
0648 I_exp = NaN;
0649
0650 if isfield(dke_out,'equil') && isfield(dke_out.equil,'icd'),
0651
0652 if p_e_2piRp == 0,
0653
0654 I_exp = dke_out.equil.icd/1e6;
0655
0656 elseif isfield(dke_out.equil,'iohm')
0657
0658 if sum(P_rf) == 0,
0659 I_exp = dke_out.equil.iohm/1e6;
0660 else
0661 I_exp = (dke_out.equil.icd + dke_out.equil.iohm)/1e6;
0662 end
0663
0664 else
0665
0666 disp('Warning : ohmic heating not available, exp. comparison discarded')
0667 end
0668 end
0669
0670 if timevol >= 1,
0671
0672 FPtime = dke_out.FPtime;
0673
0674 nit = length(FPtime);
0675
0676 if timevol == 1,
0677
0678 tJ_tot = reshape([Zcurr(end-nit+1:end).x_0_fsav],[nr,nit])*mksa.j_ref*scocurr;
0679 tP_rf = [ZP0(end-nit+1:end).x_rf_fsav]*mksa.P_ref;
0680 tP_e = reshape([ZP0(end-nit+1:end).x_e_fsav],[nr,nit])*mksa.P_ref;
0681
0682 elseif timevol == 2,
0683
0684 tJ_tot = NaN(nr,nit);
0685 tP_rf = NaN(nr,nit);
0686 tP_e = NaN(nr,nit);
0687
0688 for itn = 1:nit,
0689
0690 it_rf = length(dke_out.residue_rf{itn});
0691
0692 tJ_tot(:,itn) = Zcurr(it_rf,itn).x_0_fsav.'*mksa.j_ref*scocurr;
0693 tP_rf(:,itn) = ZP0(it_rf,itn).x_rf_fsav*mksa.P_ref;
0694 tP_e(:,itn) = ZP0(it_rf,itn).x_e_fsav.'*mksa.P_ref;
0695
0696 end
0697
0698 J_tot_rf = reshape([Zcurr(1:nit_rf,end).x_0_fsav],[nr,nit_rf])*mksa.j_ref*scocurr;
0699 I_tot_rf = sum(J_tot_rf.*repmat(equilDKE.xdA_dke.',[1,nit_rf]));
0700
0701 end
0702
0703 if ~isempty(waves)
0704 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
0705 tP_tot = tP_coll_luke + tP_e;
0706 else
0707 tP_tot = tP_rf + P_e;
0708 end
0709 else
0710 tP_tot = tP_e;
0711 end
0712
0713 tI_tot = sum(tJ_tot.*repmat(equilDKE.xdA_dke.',[1,nit]));
0714
0715 xtdV_2piRp_dke = repmat(equilDKE.xdV_2piRp_dke.',[1,nit]);
0716
0717 tp_rf_2piRp = sum(tP_rf.*xtdV_2piRp_dke);
0718 tp_e_2piRp = sum(tP_e.*xtdV_2piRp_dke);
0719 tp_tot_2piRp = sum(tP_tot.*xtdV_2piRp_dke);
0720
0721 end
0722
0723 if iscell(opt_blist),
0724
0725 nw = length(opt_blist);
0726
0727 nb = NaN(1,nw);
0728 ymask = cell(1,nw);
0729 wp0_rf_2piRp = zeros(1,nw);
0730
0731 wave = waves{1};
0732 nbw = length(wave.rays);
0733
0734 for iw = 1:nw,
0735
0736 waves{iw} = wave;
0737 waves{iw}.rays = wave.rays(opt_blist{iw});
0738
0739 nb(iw) = length(opt_blist{iw});
0740 ymask{iw} = [];
0741
0742 for iiray = 1:nb(iw),
0743
0744 iray = opt_blist{iw}(iiray);
0745
0746 wp0_rf_2piRp(iw) = wp0_rf_2piRp(iw) + waves{iw}.rays{iiray}.P0_2piRp*1e-6;
0747
0748 if iray < nbw
0749 ymaskloc = yb(iray):(yb(iray+1) - 1);
0750 else
0751 ymaskloc = yb(iray):ny;
0752 end
0753
0754 ymask{iw} = sort([ymask{iw},ymaskloc]);
0755
0756 end
0757
0758 end
0759
0760 clear wave
0761
0762 else
0763
0764 nb = NaN(1,nw);
0765 ymask = cell(1,nw);
0766 wp0_rf_2piRp = zeros(1,nw);
0767
0768 for iw = 1:nw,
0769
0770 if ~isfield(waves{iw},'model') || strcmp(waves{iw}.model,'RT'),
0771
0772 nb(iw) = length(waves{iw}.rays);
0773 for iray = 1:nb(iw),
0774
0775 wp0_rf_2piRp(iw) = wp0_rf_2piRp(iw) + waves{iw}.rays{iray}.P0_2piRp*1e-6;
0776
0777 end
0778
0779 elseif strcmp(waves{iw}.model,'FW'),
0780
0781 nb(iw) = 0;
0782 if isfield(waves{iw},'P0_2piRp'),
0783 wp0_rf_2piRp(iw) = waves{iw}.P0_2piRp*1e-6;
0784 else
0785 wp0_rf_2piRp(iw) = NaN;
0786 end
0787
0788 else
0789 error('Wave type not recognized.')
0790 end
0791
0792 ibmin = sum(nb(1:iw-1)) + 1;
0793 ibmax = sum(nb(1:iw)) + 1;
0794
0795 if iw < nw
0796 ymask{iw} = yb(ibmin):(yb(ibmax) - 1);
0797 else
0798 ymask{iw} = yb(ibmin):ny;
0799 end
0800
0801 end
0802 end
0803
0804
0805
0806 if isfield(dke_out,'yP0_2piRp'),
0807 p0_rf_2piRp = sum(dke_out.yP0_2piRp)/1e6;
0808 else
0809 p0_rf_2piRp = sum(wp0_rf_2piRp);
0810 end
0811
0812 if nr > 1,
0813
0814 if opt.wave == 'y',
0815
0816
0817
0818 wxP_rf_lin = zeros(nw,nr);
0819 wxP_rf_lin_coll = zeros(nw,nr);
0820 wxP_rf_coll_luke = zeros(nw,nr);
0821
0822 wP_rf_rhoG_lin = zeros(1,nw);
0823 wP_rf_rhoP_lin = zeros(1,nw);
0824 wP_rf_rhoT_lin = zeros(1,nw);
0825 wP_rf_rhoV_lin = zeros(1,nw);
0826
0827 wP_rf_rhoG_lin_coll = zeros(1,nw);
0828 wP_rf_rhoP_lin_coll = zeros(1,nw);
0829 wP_rf_rhoT_lin_coll = zeros(1,nw);
0830 wP_rf_rhoV_lin_coll = zeros(1,nw);
0831
0832 wP_rf_rhoG_coll_luke = zeros(1,nw);
0833 wP_rf_rhoP_coll_luke = zeros(1,nw);
0834 wP_rf_rhoT_coll_luke = zeros(1,nw);
0835 wP_rf_rhoV_coll_luke = zeros(1,nw);
0836
0837 for iw = 1:nw,
0838
0839 if ~isfield(waves{iw},'model') || strcmp(waves{iw}.model,'RT'),
0840
0841 for iray = 1:nb(iw),
0842
0843 ray = waves{iw}.rays{iray};
0844
0845 if isfield(ray,'sP_2piRp_lin') && isfield(ray,'srho'),
0846 wxP_rf_lin(iw,:) = wxP_rf_lin(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin,xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0847 else
0848 wxP_rf_lin(iw,:) = NaN;
0849 end
0850
0851 if isfield(ray,'sP_2piRp_lin_coll') && isfield(ray,'srho'),
0852 wxP_rf_lin_coll(iw,:) = wxP_rf_lin_coll(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin_coll,xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0853 else
0854 wxP_rf_lin_coll(iw,:) = NaN;
0855 end
0856
0857 if isfield(ray,'sP_2piRp_coll_luke') && isfield(ray,'srho'),
0858 wxP_rf_coll_luke(iw,:) = wxP_rf_coll_luke(iw,:) + 1e-6*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_coll_luke{1},xrhoG_S).*xdrhoG./equilDKE.xdV_2piRp_dke;
0859 else
0860 wxP_rf_coll_luke(iw,:) = NaN;
0861 end
0862
0863 end
0864
0865 elseif strcmp(waves{iw}.model,'FW'),
0866
0867 nr_fw = length(waves{iw}.rho);
0868 xwmask = zeros(nr,nr_fw);
0869
0870 for ir_fw = 1:nr_fw,
0871 if waves{iw}.rhotype == 'g',
0872 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoG_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoG_S(1:end-1)]))]).'./xdrhoG.';
0873 elseif waves{iw}.rhotype == 'p',
0874 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoP_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoP_S(1:end-1)]))]).'./xdrhoP.';
0875 elseif waves{iw}.rhotype == 't',
0876 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoT_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoT_S(1:end-1)]))]).'./xdrhoT.';
0877 elseif waves{iw}.rhotype == 'v',
0878 xwmask(:,ir_fw) = max([zeros(1,nr);(min([repmat(waves{iw}.rho_S(ir_fw + 1),[1,nr]);xrhoV_S(2:end)]) - max([repmat(waves{iw}.rho_S(ir_fw),[1,nr]);xrhoV_S(1:end-1)]))]).'./xdrhoV.';
0879 else
0880 error('rho type not recognized')
0881 end
0882 end
0883
0884 wpe_eve = 1e-6*(waves{iw}.pe_eve(1:end - 1) + waves{iw}.pe_eve(2:end))/2;
0885
0886 wxP_rf_lin(iw,:) = (xwmask*wpe_eve).';
0887
0888 else
0889 error('Wave type not recognized.')
0890 end
0891
0892 [wP_rf_rhoG_lin(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0893 [wP_rf_rhoP_lin(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0894 [wP_rf_rhoT_lin(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0895 [wP_rf_rhoV_lin(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_lin(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0896
0897 [wP_rf_rhoG_lin_coll(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0898 [wP_rf_rhoP_lin_coll(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0899 [wP_rf_rhoT_lin_coll(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0900 [wP_rf_rhoV_lin_coll(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_lin_coll(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0901
0902 [wP_rf_rhoG_coll_luke(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0903 [wP_rf_rhoP_coll_luke(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0904 [wP_rf_rhoT_coll_luke(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0905 [wP_rf_rhoV_coll_luke(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf_coll_luke(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0906 end
0907
0908 wxp_rf_2piRp_lin = wxP_rf_lin.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0909 wp_rf_2piRp_lin = sum(wxp_rf_2piRp_lin,2);
0910 xp_rf_2piRp_lin = sum(wxp_rf_2piRp_lin,1);
0911 p_rf_2piRp_lin = sum(wp_rf_2piRp_lin);
0912
0913 wxp_rf_2piRp_lin_coll = wxP_rf_lin_coll.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0914 wp_rf_2piRp_lin_coll = sum(wxp_rf_2piRp_lin_coll,2);
0915 xp_rf_2piRp_lin_coll = sum(wxp_rf_2piRp_lin_coll,1);
0916 p_rf_2piRp_lin_coll = sum(wp_rf_2piRp_lin_coll);
0917
0918 wxp_rf_2piRp_coll_luke = wxP_rf_coll_luke.*repmat(equilDKE.xdV_2piRp_dke,[nw,1]);
0919 wp_rf_2piRp_coll_luke = sum(wxp_rf_2piRp_coll_luke,2);
0920 xp_rf_2piRp_coll_luke = sum(wxp_rf_2piRp_coll_luke,1);
0921 p_rf_2piRp_coll_luke = sum(wp_rf_2piRp_coll_luke);
0922
0923 if ~all(isnan(wP_rf_rhoG_lin)),
0924 wTe_lin = interp1(xrhoG,xTe,wP_rf_rhoG_lin);
0925 wne_lin = interp1(xrhoG,xne,wP_rf_rhoG_lin);
0926 else
0927 wTe_lin = NaN(1,nw);
0928 wne_lin = NaN(1,nw);
0929 end
0930
0931 if ~all(isnan(wP_rf_rhoG_lin_coll)),
0932 wTe_lin_coll = interp1(xrhoG,xTe,wP_rf_rhoG_lin_coll);
0933 wne_lin_coll = interp1(xrhoG,xne,wP_rf_rhoG_lin_coll);
0934 else
0935 wTe_lin_coll = NaN(1,nw);
0936 wne_lin_coll = NaN(1,nw);
0937 end
0938
0939 if ~all(isnan(wP_rf_rhoG_coll_luke)),
0940 wTe_coll_luke = interp1(xrhoG,xTe,wP_rf_rhoG_coll_luke);
0941 wne_coll_luke = interp1(xrhoG,xne,wP_rf_rhoG_coll_luke);
0942 else
0943 wTe_coll_luke = NaN(1,nw);
0944 wne_coll_luke = NaN(1,nw);
0945 end
0946
0947
0948
0949 wxP_rf = zeros(nw,nr);
0950
0951 wxp_rf_2piRp = zeros(nw,nr);
0952 wnp_rf_2piRp = zeros(nw,nn_rf);
0953
0954 wP_rf_rhoG = zeros(1,nw);
0955 wP_rf_rhoP = zeros(1,nw);
0956 wP_rf_rhoT = zeros(1,nw);
0957 wP_rf_rhoV = zeros(1,nw);
0958
0959 wids = cell(1,nw);
0960
0961 wrays = cell(1,nw);
0962
0963 for iw = 1:nw,
0964
0965 wids{iw} = waves{iw}.id;
0966
0967 wxP_rf(iw,:) = sum(yP_rf(ymask{iw},:),1);
0968
0969 wxp_rf_2piRp(iw,:) = sum(yxp_rf_2piRp(ymask{iw},:),1);
0970 wnp_rf_2piRp(iw,:) = reshape(sum(sum(yxnp_rf_2piRp(:,ymask{iw},:),2),1),[1,nn_rf]);
0971
0972 [wP_rf_rhoG(iw)] = calc_peak_jd(xrhoG,xdrhoG,wxP_rf(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0973 [wP_rf_rhoP(iw)] = calc_peak_jd(xrhoP,xdrhoP,wxP_rf(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0974 [wP_rf_rhoT(iw)] = calc_peak_jd(xrhoT,xdrhoT,wxP_rf(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0975 [wP_rf_rhoV(iw)] = calc_peak_jd(xrhoV,xdrhoV,wxP_rf(iw,:),equilDKE.xdV_2piRp_dke,opt.jpeak,opt.peakmode);
0976
0977 if ~isfield(waves{iw},'model') || strcmp(waves{iw}.model,'RT'),
0978
0979 wrays{iw} = cell(1,nb(iw));
0980
0981 for iray = 1:nb(iw),
0982
0983 wrays{iw}{iray} = proc_ray_jd(waves{iw}.rays{iray});
0984
0985 end
0986
0987 elseif strcmp(waves{iw}.model,'FW'),
0988
0989 else
0990 error('Wave type not recognized.')
0991 end
0992 end
0993
0994 if ~all(isnan(wP_rf_rhoG)),
0995 wTe = interp1(xrhoG,xTe,wP_rf_rhoG);
0996 wne = interp1(xrhoG,xne,wP_rf_rhoG);
0997 else
0998 wTe = NaN(1,nw);
0999 wne = NaN(1,nw);
1000 end
1001
1002 wp_rf_2piRp = sum(wxp_rf_2piRp,2);
1003
1004 end
1005
1006 if opt.peak == 'y',
1007
1008
1009
1010
1011 dep = 1;
1012 xP_rf = P_rf;
1013 xP_rfm = 0;
1014
1015 xP_rf(xP_rf < 0) = 0;
1016
1017 p_rf_2piRp_peak = 0;
1018 I_tot_peak = 0;
1019
1020 peaks.ir = [];
1021 peaks.p_rf_2piRp = [];
1022 peaks.I_tot = [];
1023
1024 for ir = 1:nr,
1025 if (xP_rf(ir) - xP_rfm)*dep < 0,
1026 if dep == 1,
1027 peaks.ir = [peaks.ir,ir-1];
1028 dep = -1;
1029 else
1030
1031 p_rf_2piRp_peak = p_rf_2piRp_peak - xp_rf_2piRp(ir-1)/2;
1032 I_tot_peak = I_tot_peak - xI(ir-1)/2;
1033
1034 peaks.p_rf_2piRp = [peaks.p_rf_2piRp,p_rf_2piRp_peak];
1035 peaks.I_tot = [peaks.I_tot,I_tot_peak];
1036
1037 p_rf_2piRp_peak = xp_rf_2piRp(ir-1)/2;
1038 I_tot_peak = xI(ir-1)/2;
1039
1040 dep = 1;
1041 end
1042 end
1043
1044 xP_rfm = xP_rf(ir);
1045 p_rf_2piRp_peak = p_rf_2piRp_peak + xp_rf_2piRp(ir);
1046 I_tot_peak = I_tot_peak + xI(ir);
1047 end
1048
1049 if dep == 1,
1050 peaks.ir = [peaks.ir,nr];
1051 end
1052 peaks.p_rf_2piRp = [peaks.p_rf_2piRp,p_rf_2piRp_peak];
1053 peaks.I_tot = [peaks.I_tot,I_tot_peak];
1054
1055 peaks.rhoG = xrhoG(peaks.ir);
1056 peaks.rhoP = xrhoP(peaks.ir);
1057 peaks.rhoT = xrhoT(peaks.ir);
1058 peaks.rhoV = xrhoV(peaks.ir);
1059 peaks.Te = xTe(peaks.ir);
1060 peaks.ne = xne(peaks.ir);
1061 peaks.EtaRp = peaks.I_tot./(peaks.p_rf_2piRp*2*pi);
1062 peaks.ksirf = 32.7*peaks.I_tot.*peaks.ne*1e-20./(peaks.Te*2*pi.*peaks.p_rf_2piRp);
1063
1064 else
1065 peaks = '';
1066 end
1067
1068
1069
1070 data_proc.simul = dke_out.simul;
1071 data_proc.simul.residue_rf = residue_rf;
1072
1073
1074
1075 data_proc.scalar.p_tot_2piRp = p_tot_2piRp;
1076
1077 data_proc.scalar.p0_rf_2piRp = p0_rf_2piRp;
1078 data_proc.scalar.p_rf_2piRp = p_rf_2piRp;
1079 data_proc.scalar.p_rf_2piRp_lin = p_rf_2piRp_lin;
1080 data_proc.scalar.p_rf_2piRp_lin_coll = p_rf_2piRp_lin_coll;
1081 data_proc.scalar.p_rf_2piRp_coll_luke = p_rf_2piRp_coll_luke;
1082
1083 data_proc.scalar.np_rf_2piRp = np_rf_2piRp;
1084 data_proc.scalar.n_rf_list = n_rf_list;
1085
1086 data_proc.scalar.p_e_2piRp = p_e_2piRp;
1087 data_proc.scalar.I_tot = I_tot;
1088 data_proc.scalar.I_exp = I_exp;
1089 data_proc.scalar.EtaRp = EtaRp;
1090 data_proc.scalar.EtaRFRp = EtaRFRp;
1091 data_proc.scalar.EtaRFRp_nebar0 = EtaRFRp_nebar0;
1092 data_proc.scalar.ksirf = ksirf;
1093 data_proc.scalar.nfrags = nfrags;
1094 data_proc.scalar.Rp = equilDKE.Rp;
1095 data_proc.scalar.nebar0 = ne_bar0;
1096
1097 if isfield(dke_out,'xepsi_init'),
1098 data_proc.scalar.V_loop_init = 2*pi*equilDKE.Rp*mksa.Edreicer_ref*dke_out.xepsi_init(end);
1099 if isfield(data.Zbouncecoef,'XXlambda_b_p1m2p2')
1100 data_proc.scalar.V_loop_init_fsav = 2*pi*equilDKE.Rp*mksa.Edreicer_ref*dke_out.xepsi_init(end)/(data.Zbouncecoef.xqhat(end)/data.Zbouncecoef.xqtilde(end))/data.Zbouncecoef.XXlambda_b_p1m2p2(1,1,end);
1101 else
1102 data_proc.scalar.V_loop_init_fsav = NaN;
1103 end
1104 else
1105 data_proc.scalar.V_loop_init = NaN;
1106 data_proc.scalar.V_loop_init_fsav = NaN;
1107 end
1108
1109 data_proc.scalar.jrhoG = jrhoG;
1110 data_proc.scalar.jrhoP = jrhoP;
1111 data_proc.scalar.jrhoT = jrhoT;
1112 data_proc.scalar.jrhoV = jrhoV;
1113 data_proc.scalar.jdrhoG = jdrhoG;
1114 data_proc.scalar.jdrhoP = jdrhoP;
1115 data_proc.scalar.jdrhoT = jdrhoT;
1116 data_proc.scalar.jdrhoV = jdrhoV;
1117 data_proc.scalar.jTe = jTe;
1118 data_proc.scalar.jne = jne;
1119
1120 data_proc.scalar.P_rf_rhoG = P_rf_rhoG;
1121 data_proc.scalar.P_rf_rhoP = P_rf_rhoP;
1122 data_proc.scalar.P_rf_rhoT = P_rf_rhoT;
1123 data_proc.scalar.P_rf_rhoV = P_rf_rhoV;
1124 data_proc.scalar.P_rf_drhoG = P_rf_drhoG;
1125 data_proc.scalar.P_rf_drhoP = P_rf_drhoP;
1126 data_proc.scalar.P_rf_drhoT = P_rf_drhoT;
1127 data_proc.scalar.P_rf_drhoV = P_rf_drhoV;
1128 data_proc.scalar.P_rf_Te = P_rf_Te;
1129 data_proc.scalar.P_rf_ne = P_rf_ne;
1130
1131 data_proc.scalar.P_rf_rhoG_lin = P_rf_rhoG_lin;
1132 data_proc.scalar.P_rf_rhoP_lin = P_rf_rhoP_lin;
1133 data_proc.scalar.P_rf_rhoT_lin = P_rf_rhoT_lin;
1134 data_proc.scalar.P_rf_rhoV_lin = P_rf_rhoV_lin;
1135 data_proc.scalar.P_rf_drhoG_lin = P_rf_drhoG_lin;
1136 data_proc.scalar.P_rf_drhoP_lin = P_rf_drhoP_lin;
1137 data_proc.scalar.P_rf_drhoT_lin = P_rf_drhoT_lin;
1138 data_proc.scalar.P_rf_drhoV_lin = P_rf_drhoV_lin;
1139 data_proc.scalar.P_rf_Te_lin = P_rf_Te_lin;
1140 data_proc.scalar.P_rf_ne_lin = P_rf_ne_lin;
1141
1142 data_proc.scalar.P_rf_rhoG_lin_coll = P_rf_rhoG_lin_coll;
1143 data_proc.scalar.P_rf_rhoP_lin_coll = P_rf_rhoP_lin_coll;
1144 data_proc.scalar.P_rf_rhoT_lin_coll = P_rf_rhoT_lin_coll;
1145 data_proc.scalar.P_rf_rhoV_lin_coll = P_rf_rhoV_lin_coll;
1146 data_proc.scalar.P_rf_drhoG_lin_coll = P_rf_drhoG_lin_coll;
1147 data_proc.scalar.P_rf_drhoP_lin_coll = P_rf_drhoP_lin_coll;
1148 data_proc.scalar.P_rf_drhoT_lin_coll = P_rf_drhoT_lin_coll;
1149 data_proc.scalar.P_rf_drhoV_lin_coll = P_rf_drhoV_lin_coll;
1150 data_proc.scalar.P_rf_Te_lin_coll = P_rf_Te_lin_coll;
1151 data_proc.scalar.P_rf_ne_lin_coll = P_rf_ne_lin_coll;
1152
1153 data_proc.scalar.P_rf_rhoG_coll_luke = P_rf_rhoG_coll_luke;
1154 data_proc.scalar.P_rf_rhoP_coll_luke = P_rf_rhoP_coll_luke;
1155 data_proc.scalar.P_rf_rhoT_coll_luke = P_rf_rhoT_coll_luke;
1156 data_proc.scalar.P_rf_rhoV_coll_luke = P_rf_rhoV_coll_luke;
1157 data_proc.scalar.P_rf_drhoG_coll_luke = P_rf_drhoG_coll_luke;
1158 data_proc.scalar.P_rf_drhoP_coll_luke = P_rf_drhoP_coll_luke;
1159 data_proc.scalar.P_rf_drhoT_coll_luke = P_rf_drhoT_coll_luke;
1160 data_proc.scalar.P_rf_drhoV_coll_luke = P_rf_drhoV_coll_luke;
1161 data_proc.scalar.P_rf_Te_coll_luke = P_rf_Te_coll_luke;
1162 data_proc.scalar.P_rf_ne_coll_luke = P_rf_ne_coll_luke;
1163
1164
1165
1166 data_proc.radial.xrhoG = xrhoG;
1167 data_proc.radial.xrhoP = xrhoP;
1168 data_proc.radial.xrhoT = xrhoT;
1169 data_proc.radial.xrhoV = xrhoV;
1170 data_proc.radial.xdrhoG = xdrhoG;
1171 data_proc.radial.xdrhoP = xdrhoP;
1172 data_proc.radial.xdrhoT = xdrhoT;
1173 data_proc.radial.xdrhoV = xdrhoV;
1174
1175 data_proc.radial.xdA_dke = equilDKE.xdA_dke;
1176 data_proc.radial.xdV_2piRp_dke = equilDKE.xdV_2piRp_dke;
1177
1178 if isfield(dke_out,'xepsi'),
1179 data_proc.radial.xepsi = dke_out.xepsi;
1180 end
1181
1182 data_proc.radial.J_tot = J_tot;
1183 data_proc.radial.P_tot = P_tot;
1184 data_proc.radial.eta = eta;
1185 data_proc.radial.P_rf = P_rf;
1186 data_proc.radial.P_e = P_e;
1187 data_proc.radial.xdpdt_rf = xdpdt_rf;
1188
1189 data_proc.radial.xnp_rf_2piRp = xnp_rf_2piRp;
1190
1191 data_proc.radial.xp_rf_2piRp = xp_rf_2piRp;
1192 data_proc.radial.xp_rf_2piRp_lin = xp_rf_2piRp_lin;
1193 data_proc.radial.xp_rf_2piRp_lin_coll = xp_rf_2piRp_lin_coll;
1194 data_proc.radial.xp_rf_2piRp_coll_luke = xp_rf_2piRp_coll_luke;
1195
1196 data_proc.scalar.li = internal_inductance_yp(data_proc.radial.xrhoG,data_proc.radial.J_tot);
1197
1198
1199
1200 if timevol >= 1,
1201 data_proc.timevol.FPtime = FPtime;
1202 data_proc.timevol.tJ_tot = tJ_tot;
1203 data_proc.timevol.tP_rf = tP_rf;
1204 data_proc.timevol.tP_e = tP_e;
1205 data_proc.timevol.tP_tot = tP_tot;
1206 data_proc.timevol.tI_tot = tI_tot;
1207 data_proc.timevol.tp_rf_2piRp = tp_rf_2piRp;
1208 data_proc.timevol.tp_e_2piRp = tp_e_2piRp;
1209 data_proc.timevol.tp_tot_2piRp = tp_tot_2piRp;
1210 end
1211
1212 if timevol >= 2,
1213 data_proc.timevol.J_tot_rf = J_tot_rf;
1214 data_proc.timevol.I_tot_rf = I_tot_rf;
1215 end
1216
1217
1218
1219 if opt.wave == 'y',
1220 data_proc.wave.wids = wids;
1221
1222 data_proc.wave.wxP_rf = wxP_rf;
1223 data_proc.wave.wp0_rf_2piRp = wp0_rf_2piRp;
1224 data_proc.wave.wp_rf_2piRp = wp_rf_2piRp;
1225 data_proc.wave.wnp_rf_2piRp = wnp_rf_2piRp;
1226 data_proc.wave.wn_rf_list = n_rf_list;
1227 data_proc.wave.wxp_rf_2piRp = wxp_rf_2piRp;
1228
1229 data_proc.wave.wP_rf_rhoG = wP_rf_rhoG;
1230 data_proc.wave.wP_rf_rhoP = wP_rf_rhoP;
1231 data_proc.wave.wP_rf_rhoT = wP_rf_rhoT;
1232 data_proc.wave.wP_rf_rhoV = wP_rf_rhoV;
1233
1234 data_proc.wave.wTe = wTe;
1235 data_proc.wave.wne = wne;
1236
1237 data_proc.wave.wrays = wrays;
1238
1239
1240
1241 data_proc.wave.wxP_rf_lin = wxP_rf_lin;
1242 data_proc.wave.wp_rf_2piRp_lin = wp_rf_2piRp_lin;
1243 data_proc.wave.wxp_rf_2piRp_lin = wxp_rf_2piRp_lin;
1244
1245 data_proc.wave.wP_rf_rhoG_lin = wP_rf_rhoG_lin;
1246 data_proc.wave.wP_rf_rhoP_lin = wP_rf_rhoP_lin;
1247 data_proc.wave.wP_rf_rhoT_lin = wP_rf_rhoT_lin;
1248 data_proc.wave.wP_rf_rhoV_lin = wP_rf_rhoV_lin;
1249
1250 data_proc.wave.wTe_lin = wTe_lin;
1251 data_proc.wave.wne_lin = wne_lin;
1252
1253
1254
1255 data_proc.wave.wxP_rf_lin_coll = wxP_rf_lin_coll;
1256 data_proc.wave.wp_rf_2piRp_lin_coll = wp_rf_2piRp_lin_coll;
1257 data_proc.wave.wxp_rf_2piRp_lin_coll = wxp_rf_2piRp_lin_coll;
1258
1259 data_proc.wave.wP_rf_rhoG_lin_coll = wP_rf_rhoG_lin_coll;
1260 data_proc.wave.wP_rf_rhoP_lin_coll = wP_rf_rhoP_lin_coll;
1261 data_proc.wave.wP_rf_rhoT_lin_coll = wP_rf_rhoT_lin_coll;
1262 data_proc.wave.wP_rf_rhoV_lin_coll = wP_rf_rhoV_lin_coll;
1263
1264 data_proc.wave.wTe_lin_coll = wTe_lin_coll;
1265 data_proc.wave.wne_lin_coll = wne_lin_coll;
1266
1267 data_proc.wave.wxP_rf_coll_luke = wxP_rf_coll_luke;
1268 data_proc.wave.wp_rf_2piRp_coll_luke = wp_rf_2piRp_coll_luke;
1269 data_proc.wave.wxp_rf_2piRp_coll_luke = wxp_rf_2piRp_coll_luke;
1270
1271 data_proc.wave.wP_rf_rhoG_coll_luke = wP_rf_rhoG_coll_luke;
1272 data_proc.wave.wP_rf_rhoP_coll_luke = wP_rf_rhoP_coll_luke;
1273 data_proc.wave.wP_rf_rhoT_coll_luke = wP_rf_rhoT_coll_luke;
1274 data_proc.wave.wP_rf_rhoV_coll_luke = wP_rf_rhoV_coll_luke;
1275
1276 data_proc.wave.wTe_coll_luke = wTe_coll_luke;
1277 data_proc.wave.wne_coll_luke = wne_coll_luke;
1278
1279 else
1280 data_proc.wave = '';
1281 end
1282
1283
1284
1285 data_proc.peaks = peaks;
1286
1287
1288
1289 data_proc.opt = opt;
1290
1291 else
1292
1293 data_proc = struct;
1294
1295 end
1296 end
1297
1298 if isfield(data_proc,'radial')
1299 if opt.rho == 'g',
1300 xrho = data_proc.radial.xrhoG;
1301 dxrho = data_proc.radial.xdrhoG;
1302 xlab_rho = '\rho_G';
1303 elseif opt.rho == 'p',
1304 xrho = data_proc.radial.xrhoP;
1305 dxrho = data_proc.radial.xdrhoP;
1306 xlab_rho = '\rho_P';
1307 elseif opt.rho == 't',
1308 xrho = data_proc.radial.xrhoT;
1309 dxrho = data_proc.radial.xdrhoT;
1310 xlab_rho = '\rho_T';
1311 elseif opt.rho == 'v',
1312 xrho = data_proc.radial.xrhoV;
1313 dxrho = data_proc.radial.xdrhoV;
1314 xlab_rho = '\rho_V';
1315 end
1316 else
1317 if opt.rho == 'g',
1318 xrho = 0;
1319 dxrho = 0;
1320 xlab_rho = '\rho_G';
1321 elseif opt.rho == 'p',
1322 xrho = 0;
1323 dxrho = 0;
1324 xlab_rho = '\rho_P';
1325 elseif opt.rho == 't',
1326 xrho = 0;
1327 dxrho = 0;
1328 xlab_rho = '\rho_T';
1329 elseif opt.rho == 'v',
1330 xrho = 0;
1331 dxrho = 0;
1332 xlab_rho = '\rho_V';
1333 end
1334 end
1335
1336
1337 if isnan(ir_display) || ir_display >= 0,
1338
1339 if isfield(opt,'diaryname') && ischar(opt.diaryname),
1340 diaryname = opt.diaryname;
1341 else
1342 if isfield(data_proc,'simul'),
1343 diaryname = ['res_',data_proc.simul.id,'.log'];
1344 else
1345 diaryname = '';
1346 end
1347
1348 if exist(diaryname,'file') && gettime_jd(diaryname) > datenum(timeid_jd(data_proc.simul.timeid),'mm/dd/yyyy HH:MM:SS'),
1349 diaryname = '';
1350 end
1351
1352 diaryname = input_dke_yp('Provide a name for the log file? (return empty string if no log file)',diaryname);
1353 end
1354
1355 if ~isempty(diaryname),
1356 if exist(diaryname,'file'),
1357 delete(diaryname)
1358 end
1359 diary(diaryname)
1360 end
1361
1362 if isfield(data_proc,'simul'),
1363 disp(['Simulation id: ',data_proc.simul.id]);
1364 disp(['Simulation date: ',timeid_jd(data_proc.simul.timeid)]);
1365 disp(['User id: ',data_proc.simul.userid]);
1366 disp(['LUKE Version: ',data_proc.simul.LUKEversion]);
1367 if all(isfield(data_proc.simul,{'matver','hostname'})),
1368 disp(['MATLAB Version: ',data_proc.simul.matver]);
1369 disp(['MACHINE: ',data_proc.simul.hostname]);
1370 end
1371 disp(['Final RF convergence parameter: ',num2str(data_proc.simul.residue_rf(end))]);
1372 disp(' ');
1373 else
1374 data_proc.simul.id = '';
1375 data_proc.simul.timeid = '';
1376 data_proc.simul.userid = '';
1377 data_proc.simul.LUKEversion = '';
1378 data_proc.simul.matver = '';
1379 data_proc.simul.hostname = '';
1380 data_proc.simul.residue_rf = '';
1381
1382 red = 1;
1383 lspace = NaN;
1384 lspace2 = NaN;
1385 bspace = NaN;
1386 bspace2 = NaN;
1387 end
1388
1389 if nr > 1,
1390
1391 disp(['Total current (I): ',num2str(data_proc.scalar.I_tot),' MA (positive for co-current)'])
1392 disp(['Exp. current (I): ',num2str(data_proc.scalar.I_exp),' MA (positive for co-current)'])
1393
1394 if opt.rho == 'g',
1395 jrho = data_proc.scalar.jrhoG;
1396 jdrho = data_proc.scalar.jdrhoG;
1397
1398 if ~isempty(waves),
1399 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
1400 P_rf_rho_lin = data_proc.scalar.P_rf_rhoG_lin_coll;
1401 P_rf_drho_lin = data_proc.scalar.P_rf_drhoG_lin_coll;
1402 P_rf_rho = data_proc.scalar.P_rf_rhoG_coll_luke;
1403 P_rf_drho = data_proc.scalar.P_rf_drhoG_coll_luke;
1404 else
1405 P_rf_rho_lin = data_proc.scalar.P_rf_rhoG_lin;
1406 P_rf_drho_lin = data_proc.scalar.P_rf_drhoG_lin;
1407 P_rf_rho = data_proc.scalar.P_rf_rhoG;
1408 P_rf_drho = data_proc.scalar.P_rf_drhoG;
1409 end
1410 end
1411
1412 elseif opt.rho == 'p',
1413 jrho = data_proc.scalar.jrhoP;
1414 jdrho = data_proc.scalar.jdrhoP;
1415
1416 if ~isempty(waves),
1417 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
1418 P_rf_rho_lin = data_proc.scalar.P_rf_rhoP_lin_coll;
1419 P_rf_drho_lin = data_proc.scalar.P_rf_drhoP_lin_coll;
1420 P_rf_rho = data_proc.scalar.P_rf_rhoP_coll_luke;
1421 P_rf_drho = data_proc.scalar.P_rf_drhoP_coll_luke;
1422 else
1423 P_rf_rho_lin = data_proc.scalar.P_rf_rhoP_lin;
1424 P_rf_drho_lin = data_proc.scalar.P_rf_drhoP_lin;
1425 P_rf_rho = data_proc.scalar.P_rf_rhoP;
1426 P_rf_drho = data_proc.scalar.P_rf_drhoP;
1427 end
1428 end
1429
1430 elseif opt.rho == 't',
1431 jrho = data_proc.scalar.jrhoT;
1432 jdrho = data_proc.scalar.jdrhoT;
1433
1434 if ~isempty(waves),
1435 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
1436 P_rf_rho_lin = data_proc.scalar.P_rf_rhoT_lin_coll;
1437 P_rf_drho_lin = data_proc.scalar.P_rf_drhoT_lin_coll;
1438 P_rf_rho = data_proc.scalar.P_rf_rhoT_coll_luke;
1439 P_rf_drho = data_proc.scalar.P_rf_drhoT_coll_luke;
1440 else
1441 P_rf_rho_lin = data_proc.scalar.P_rf_rhoT_lin;
1442 P_rf_drho_lin = data_proc.scalar.P_rf_drhoT_lin;
1443 P_rf_rho = data_proc.scalar.P_rf_rhoT;
1444 P_rf_drho = data_proc.scalar.P_rf_drhoT;
1445 end
1446 end
1447
1448 elseif opt.rho == 'v',
1449 jrho = data_proc.scalar.jrhoV;
1450 jdrho = data_proc.scalar.jdrhoV;
1451
1452 if ~isempty(waves),
1453 if isfield(waves{1}.rayparam,'colldamp') && waves{1}.rayparam.colldamp == 1,
1454 P_rf_rho_lin = data_proc.scalar.P_rf_rhoV_lin_coll;
1455 P_rf_drho_lin = data_proc.scalar.P_rf_drhoV_lin_coll;
1456 P_rf_rho = data_proc.scalar.P_rf_rhoV_coll_luke;
1457 P_rf_drho = data_proc.scalar.P_rf_drhoV_coll_luke;
1458 else
1459 P_rf_rho_lin = data_proc.scalar.P_rf_rhoV_lin;
1460 P_rf_drho_lin = data_proc.scalar.P_rf_drhoV_lin;
1461 P_rf_rho = data_proc.scalar.P_rf_rhoV;
1462 P_rf_drho = data_proc.scalar.P_rf_drhoV;
1463 end
1464 end
1465
1466 end
1467 disp(['Location of current peak (rho',opt.rho,'): ',num2str(jrho),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1468 disp(['Current FWHM (drho',opt.rho,'): ',num2str(jdrho),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1469 if sum(data_proc.radial.P_rf) ~= 0,
1470 disp(['RF efficiency (ne@jpeak*I*Rp/P): ',num2str(data_proc.scalar.EtaRFRp),' 10^19 A/(W.m^2)'])
1471 disp(['RF efficiency (ne_bar*I*Rp/P): ',num2str(data_proc.scalar.EtaRFRp_nebar0),' 10^19 A/(W.m^2)'])
1472 disp(['Location of RF power deposition peak (rho',opt.rho,'): ',num2str(P_rf_rho),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1473 disp(['RF power deposition FWHM (drho',opt.rho,'): ',num2str(P_rf_drho),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1474 disp(['Location of linear RF power deposition peak (rho',opt.rho,'): ',num2str(P_rf_rho_lin),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1475 disp(['Linear RF power deposition FWHM (drho',opt.rho,'): ',num2str(P_rf_drho_lin),' (peak mode = ',num2str(opt.peakmode),'; jpeak = ',num2str(opt.jpeak),')'])
1476 disp(['Number of RF ray fragments: ',num2str(data_proc.scalar.nfrags)])
1477
1478 end
1479
1480 disp(['Internal inductance : ',num2str(data_proc.scalar.li)])
1481 if isfield(dke_out.equil,'Te_vol')
1482 disp(['Volume averaged <Te>: ',num2str(dke_out.equil.Te_vol),' (keV)'])
1483 end
1484 disp(['Line-averaged density (Bmin) ne_bar: ',num2str(ne_bar0),' (10^19 m-3)'])
1485
1486 if isinf(data_proc.scalar.Rp) && sum(data_proc.radial.P_rf) ~= 0,
1487 disp(['Launched RF power (P/2piRp): ',num2str(data_proc.scalar.p0_rf_2piRp),' MW/m'])
1488 disp(['RF power absorbed (P/2piRp): ',num2str(data_proc.scalar.p_rf_2piRp),' MW/m ---> ',num2str(100*data_proc.scalar.p_rf_2piRp/data_proc.scalar.p0_rf_2piRp),' %.'])
1489 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1490 if dke_out.waves{1}.rayparam.colldamp == 1,
1491 disp(['RF power absorbed with NRCA (P/2piRp): ',num2str(data_proc.scalar.p_rf_2piRp_coll_luke),' MW/m ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_coll_luke/data_proc.scalar.p0_rf_2piRp),' %.'])
1492 else
1493 disp('Non-resonant collisional absorption (NRCA) of RF waves is not considered.')
1494 end
1495 else
1496 disp('Non-resonant collisional absorption (NRCA) of RF waves is not considered.')
1497 end
1498
1499 for in_rf = 1:nn_rf,
1500 disp(['RF power absorbed on harm. ',num2str(data_proc.scalar.n_rf_list(in_rf)),' (P/2piRp): ',num2str(data_proc.scalar.np_rf_2piRp(in_rf)),' MW/m ---> ',num2str(100*data_proc.scalar.np_rf_2piRp(in_rf)/data_proc.scalar.p0_rf_2piRp),' %.'])
1501 end
1502 if opt.wave == 'y',
1503 disp(['Linear RF power absorbed (P/2piRp): ',num2str(data_proc.scalar.p_rf_2piRp_lin),' MW/m ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_lin/data_proc.scalar.p0_rf_2piRp),' %.'])
1504 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1505 if dke_out.waves{1}.rayparam.colldamp == 1,
1506 disp(['Linear RF power absorbed with NRCA (P/2piRp): ',num2str(data_proc.scalar.p_rf_2piRp_lin_coll),' MW/m ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_lin_coll/data_proc.scalar.p0_rf_2piRp),' %.'])
1507 end
1508 end
1509 end
1510 disp(['Ohmic power deposited (P/2piRp): ',num2str(data_proc.scalar.p_e_2piRp),' MW/m'])
1511 disp(['Total power deposited (P/2piRp): ',num2str(data_proc.scalar.p_tot_2piRp),' MW/m'])
1512 disp(['CD efficiency (I*Rp/P): ',num2str(data_proc.scalar.EtaRp),' A.m/W'])
1513
1514 disp(['Normalized CD efficiency xi: ',num2str(data_proc.scalar.ksirf)])
1515
1516 if opt.peak == 'y',
1517
1518 for ipeak = 1:length(data_proc.peaks.ir),
1519 disp(' ');
1520 disp('----------------------------------------');
1521 disp(' ');
1522 disp(['RF Peak # ',num2str(ipeak)])
1523 disp(['Position rho G: ',num2str(data_proc.peaks.rhoG(ipeak))])
1524 disp(['Position rho P: ',num2str(data_proc.peaks.rhoP(ipeak))])
1525 disp(['Position rho T: ',num2str(data_proc.peaks.rhoT(ipeak))])
1526 disp(['Position rho V: ',num2str(data_proc.peaks.rhoV(ipeak))])
1527 disp(['Temperature Te: ',num2str(data_proc.peaks.Te(ipeak)),' keV'])
1528 disp(['Density ne: ',num2str(data_proc.peaks.ne(ipeak)/1e20),' *10^20 m^-3'])
1529 disp(' ');
1530 disp(['RF power (P/2piRp): ',num2str(data_proc.peaks.p_rf_2piRp(ipeak)),' MW/m'])
1531 disp(['Total current (I): ',num2str(data_proc.peaks.I_tot(ipeak)),' MA'])
1532 disp(['CD efficiency (I*Rp/P): ',num2str(data_proc.peaks.EtaRp(ipeak)),' A.m/W'])
1533 disp(['Normalized CD efficiency xi: ',num2str(data_proc.peaks.ksirf(ipeak))])
1534 end
1535 end
1536
1537 if opt.wave == 'y',
1538
1539 for iw = 1:nw,
1540 disp(' ');
1541 disp('----------------------------------------');
1542 disp(' ');
1543 disp(['RF Wave # ',num2str(iw)])
1544 disp(['Position rho G: ',num2str(data_proc.wave.wP_rf_rhoG(iw))])
1545 disp(['Position rho P: ',num2str(data_proc.wave.wP_rf_rhoP(iw))])
1546 disp(['Position rho T: ',num2str(data_proc.wave.wP_rf_rhoT(iw))])
1547 disp(['Position rho V: ',num2str(data_proc.wave.wP_rf_rhoV(iw))])
1548 disp(['Temperature Te: ',num2str(data_proc.wave.wTe(iw)),' keV'])
1549 disp(['Density ne: ',num2str(data_proc.wave.wne(iw)/1e20),' *10^20 m^-3'])
1550 disp(' ');
1551 disp(['Launched RF power (P/2piRp): ',num2str(data_proc.wave.wp0_rf_2piRp(iw)),'MW/m'])
1552 disp(['RF power absorbed (P/2piRp): ',num2str(data_proc.wave.wp_rf_2piRp(iw)),' MW/m ---> ',num2str(100*data_proc.wave.wp_rf_2piRp(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1553 for in_rf = 1:nn_rf,
1554 disp(['RF power absorbed on harm. ',num2str(data_proc.wave.wn_rf_list(in_rf)),' (P/2piRp): ',num2strdata_proc.wave.(wnp_rf_2piRp(iw,in_rf)),' MW/m ---> ',num2str(100*data_proc.wave.wnp_rf_2piRp(iw,in_rf)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1555 end
1556 disp(['Linear RF power absorbed (P/2piRp): ',num2str(data_proc.wave.wp_rf_2piRp_lin(iw)),' MW/m ---> ',num2str(100*data_proc.wave.wp_rf_2piRp_lin(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1557 end
1558 end
1559
1560 else
1561 disp(['Launched RF power (P): ',num2str(data_proc.scalar.p0_rf_2piRp*(2*pi*data_proc.scalar.Rp)),' MW'])
1562 if sum(data_proc.radial.P_rf) ~= 0,
1563 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1564 if dke_out.waves{1}.rayparam.colldamp == 1,
1565 disp(['RF absorbed power (RLA + NRCA) (P): ',num2str(data_proc.scalar.p_rf_2piRp_coll_luke*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_coll_luke/data_proc.scalar.p0_rf_2piRp),' %.'])
1566 disp(['Linear RF absorbed power (RLA + NRCA) (P): ',num2str(data_proc.scalar.p_rf_2piRp_lin*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_lin/data_proc.scalar.p0_rf_2piRp),' %.']);
1567 else
1568 disp(['RF absorbed power (RLA) (P): ',num2str(data_proc.scalar.p_rf_2piRp*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp/data_proc.scalar.p0_rf_2piRp),' %.'])
1569 disp(['Linear RF absorbed power (RLA) (P): ',num2str(data_proc.scalar.p_rf_2piRp_lin*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_lin/data_proc.scalar.p0_rf_2piRp),' %.']);
1570 end
1571 else
1572 disp(['RF absorbed power (RLA) (P): ',num2str(data_proc.scalar.p_rf_2piRp*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp/data_proc.scalar.p0_rf_2piRp),' %.'])
1573 disp(['Linear RF absorbed power (RLA) (P): ',num2str(data_proc.scalar.p_rf_2piRp_lin*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.p_rf_2piRp_lin/data_proc.scalar.p0_rf_2piRp),' %.']);
1574 end
1575
1576 for in_rf = 1:nn_rf,
1577 disp(['RF power absorbed on harm. ',num2str(data_proc.scalar.n_rf_list(in_rf)),' (P): ',num2str(data_proc.scalar.np_rf_2piRp(in_rf)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.scalar.np_rf_2piRp(in_rf)/data_proc.scalar.p0_rf_2piRp),' %.'])
1578 end
1579
1580 end
1581 disp(['Ohmic power deposited (P): ',num2str(data_proc.scalar.p_e_2piRp*(2*pi*data_proc.scalar.Rp)),' MW'])
1582 disp(['Total power deposited (P): ',num2str(data_proc.scalar.p_tot_2piRp*(2*pi*data_proc.scalar.Rp)),' MW'])
1583 disp(['CD efficiency (I/P): ',num2str(data_proc.scalar.EtaRp/data_proc.scalar.Rp),' A/W'])
1584 disp(['Initial loop voltage (V) @ Bmin: ',num2str(data_proc.scalar.V_loop_init),' V'])
1585 disp(['Initial flux-surface averaged loop voltage (V): ',num2str(data_proc.scalar.V_loop_init_fsav),' V'])
1586
1587 if opt.peak == 'y',
1588
1589 for ipeak = 1:length(data_proc.peaks.ir),
1590 disp(' ');
1591 disp('----------------------------------------');
1592 disp(' ');
1593 disp(['RF Peak # ',num2str(ipeak)])
1594 disp(['Position rho G: ',num2str(data_proc.peaks.rhoG(ipeak))])
1595 disp(['Position rho P: ',num2str(data_proc.peaks.rhoP(ipeak))])
1596 disp(['Position rho T: ',num2str(data_proc.peaks.rhoT(ipeak))])
1597 disp(['Position rho V: ',num2str(data_proc.peaks.rhoV(ipeak))])
1598 disp(['Temperature Te: ',num2str(data_proc.peaks.Te(ipeak)),' keV'])
1599 disp(['Density ne: ',num2str(data_proc.peaks.ne(ipeak)/1e20),' *10^20 m^-3'])
1600 disp(' ');
1601 disp(['RF power (P): ',num2str(data_proc.peaks.p_rf_2piRp(ipeak)*(2*pi*data_proc.scalar.Rp)),' MW'])
1602 disp(['Total current (I): ',num2str(data_proc.peaks.I_tot(ipeak)),' MA'])
1603 disp(['CD efficiency (I/P): ',num2str(data_proc.peaks.EtaRp(ipeak)/data_proc.scalar.Rp),' A/W'])
1604 disp(['Normalized CD efficiency xi: ',num2str(data_proc.peaks.ksirf(ipeak))])
1605 end
1606 end
1607
1608 if opt.wave == 'y',
1609
1610 for iw = 1:nw,
1611 disp(' ');
1612 disp('----------------------------------------');
1613 disp(' ');
1614 disp(['RF Wave # ',num2str(iw),' : ',data_proc.wave.wids{iw}])
1615 disp(['Position rho G: ',num2str(data_proc.wave.wP_rf_rhoG(iw))])
1616 disp(['Position rho P: ',num2str(data_proc.wave.wP_rf_rhoP(iw))])
1617 disp(['Position rho T: ',num2str(data_proc.wave.wP_rf_rhoT(iw))])
1618 disp(['Position rho V: ',num2str(data_proc.wave.wP_rf_rhoV(iw))])
1619 disp(['Temperature Te: ',num2str(data_proc.wave.wTe(iw)),' keV'])
1620 disp(['Density ne: ',num2str(data_proc.wave.wne(iw)/1e20),' *10^20 m^-3'])
1621 disp(' ');
1622 disp(['Launched RF power (P): ',num2str(data_proc.wave.wp0_rf_2piRp(iw)*(2*pi*data_proc.scalar.Rp)),'MW'])
1623 disp(['RF power absorbed (P): ',num2str(data_proc.wave.wp_rf_2piRp(iw)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.wave.wp_rf_2piRp(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1624 disp(['RF power absorbed with NRCA (P): ',num2str(data_proc.wave.wp_rf_2piRp_coll_luke(iw)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.wave.wp_rf_2piRp_coll_luke(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1625 for in_rf = 1:nn_rf,
1626 disp(['RF power absorbed on harm. ',num2str(data_proc.wave.wn_rf_list(in_rf)),' (P): ',num2str(data_proc.wave.wnp_rf_2piRp(iw,in_rf)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.wave.wnp_rf_2piRp(iw,in_rf)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1627 end
1628 disp(['Linear RF power absorbed (P): ',num2str(data_proc.wave.wp_rf_2piRp_lin(iw)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.wave.wp_rf_2piRp_lin(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1629 disp(['Linear RF power absorbed with NRCA (P): ',num2str(data_proc.wave.wp_rf_2piRp_lin_coll(iw)*(2*pi*data_proc.scalar.Rp)),' MW ---> ',num2str(100*data_proc.wave.wp_rf_2piRp_lin_coll(iw)/data_proc.wave.wp0_rf_2piRp(iw)),' %.'])
1630 end
1631 end
1632 end
1633
1634 diary off
1635
1636
1637
1638 if ~isfield(opt,'axs'),
1639 red = 0.9;
1640 lspace = 0.7;
1641 lspace2 = 0.5;
1642 bspace = 0.7;
1643 bspace2 = 0.5;
1644
1645 else
1646 red = 1;
1647 lspace = NaN;
1648 lspace2 = NaN;
1649 bspace = NaN;
1650 bspace2 = NaN;
1651
1652 end
1653
1654 tit = ['Launched RF power (P): ',num2str(data_proc.scalar.p0_rf_2piRp*(2*pi*data_proc.scalar.Rp)),' MW'];
1655
1656 xlim_rho = [0,1];
1657 xtick_rho = 0:0.2:1;
1658
1659 if isfield(dke_out,'waves') && iscell(dke_out.waves) && ~isempty(dke_out.waves),
1660 if ~isfield(opt,'axs'),
1661
1662 figure(1),clf,set(1,'Name',['RF power density profile - ',data_proc.simul.id])
1663
1664 ylab = 'P_{RF} (MW/m^3)';
1665
1666 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1667 if dke_out.waves{1}.rayparam.colldamp == 1,
1668 if isnan(data_proc.radial.xp_rf_2piRp_coll_luke)
1669 graph1D_jd(xrho,data_proc.radial.P_rf,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,'-',marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1670 else
1671 graph1D_jd(xrho,data_proc.radial.P_rf,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,'--',marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1672 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,'-',marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1673 legend('LUKE RF absorbed power (RLA)','LUKE RF absorbed power (RLA + NRCA)')
1674 end
1675 else
1676 graph1D_jd(xrho,data_proc.radial.P_rf,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,'-',marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1677 end
1678 else
1679 graph1D_jd(xrho,data_proc.radial.P_rf,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,'-',marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1680 end
1681
1682 set(gca,'xtick',xtick_rho);
1683
1684 print_jd(p_opt,[data_proc.simul.id,'_Prf_profile'],'figures_LUKE');
1685
1686 end
1687
1688
1689
1690 if ~isfield(opt,'axs'),
1691 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1692 if dke_out.waves{1}.rayparam.colldamp == 1,
1693 figure(101),clf,set(101,'Name',['RF power deposition profile (RLA + NRCA, LUKE and linear) - ',data_proc.simul.id]);
1694 else
1695 figure(101),clf,set(101,'Name',['RF power deposition profile (RLA, LUKE and linear) - ',data_proc.simul.id]);
1696 end
1697 else
1698 figure(101),clf,set(101,'Name',['RF power deposition profile (RLA, LUKE and linear) - ',data_proc.simul.id]);
1699 end
1700 ax = gca;
1701 else
1702 ax = opt.axs(1);
1703 end
1704
1705 ylab = 'P_{RF}/d\rho/(2\piR_p) (MW/m)';
1706
1707 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1708 if dke_out.waves{1}.rayparam.colldamp == 1,
1709 leg = {'C3PO (RLA + NRCA)','LUKE (RLA + NRCA)','C3PO (RLA)','LUKE (RLA)'};
1710 else
1711 leg = {'C3PO (RLA)','LUKE (RLA)'};
1712 end
1713 else
1714 leg = {'C3PO (RLA)','LUKE (RLA)'};
1715 end
1716
1717 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1718 if dke_out.waves{1}.rayparam.colldamp == 1,
1719 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_lin_coll./dxrho,0,0,'','','',NaN,xlim_rho,NaN,'-','+','b',width2,opt.font,ax);
1720 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_coll_luke./dxrho,0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,'-','o','r',width2,opt.font,ax);
1721 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_lin./dxrho,0,0,'','','',NaN,xlim_rho,NaN,'--','+','b',width2,opt.font,ax);
1722 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp./dxrho,0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,'--','o','r',width2,opt.font,ax,red,lspace,bspace);
1723 else
1724 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_lin./dxrho,0,0,'','','',NaN,xlim_rho,NaN,'--','+','b',width2,opt.font,ax);
1725 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp./dxrho,0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,'-','o','r',width2,opt.font,ax,red,lspace,bspace);
1726 end
1727 else
1728 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp_lin./dxrho,0,0,'','','',NaN,xlim_rho,NaN,'--','+','b',width2,opt.font,ax);
1729 graph1D_jd(xrho,data_proc.radial.xp_rf_2piRp./dxrho,0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,'-','o','r',width2,opt.font,ax,red,lspace,bspace);
1730 end
1731
1732 set(ax,'xtick',xtick_rho);
1733
1734 print_jd(p_opt,[data_proc.simul.id,'_Prf_lincomp_profile'],'figures_LUKE');
1735
1736
1737
1738 if ~isfield(opt,'axs') & opt.wave == 'y' & nw > 1,
1739
1740 figure(102),clf,set(102,'Name',['RF power deposition profile per wave # - ',data_proc.simul.id])
1741
1742 leg = cell(1,2*nw);
1743
1744 for iw = 1:nw,
1745
1746 if ~isempty(data_proc.wave.wids{iw}),
1747 leg{nw+iw} = data_proc.wave.wids{iw};
1748 else
1749 leg{nw+iw} = ['# ',num2str(iw)];
1750 end
1751 leg{iw} = [leg{nw+iw},' (linear)'];
1752 end
1753
1754 for iw = 1:nw,
1755 graph1D_jd(xrho,data_proc.wave.wxp_rf_2piRp_lin(iw,:)./dxrho,0,0,'','','',NaN,xlim_rho,NaN,style2,marker,colors{iw},width,opt.font);
1756 end
1757
1758 for iw = 1:nw-1,
1759 graph1D_jd(xrho,data_proc.wave.wxp_rf_2piRp(iw,:)./dxrho,0,0,'','','',NaN,xlim_rho,NaN,style,marker1,colors{iw},width2,opt.font);
1760 end
1761 graph1D_jd(xrho,data_proc.wave.wxp_rf_2piRp(nw,:)./dxrho,0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,style,marker1,colors{nw},width2,opt.font,gca,red,lspace,bspace);
1762
1763 set(gca,'xtick',xtick_rho);
1764
1765 print_jd(p_opt,[data_proc.simul.id,'_Prf_wave_profile'],'figures_LUKE');
1766
1767 end
1768
1769
1770
1771 if ~isfield(opt,'axs') && nn_rf > 1,
1772
1773 figure(103),clf,set(103,'Name',['RF power deposition profile per harmonic # - ',data_proc.simul.id])
1774
1775 leg = cell(1,nn_rf);
1776
1777 for in_rf = 1:nn_rf,
1778 leg{in_rf} = ['n = ',num2str(data_proc.scalar.n_rf_list(in_rf))];
1779 end
1780
1781 graph1D_jd(xrho,data_proc.radial.xnp_rf_2piRp./repmat(dxrho.',[1,nn_rf]),0,0,xlab_rho,ylab,tit,leg,xlim_rho,NaN,style,marker1,NaN,width2,opt.font,gca,red,lspace,bspace);
1782
1783 set(gca,'xtick',xtick_rho);
1784
1785 print_jd(p_opt,[data_proc.simul.id,'_Prf_nrf_profile'],'figures_LUKE');
1786
1787 end
1788 end
1789
1790
1791
1792 if ~isfield(opt,'axs'),
1793 figure(2),clf,set(2,'Name',['Ohmic power density profile - ',data_proc.simul.id])
1794 ax = gca;
1795 else
1796 ax = opt.axs(2);
1797 end
1798
1799 ylab = 'P_{OHM} (MW/m^3)';
1800
1801 graph1D_jd(xrho,data_proc.radial.P_e,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,style,marker1,color2,width2,opt.font,ax,red,lspace,bspace);
1802
1803 set(ax,'xtick',xtick_rho);
1804
1805 print_jd(p_opt,[data_proc.simul.id,'_Pe_profile'],'figures_LUKE');
1806
1807
1808
1809 if ~isfield(opt,'axs'),
1810
1811 figure(3),clf,set(3,'Name',['Total power density profile - ',data_proc.simul.id])
1812
1813 ylab = 'P_{TOT} (MW/m^3)';
1814
1815 graph1D_jd(xrho,data_proc.radial.P_tot,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1816
1817 set(gca,'xtick',xtick_rho);
1818
1819 print_jd(p_opt,[data_proc.simul.id,'_Ptot_profile'],'figures_LUKE');
1820
1821 end
1822
1823
1824
1825 if ~isfield(opt,'axs'),
1826 figure(4),clf,set(4,'Name',['Co-current density profile - ',data_proc.simul.id])
1827 ax = gca;
1828 else
1829 ax = opt.axs(3);
1830 end
1831
1832 ylab = 'J (MA/m^2)';
1833
1834 graph1D_jd(xrho,data_proc.radial.J_tot,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,style,marker1,color2,width2,opt.font,ax,red,lspace,bspace);
1835
1836 set(ax,'xtick',xtick_rho);
1837 set(ax,'ytickmode','auto');drawnow;
1838
1839 print_jd(p_opt,[data_proc.simul.id,'_J_profile'],'figures_LUKE');
1840
1841
1842
1843 if ~isfield(opt,'axs'),
1844
1845 figure(401),clf,set(401,'Name',['Parallel Momentum transfer rate density profile - ',data_proc.simul.id])
1846
1847 ylab = 'dp_{RF}/dt (N/m^3)';
1848
1849 graph1D_jd(xrho,data_proc.radial.xdpdt_rf,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,NaN,style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1850
1851 set(gca,'xtick',xtick_rho);
1852
1853 print_jd(p_opt,[data_proc.simul.id,'_Ptrans_profile'],'figures_LUKE');
1854
1855 end
1856
1857
1858
1859 if ~isfield(opt,'axs') && isfield(data_proc,'timevol') && isfield(data_proc.timevol,'I_tot_rf'),
1860
1861 figure(402),clf,set(402,'Name',['Current convergence - ',data_proc.simul.id])
1862
1863 xlab = 'RF iterations';
1864 ylab = 'I (MA)';
1865 xlim = [0,nit_rf + 1];
1866
1867 if ~isnan(I_exp),
1868 graph1D_jd(xlim,I_exp*[1,1],0,0,'','','',NaN,xlim,NaN,style,marker,color1,width2,opt.font);
1869 end
1870 graph1D_jd(1:nit_rf,data_proc.timevol.I_tot_rf,0,0,xlab,ylab,tit,NaN,xlim,0,style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1871
1872 print_jd(p_opt,[data_proc.simul.id,'_I_convergence'],'figures_LUKE');
1873 end
1874
1875
1876
1877 if ~isfield(opt,'axs') && ~all(isnan(data_proc.simul.residue_rf)),
1878
1879 figure(403),clf,set(403,'Name',['RF convergence - ',data_proc.simul.id])
1880
1881 xlab = 'RF iterations';
1882 ylab = 'Conv. Param.';
1883 xlim = [0,nit_rf + 1];
1884
1885 graph1D_jd(1:nit_rf,data_proc.simul.residue_rf,0,1,xlab,ylab,tit,NaN,xlim,NaN,style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1886
1887 print_jd(p_opt,[data_proc.simul.id,'_RF_convergence'],'figures_LUKE');
1888 end
1889
1890
1891
1892 if ~isfield(opt,'axs'),
1893
1894 figure(5),clf,set(5,'Name',['Current efficiency profile - ',data_proc.simul.id])
1895
1896 ylab = '\eta (A.m/W)';
1897
1898 if ~isnan(data_proc.scalar.P_rf_rhoG),
1899 eta_peak = interp1(data_proc.radial.xrhoG,data_proc.radial.eta,data_proc.scalar.P_rf_rhoG);
1900 elseif ~isnan(data_proc.scalar.jrhoG),
1901 eta_peak = interp1(data_proc.radial.xrhoG,data_proc.radial.eta,data_proc.scalar.jrhoG);
1902 else
1903 eta_peak = max(abs(data_proc.radial.eta));
1904 end
1905 ylim = 3*[min([0,eta_peak]),max([0,eta_peak])];
1906
1907 graph1D_jd(xrho,data_proc.radial.eta,0,0,xlab_rho,ylab,tit,NaN,xlim_rho,ylim,style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);
1908 set(gca,'xtick',xtick_rho);
1909
1910 print_jd(p_opt,[data_proc.simul.id,'_eta_profile'],'figures_LUKE');
1911
1912 end
1913
1914
1915
1916 if ~isfield(opt,'axs')
1917 if sum(data_proc.radial.P_rf) ~= 0,
1918
1919 figure(6),clf,set(6,'Name',['Current and RF power density profiles - ',data_proc.simul.id])
1920
1921 ylab1 = 'J (MA/m^2)';
1922 ylab2 = 'P_{RF} (MW/m^3)';
1923 ylim1 = NaN;
1924 ylim2 = NaN;
1925
1926 if isfield(dke_out.waves{1}.rayparam,'colldamp')
1927 if dke_out.waves{1}.rayparam.colldamp == 1,
1928 if isnan(data_proc.radial.xp_rf_2piRp_coll_luke)
1929 graph1D2jd(xrho,data_proc.radial.J_tot,xrho,data_proc.radial.P_rf,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,style,marker1,color2,width2,style,marker1,color3,width2,opt.font,red,lspace2,bspace);
1930 else
1931 graph1D2jd(xrho,data_proc.radial.J_tot,xrho,data_proc.radial.xp_rf_2piRp_coll_luke./equilDKE.xdV_2piRp_dke,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,style,marker1,color2,width2,style,marker1,color3,width2,opt.font,red,lspace2,bspace);
1932 end
1933 else
1934 graph1D2jd(xrho,data_proc.radial.J_tot,xrho,data_proc.radial.P_rf,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,style,marker1,color2,width2,style,marker1,color3,width2,opt.font,red,lspace2,bspace);
1935 end
1936 else
1937 graph1D2jd(xrho,data_proc.radial.J_tot,xrho,data_proc.radial.P_rf,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,style,marker1,color2,width2,style,marker1,color3,width2,opt.font,red,lspace2,bspace);
1938 end
1939
1940 print_jd(p_opt,[data_proc.simul.id,'_J_Prf_profile'],'figures_LUKE');
1941
1942 else
1943
1944 figure(6),clf,set(6,'Name',['Current and Ohmic power density profiles - ',data_proc.simul.id])
1945
1946 ylab1 = 'J (MA/m^2)';
1947 ylab2 = 'P_{Ohm} (MW/m^3)';
1948 ylim1 = NaN;
1949 ylim2 = NaN;
1950
1951 graph1D2jd(xrho,data_proc.radial.J_tot,xrho,data_proc.radial.P_tot,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,...
1952 style,marker1,color2,width2,style,marker1,color3,width2,opt.font,red,lspace2,bspace);
1953
1954 print_jd(p_opt,[data_proc.simul.id,'_J_Pohm_profile'],'figures_LUKE');
1955
1956 end
1957 end
1958
1959
1960
1961 if ~isfield(opt,'axs') && exist('Zmripple','var'),
1962 if ~isempty(Zmripple) && sum(sum(Zmripple.Xnhuloss)) ~= 0,
1963 figure(7),clf,set(7,'Name',['Magnetic ripple losses (rho) - ',data_proc.simul.id])
1964 [ax] = graph1D_jd(xrho,dke_out.xMRR_flux,0,0,'\rho','Ripple losses (m-3.s-1)',subtitle,NaN,[0,1],[0,max(max([dke_out.xMRR_flux(:),dke_out.xMRR_tau(:)]))],style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);hold on
1965 [ax] = graph1D_jd(xrho,dke_out.xMRR_tau,0,0,'\rho','Ripple losses (m-3.s-1)',subtitle,NaN,[0,1],[0,max(max([dke_out.xMRR_flux(:),dke_out.xMRR_tau(:)]))],style,marker1,color3,width2,opt.font,gca,red,lspace,bspace);hold on
1966 legend('Boundary losses','Volumic losses');
1967 print_jd(p_opt,[data_proc.simul.id,'_RippleLosses_profile'],'figures_LUKE');
1968
1969 figure(70),clf,set(70,'Name',['Magnetic ripple losses (R) - ',data_proc.simul.id])
1970 [ax] = graph1D_jd(equilDKE.Rp+equilDKE.Xx(:,1),dke_out.xMRR_flux,0,0,'R(m)','Ripple losses (m-3.s-1)',subtitle,NaN,[equilDKE.Rp,equilDKE.Rp+max(max(equilDKE.Xx))],[0,max(max([dke_out.xMRR_flux(:),dke_out.xMRR_tau(:)]))],style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);hold on
1971 [ax] = graph1D_jd(equilDKE.Rp+equilDKE.Xx(:,1),dke_out.xMRR_tau,0,0,'R(m)','Ripple losses (m-3.s-1)',subtitle,NaN,[equilDKE.Rp,equilDKE.Rp+max(max(equilDKE.Xx))],[0,max(max([dke_out.xMRR_flux(:),dke_out.xMRR_tau(:)]))],style,marker1,color3,width2,opt.font,gca,red,lspace,bspace);hold on
1972 legend('Boundary losses','Volumic losses');
1973 print_jd(p_opt,[data_proc.simul.id,'_RippleLosses_profile'],'figures_LUKE');
1974
1975 figure(8),clf,set(8,'Name',['Cones angles - ',data_proc.simul.id])
1976 [ax] = graph1D_jd(xrho,Zmripple.xthetabanana,0,0,'\rho','Angle (Deg.)',subtitle,NaN,[0,1],[0,60],style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);hold on
1977 [ax] = graph1D_jd(xrho,Zmripple.xthetaloss,0,0,'\rho','Angle (Deg.)',subtitle,NaN,[0,1],[0,60],style,marker1,color3,width2,opt.font,gca,red,lspace,bspace);hold on
1978 legend('Banana cone angle','Loss cone angle');
1979 print_jd(p_opt,[data_proc.simul.id,'_RippleConeAngles_profile'],'figures_LUKE');
1980
1981 ylab1 = 'p/p_{Te}';
1982 ylab2 = 'p/p_{Te(loc.)}';
1983 ytick = [0:2:10];
1984 ylim1 = [0,10];
1985 ylim2 = [0,40];
1986 figure(9),clf,set(9,'Name',['Collision detrapping threshold - ',data_proc.simul.id])
1987 [ax] = graph1D_jd(xrho,Zmripple.xpndetrap,0,0,'\rho','',subtitle,NaN,[0,1],'',style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);hold on
1988 [ax] = graph1D_jd(xrho,Zmripple.xpnth,0,0,'\rho','',subtitle,NaN,[0,1],'',style,marker1,color2,width2,opt.font,gca,red,lspace,bspace);hold on
1989 [ax] = graph1D2jd(xrho,Zmripple.xpnth,xrho,Zmripple.xpndetrap./Zmripple.xpnth,xlab_rho,ylab1,ylab2,tit,xlim_rho,ylim1,ylim2,style,marker1,color3,width2,style,marker1,color4,width2,opt.font,red,lspace2,bspace);
1990 set(gca,'ytick',ytick);
1991 legend('p_{c}/p_{Te}','p_{Te(loc.)}/p_{Te}','p_{c}/p_{Te(loc.)}');
1992 print_jd(p_opt,[data_proc.simul.id,'_RippleDetrapping_profile'],'figures_LUKE');
1993 end
1994 end
1995
1996 if ~isfield(opt,'axs') && isfield(data_proc,'timevol') && isfield(data_proc.timevol,'I_tot_rf') && isfield(opt,'its_rf') && any(opt.its_rf),
1997
1998 if islogical(opt.its_rf),
1999 its_rf = input_dke_yp('RF interation indices to compare',[1,nit_rf],[1;nit_rf]);
2000 else
2001 its_rf = opt.its_rf;
2002 end
2003
2004 nits_rf = length(its_rf);
2005
2006
2007
2008 figure(404),clf,set(404,'Name',['Current density profile it_rf comp - ',data_proc.simul.id])
2009
2010 ylab = 'J (MA/m^2)';
2011 leg = cell(1,nits_rf);
2012
2013 for iit_rf = 1:nits_rf-1,
2014 graph1D_jd(xrho,Zcurr(its_rf(iit_rf),end).x_0_fsav*mksa.j_ref*scocurr,0,0,'','','',NaN,xlim_rho,NaN,style,marker1,colors{iit_rf},width2,opt.font);
2015 leg{iit_rf} = ['Iteration # ',num2str(its_rf(iit_rf))];
2016 end
2017 leg{nits_rf} = ['Iteration # ',num2str(its_rf(nits_rf))];
2018 graph1D_jd(xrho,Zcurr(its_rf(nits_rf),end).x_0_fsav*mksa.j_ref*scocurr,0,0,xlab,ylab,tit,leg,xlim_rho,NaN,style,marker1,colors{nits_rf},width2,opt.font,gca,red,lspace,bspace);
2019
2020 set(ax,'xtick',xtick_rho);
2021 set(ax,'ytickmode','auto')
2022
2023 print_jd(p_opt,[data_proc.simul.id,'_J_rfcomp_profile'],'figures_LUKE');
2024
2025
2026
2027 figure(405),clf,set(404,'Name',['Power density profile it_rf comp - ',data_proc.simul.id])
2028
2029 ylab = 'P_{RF} (MW/m^3)';
2030
2031 for iit_rf = 1:nits_rf-1,
2032 graph1D_jd(xrho,ZP0(its_rf(iit_rf),end).x_rf_fsav.'*mksa.P_ref,0,0,'','','',NaN,xlim_rho,NaN,style,marker1,colors{iit_rf},width2,opt.font);
2033 end
2034 graph1D_jd(xrho,ZP0(its_rf(nits_rf),end).x_rf_fsav.'*mksa.P_ref,0,0,xlab,ylab,tit,leg,xlim_rho,NaN,style,marker1,colors{nits_rf},width2,opt.font,gca,red,lspace,bspace);
2035
2036 set(ax,'xtick',xtick_rho);
2037 set(ax,'ytickmode','auto')
2038
2039 print_jd(p_opt,[data_proc.simul.id,'_P_rfcomp_profile'],'figures_LUKE');
2040
2041 end
2042
2043 if opt.spec == 1 && sum(data_proc.radial.P_rf) ~= 0,
2044 disp('******** power absorption along specific rays ********');
2045
2046 iw = NaN;
2047 ibb = NaN;
2048 ib = 0;
2049 iwdef = 0;
2050
2051 while iw ~= 0 && ibb ~= 0
2052
2053 ib = ib + 1;
2054
2055 iw = input_dke_yp('wave index (0 to return)',iwdef,0:nw,'',[1,1]);
2056
2057 if iw == 0,
2058 continue
2059 end
2060
2061 ibb = input_dke_yp('ray index (0 to return)',1,0:nb(iw),'',[1,1]);
2062
2063 if ibb == 0,
2064 continue
2065 end
2066
2067 iwdef = 0;
2068
2069 iy = sum(nb(1:iw - 1)) + ibb;
2070
2071 ray = waves{iw}.rays{ibb};
2072
2073 sS_c3po = ray.ss;
2074 if isfield(waves{iw}.rayparam,'colldamp') && waves{iw}.rayparam.colldamp == 1,
2075 sP_c3po = ray.sP_2piRp_lin_coll*2*pi*data_proc.scalar.Rp;
2076 else
2077 sP_c3po = ray.sP_2piRp_lin*2*pi*data_proc.scalar.Rp;
2078 end
2079
2080 sS_luke = zS{iy};
2081 if isfield(waves{iw}.rayparam,'colldamp') && waves{iw}.rayparam.colldamp == 1,
2082 sP_luke = zP_0_2piRp_mod_coll{end}{iy}*2*pi*data_proc.scalar.Rp;
2083 else
2084 sP_luke = zP_0_2piRp_mod{end}{iy}*2*pi*data_proc.scalar.Rp;
2085 end
2086
2087 if isfield(waves{iw}.rayparam,'colldamp') && waves{iw}.rayparam.colldamp == 1,
2088 leg = {'Ray-Tracing (NRCA + RLA)','Fokker-Planck (NRCA + RLA)'};
2089 else
2090 leg = {'Ray-Tracing (RLA)','Fokker-Planck (RLA)'};
2091 end
2092
2093
2094
2095 figlab = ['wave #',num2str(iw),'; ray #',num2str(ibb)];
2096
2097 figure(14+iy),clf;set(gcf,'name',figlab);
2098
2099 xlim = [0,max([sS_c3po,sS_luke])];
2100 ylim = [0,1.2*max([sP_c3po,sP_luke])/1e3];
2101 xlab = 's (m)';
2102 ylab = 'P (kW)';
2103
2104 graph1D_jd(sS_c3po,sP_c3po/1e3,0,0,xlab,ylab,figlab,NaN,xlim,ylim,style2,marker,color3,width2,opt.font,gca,red,lspace,bspace2);
2105 graph1D_jd(sS_luke,sP_luke/1e3,0,0,'','','',leg,xlim,ylim,style,marker2,color2,width2,opt.font);
2106
2107 if timevol == 2 && isfield(opt,'its_rf') && any(opt.its_rf) && length(zP_0_2piRp_mod) == nit_rf;
2108
2109 figure(414+iy),clf,set(404,'Name',['Power it_rf comp - ',figlab])
2110
2111 for iit_rf = 1:nits_rf-1,
2112 graph1D_jd(sS_luke,zP_0_2piRp_mod{its_rf(iit_rf)}{iy}*2*pi*data_proc.scalar.Rp/1e3,0,0,'','','',NaN,xlim,ylim,style,marker1,colors{iit_rf},width2,opt.font);
2113 end
2114 graph1D_jd(sS_luke,zP_0_2piRp_mod{its_rf(nits_rf)}{iy}*2*pi*data_proc.scalar.Rp/1e3,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker1,colors{nits_rf},width2,opt.font,gca,red,lspace,bspace);
2115
2116 set(ax,'ytickmode','auto')
2117
2118 end
2119
2120 end
2121
2122 end
2123
2124 else
2125
2126 disp(['normalized current: ',num2str(Zcurr(end).x_0_fsav)])
2127 disp(['normalized RF power: ',num2str(ZP0(end).x_rf_fsav)])
2128 disp(['normalized Ohmic power: ',num2str(ZP0(end).x_e_fsav)])
2129 disp(['normalized Total power: ',num2str(ZP0(end).x_e_fsav+ZP0(end).x_rf_fsav)])
2130 disp(['normalized CD efficiency: ',num2str(Zcurr(end).x_0_fsav./(ZP0(end).x_e_fsav+ZP0(end).x_rf_fsav))])
2131
2132 diary off
2133
2134 end
2135
2136 end
2137
2138 if isnan(ir_display),
2139
2140 if isfield(opt,'axs'),
2141 recset_jd(opt.axs(4:5),'Visible','off',true)
2142 end
2143
2144 sir = cell(1,nr+1);
2145 for ir = 0:nr,
2146 sir{ir+1} = num2str(ir);
2147 end
2148 ir_display = iselect_jd(sir,'Select radial position for display (0 to return)',opt.gui,opt.style.distselectstyle,1,'popupmenu') - 1;
2149
2150 if isfield(opt,'axs'),
2151 recset_jd(opt.axs(4:5),'Visible','on',true)
2152 end
2153
2154 end
2155
2156 if ir_display > 0 && ir_display <= nr,
2157
2158 tir_display = num2str(ir_display);
2159
2160 disp(['normalized radius for distribution plots: ',num2str(xrho(ir_display))])
2161
2162 if isfield(data,'savefile'),
2163
2164 if ~exist(data.savefile,'file'),
2165 disp(['WARNING : LUKE simulation file ',data.savefile,' cannot be found']);
2166 ir_display = 0;
2167 return
2168 end
2169
2170 lukedata = load(data.savefile);
2171
2172 dke_out = lukedata.lukestructs.output.dke_out;
2173 equilDKE = lukedata.lukestructs.output.equilDKE;
2174 momentumDKE = lukedata.lukestructs.output.momentumDKE;
2175 gridDKE = lukedata.lukestructs.output.gridDKE;
2176 Zmomcoef = lukedata.lukestructs.output.Zmomcoef;
2177 Zbouncecoef = lukedata.lukestructs.output.Zbouncecoef;
2178 Zmripple = lukedata.lukestructs.output.Zmripple;
2179 mksa = lukedata.lukestructs.output.mksa;
2180
2181 scocurr = sign(equilDKE.psia_apRp)*sign(equilDKE.XBphi(1,1));
2182 elseif isfield(data,'data_proc'),
2183 dke_out = data.dke_out;
2184 equilDKE = data.equilDKE;
2185 momentumDKE = data.momentumDKE;
2186 gridDKE = data.gridDKE;
2187 Zmomcoef = data.Zmomcoef;
2188 Zbouncecoef = data.Zbouncecoef;
2189 Zmripple = data.Zmripple;
2190 mksa = data.mksa;
2191
2192 scocurr = sign(equilDKE.psia_apRp)*sign(equilDKE.XBphi(1,1));
2193 end
2194
2195 pn = momentumDKE.pn;
2196 dpn = momentumDKE.dpn;
2197 mhu = momentumDKE.mhu;
2198 dmhu = momentumDKE.dmhu;
2199 gamma = Zmomcoef.gamma;
2200 Ec_ref = Zmomcoef.Ec_ref;
2201 if isfield(Zmomcoef,'dEc_ref'),
2202 dEc_ref = Zmomcoef.dEc_ref;
2203 else
2204 [dummy,dummy,dummy,dummy,dummy,dummy,dummy,mc2] = pc_dke_yp;
2205 dEc_ref = diff(mc2*(sqrt(1 + momentumDKE.pn_S.*momentumDKE.pn_S*mksa.betath_ref^2) - 1));
2206 end
2207
2208 Xf0 = dke_out.XXf0{end}(:,:,ir_display);
2209 Xfinit = dke_out.XXfinit(:,:,ir_display);
2210 Xmhu = repmat(mhu,[length(pn),1]);
2211
2212 Xlambda = repmat(Zbouncecoef.Xxlambda(:,ir_display).',[length(pn),1]);
2213 qtilde = Zbouncecoef.xqtilde(ir_display);
2214 qhat = Zbouncecoef.xqhat(ir_display);
2215
2216
2217
2218 pdj0dp = scocurr*mksa.j_ref*2*pi*pn.*pn.*pn.*integral_dke_jd(dmhu,Xf0.*Xmhu,2).'./gamma;
2219
2220
2221
2222 pdN0_dEc_ref = 2*pi.*(qtilde./qhat).*integral_dke_jd(dmhu,Xlambda.*Xf0,2).'.*pn.*pn.*dpn./dEc_ref;
2223 pdNinit_dEc_ref = 2*pi.*(qtilde./qhat).*integral_dke_jd(dmhu,Xlambda.*Xfinit,2).'.*pn.*pn.*dpn./dEc_ref;
2224
2225
2226
2227 [logXf0_cyl_ref,ppar_cyl_ref,dppar_cyl_ref,pperp_cyl_ref,dpperp_cyl_ref] = s2c_dke_yp(log(Xf0),pn,mhu,dp_cyl);
2228 Xf0_cyl_ref = exp(logXf0_cyl_ref);
2229
2230 logXfinit_cyl_ref = s2c_dke_yp(log(Xfinit),pn,mhu,dp_cyl);
2231 Xfinit_cyl_ref = exp(logXfinit_cyl_ref);
2232
2233 XDparpar_cyl_ref = s2c_dke_yp(dke_out.ZXXD_rf.parpar_ij(:,:,ir_display),pn,mhu,dp_cyl);
2234 XDpp_cyl_ref = s2c_dke_yp(dke_out.ZXXD_rf.pp_ij(:,:,ir_display),pn,mhu,dp_cyl);
2235
2236 Xppar_cyl_ref = ones(length(pperp_cyl_ref),1)*ppar_cyl_ref(:)';
2237 Xpperp_cyl_ref = pperp_cyl_ref(:)*ones(1,length(ppar_cyl_ref));
2238 Xp_cyl_ref = sqrt(Xppar_cyl_ref.*Xppar_cyl_ref + Xpperp_cyl_ref.*Xpperp_cyl_ref);
2239
2240 Xf0_cyl_ref = Xf0_cyl_ref.*(Xf0_cyl_ref>0);
2241 Xf0_cyl_ref = Xf0_cyl_ref.*(Xp_cyl_ref<=max(ppar_cyl_ref));
2242 Xf0_cyl_ref(isnan(Xf0_cyl_ref)) = 0;
2243
2244 Xfinit_cyl_ref = Xfinit_cyl_ref.*(Xfinit_cyl_ref>0);
2245 Xfinit_cyl_ref = Xfinit_cyl_ref.*(Xp_cyl_ref<=max(ppar_cyl_ref));
2246 Xfinit_cyl_ref(isnan(Xfinit_cyl_ref)) = 0;
2247
2248
2249
2250 Xmask = abs(Xf0_cyl_ref) < eps;
2251
2252 pnmax_S = max(momentumDKE.pn_S);
2253 pnmax_0 = min(min(Xp_cyl_ref(Xmask)));
2254
2255 if isempty(pnmax),
2256 pnmax = iselect_jd({},'Select the maximum momentum for display',opt.gui,opt.style.editselectstyle,min([pnmax_0,pnmax_S]),'edit');
2257 end
2258
2259 if isnan(pnmax),
2260 pnmax = min([pnmax_0,pnmax_S]);
2261 end
2262
2263 Xmask = Xp_cyl_ref >= pnmax;
2264
2265 Xf0_cyl_ref(Xmask) = 0;
2266 Xfinit_cyl_ref(Xmask) = 0;
2267 XDparpar_cyl_ref(Xmask) = 0;
2268 XDpp_cyl_ref(Xmask) = 0;
2269
2270 Fpar0_cyl_ref = 2*pi*integral_dke_jd(dpperp_cyl_ref,Xpperp_cyl_ref.*Xf0_cyl_ref,1);
2271 Fparinit_cyl_ref = 2*pi*integral_dke_jd(dpperp_cyl_ref,Xpperp_cyl_ref.*Xfinit_cyl_ref,1);
2272
2273 data_proc.cyl.ppar_ref = ppar_cyl_ref;
2274 data_proc.cyl.Fpar0_ref = Fpar0_cyl_ref;
2275 data_proc.cyl.Fparinit_ref = Fparinit_cyl_ref;
2276 data_proc.cyl.Dparpar_ref = XDparpar_cyl_ref(1,:);
2277
2278 nhnorm = mean((Xfinit(1,:) - Xf0(1,:))./Xfinit(1,:));
2279
2280 disp(['Fraction of suprathermal electrons : ',num2str(nhnorm)]);
2281
2282 Xfh_cyl_ref = (Xf0_cyl_ref - (1 - nhnorm)*Xfinit_cyl_ref)/nhnorm;
2283
2284
2285
2286 radialloc = ['\rho = ',num2str(xrho(ir_display))];
2287
2288
2289
2290 if ~isfield(opt,'axs'),
2291 figure(10),clf,set(10,'Name',['dj/dp - ',data_proc.simul.id])
2292
2293 xlab = 'p/p_{Te}';
2294 ylab = 'dj/dp_n (MA/m^2)';
2295 xlim = [0,pnmax];
2296
2297 graph1D_jd(pn,pdj0dp,0,0,xlab,ylab,[radialloc,' ; dj/dp positive for co-current'],NaN,xlim,NaN,style,marker,color2,width2,opt.font,gca,red,lspace,bspace2);
2298
2299 print_jd(p_opt,[data_proc.simul.id,'_djdp_',tir_display],'figures_LUKE');
2300 end
2301
2302
2303
2304 if ~isfield(opt,'axs'),
2305 figure(11),clf,set(11,'Name',['2-D contour of f - ',data_proc.simul.id])
2306 ax = gca;
2307 else
2308 ax = opt.axs(4);
2309 end
2310
2311 x1 = ppar_cyl_ref;
2312 xlab1 = 'p_{||}/p_{Te}';
2313 xlim1 = [-pnmax,pnmax];
2314
2315 x2 = pperp_cyl_ref;
2316 xlab2 = 'p_{\perp}/p_{Te}';
2317 xlim2 = [0,pnmax];
2318
2319 leg = {'f_M','f'};
2320
2321 y1 = Xfinit_cyl_ref;
2322 y2 = Xf0_cyl_ref;
2323
2324 if max(max(XDparpar_cyl_ref)) > max(max(XDpp_cyl_ref)),
2325 y3 = XDparpar_cyl_ref;
2326 y3lab = 'D||';
2327 else
2328 y3 = XDpp_cyl_ref;
2329 y3lab = 'Dpp';
2330 end
2331 y3max = max(max(y3));
2332
2333 y9 = sqrt(pnmax^2 - x1.^2);
2334
2335 graph2D_jd(x1,x2,y1','','',radialloc,xlim1,xlim2,0,mksa.betath_ref,[0:0.5:pnmax],0,style,color3,width,opt.font,1,max(max(y1)),ax);
2336 hold(ax,'on')
2337 graph2D_jd(x1,x2,y2','','',radialloc,xlim1,xlim2,0,mksa.betath_ref,[0:0.5:pnmax],0,style,color2,width2,opt.font,1,max(max(y1)),ax);
2338
2339 if y3max > 1e-5,
2340 graph2D_jd(x1,x2,y3','','',radialloc,xlim1,xlim2,0,NaN,y3max/100,0,style,color4,width2,opt.font,1,NaN,ax);
2341 end
2342 hold(ax,'off')
2343
2344 graph1D_jd(x1,y9,0,0,xlab1,xlab2,radialloc,NaN,xlim1,xlim2,style,marker,color1,width2,opt.font,ax,red,lspace,bspace2);
2345
2346 axis(ax,'equal');axis(ax,[xlim1,xlim2]);
2347 legend(ax,leg);
2348
2349 if isfield(data_proc,'scalar') && data_proc.scalar.Rp ~= Inf,
2350 bounce_ripple_jd(1,pnmax,gridDKE.xmhubounce2,ir_display,Zmripple,opt.font,ax);
2351 end
2352
2353 print_jd(p_opt,[data_proc.simul.id,'_f_2D',tir_display],'figures_LUKE');
2354
2355 if ~isfield(opt,'axs') && max(max(abs(y3))) > 0,
2356
2357 figure(12),clf,set(12,'Name',['2-D contour of ',y3lab,' - ',data_proc.simul.id])
2358
2359 cont = input_dke_yp('contours for DQL',10);
2360
2361 graph2D_jd(x1,x2,y3','','',radialloc,xlim1,xlim2,1,NaN,cont,0,style,NaN,width,opt.font,1);
2362 graph1D_jd(x1,y9,0,0,xlab1,xlab2,radialloc,NaN,xlim1,xlim2,style,marker,color1,width2,opt.font,gca,red,lspace,bspace2);
2363 axis('equal');axis([xlim1,xlim2])
2364
2365 if isfield(data_proc,'scalar') && data_proc.scalar.Rp ~= Inf,
2366 bounce_ripple_jd(1,pnmax,gridDKE.xmhubounce2,ir_display,Zmripple);
2367 end
2368 colorbar
2369
2370 print_jd(p_opt,[data_proc.simul.id,'_DQL',tir_display],'figures_LUKE');
2371
2372 if strcmp(y3lab,'D||'),
2373
2374 figure(121),clf,set(121,'Name',['1-D plot of ',y3lab,' - ',data_proc.simul.id])
2375
2376 y1 = y3(1,:);
2377 ylab = 'D_{||} (p_\perp = 0)';
2378
2379 graph1D_jd(x1,y1,0,0,xlab1,ylab,radialloc,NaN,xlim1,NaN,style,marker,color2,width2,opt.font,gca,red,lspace,bspace2);
2380
2381 print_jd(p_opt,[data_proc.simul.id,'_DQL_1D',tir_display],'figures_LUKE');
2382
2383 end
2384 end
2385
2386
2387
2388 if ~isfield(opt,'axs'),
2389 figure(111),clf,set(11,'Name',['2-D contour of fh - ',data_proc.simul.id])
2390
2391 x1 = ppar_cyl_ref;
2392 xlab1 = 'p_{||}/p_{Te}';
2393 xlim1 = [-pnmax,pnmax];
2394
2395 x2 = pperp_cyl_ref;
2396 xlab2 = 'p_{\perp}/p_{Te}';
2397 xlim2 = [0,pnmax];
2398
2399 y = Xfh_cyl_ref;
2400
2401 if max(max(XDparpar_cyl_ref)) > max(max(XDpp_cyl_ref)),
2402 y3 = XDparpar_cyl_ref;
2403 y3lab = 'D||';
2404 else
2405 y3 = XDpp_cyl_ref;
2406 y3lab = 'Dpp';
2407 end
2408 y3max = max(max(y3));
2409
2410 y9 = sqrt(pnmax^2 - x1.^2);
2411
2412 graph2D_jd(x1,x2,y','','',radialloc,xlim1,xlim2,0,mksa.betath_ref,[0:0.5:pnmax],0,style,color3,width,opt.font,1,max(max(y1)));
2413 hold on
2414
2415 if y3max > 1e-5,
2416 graph2D_jd(x1,x2,y3','','',radialloc,xlim1,xlim2,0,NaN,y3max/100,0,style,color4,width2,opt.font);
2417 end
2418 hold off
2419
2420 graph1D_jd(x1,y9,0,0,xlab1,xlab2,radialloc,NaN,xlim1,xlim2,style,marker,color1,width2,opt.font,gca,red,lspace,bspace2);
2421
2422 axis('equal');axis([xlim1,xlim2]);
2423 legend(leg);
2424
2425 if isfield(data_proc,'scalar') && data_proc.scalar.Rp ~= Inf,
2426 bounce_ripple_jd(1,pnmax,gridDKE.xmhubounce2,ir_display,Zmripple);
2427 end
2428
2429 print_jd(p_opt,[data_proc.simul.id,'_fh_2D',tir_display],'figures_LUKE');
2430
2431 end
2432
2433
2434
2435 if ~isfield(opt,'axs'),figure(13),clf,set(13,'Name',['F||-F||M - ',data_proc.simul.id])
2436
2437 y1 = Fpar0_cyl_ref - Fparinit_cyl_ref;
2438 ylab = 'F_{||} - F_{||M}';
2439
2440 graph1D_jd(x1,y1,0,0,xlab1,ylab,radialloc,NaN,xlim1,NaN,style,marker,color2,width2,opt.font,gca,red,lspace,bspace2);
2441
2442 print_jd(p_opt,[data_proc.simul.id,'_F0mFM_1D',tir_display],'figures_LUKE');
2443 end
2444
2445
2446
2447 if ~isfield(opt,'axs'),
2448 figure(14),clf,set(14,'Name',['F|| - ',data_proc.simul.id])
2449 ax = gca;
2450 else
2451 ax = opt.axs(5);
2452 end
2453
2454 y1 = Fparinit_cyl_ref;
2455 y2 = Fpar0_cyl_ref;
2456 ylab = 'F_{||}';
2457 ylim = [eps,1];
2458
2459 graph1D_jd(x1,y1,0,1,'','',radialloc,NaN,xlim1,ylim,style,marker,color3,width,opt.font,ax,red,lspace,bspace2);
2460 graph1D_jd(x1,y2,0,1,xlab1,ylab,radialloc,NaN,xlim1,ylim,style,marker,color2,width2,opt.font,ax);
2461 legend(ax,leg);
2462
2463 print_jd(p_opt,[data_proc.simul.id,'_F0pFM_1D_',tir_display],'figures_LUKE');
2464
2465
2466
2467 if ~isfield(opt,'axs'),
2468 figure(15),clf,set(14,'Name',['dNdEc - ',data_proc.simul.id])
2469
2470 x1 = Ec_ref;
2471 xlab1 = 'E_c (keV)';
2472 xlim1 = [0,max(x1)];
2473 y1 = pdNinit_dEc_ref;
2474 y2 = pdN0_dEc_ref;
2475 ylab = 'dn_e/dE_c/n_0 (keV^{-1})';
2476
2477 graph1D_jd(x1,y1,0,1,'','',radialloc,NaN,xlim1,NaN,style2,marker,color3,width2,opt.font,gca,red,lspace,bspace2);
2478 graph1D_jd(x1,y2,0,1,xlab1,ylab,radialloc,leg,xlim1,NaN,style,marker,color2,width2,opt.font);
2479 legend(leg);
2480
2481 print_jd(p_opt,[data_proc.simul.id,'_dNdEc_',tir_display],'figures_LUKE');
2482 end
2483
2484 end
2485 end
2486
2487 function [rho,drho] = calc_peak_jd(xrho,xdrho,xQ,xdV,opt_dens,opt_peak)
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504 nr = length(xrho);
2505
2506 if opt_dens == 0,
2507 xprof = xQ.*xdV./xdrho;
2508 xint = xQ.*xdV;
2509 else
2510 xprof = xQ;
2511 xint = xQ.*xdrho;
2512 end
2513
2514 profmax = max(abs(xprof));
2515 ipeak = find(abs(xprof) == profmax,1,'first');
2516 peaksign = sign(xprof(ipeak));
2517
2518 xmaskp = ipeak:nr;
2519 xmaskm = 1:ipeak;
2520
2521 if opt_peak == 0,
2522
2523 ixp = find(diff(xprof(xmaskp))*peaksign >= 0,1,'first');
2524 xmaskp = xmaskp(1:min([ixp,length(xmaskp)]));
2525
2526 ixm = 1 + find(diff(xprof(xmaskm))*peaksign <= 0,1,'last');
2527 xmaskm = xmaskm(max([1,ixm]):end);
2528
2529 rho = xrho(ipeak);
2530 drho = interp1([peaksign*xprof(xmaskp),eps],[xrho(xmaskp),1],profmax/2) - interp1([eps,peaksign*xprof(xmaskm)],[eps,xrho(xmaskm)],profmax/2);
2531 else
2532
2533 ixp = find(xint(xmaskp)*peaksign <= 0,1,'first') - 1;
2534 xmaskp = xmaskp(1:min([ixp,length(xmaskp)]));
2535
2536 ixm = find(xint(xmaskm)*peaksign <= 0,1,'last') + 1;
2537 xmaskm = xmaskm(max([1,ixm]):end);
2538
2539 xmask = unique([xmaskm,xmaskp]);
2540
2541 rho = sum(xint(xmask).*xrho(xmask))/sum(xint(xmask));
2542 drho = 2*sqrt(2*log(2))*sqrt(sum(xint(xmask).*(xrho(xmask) - rho).^2)/sum(xint(xmask)));
2543 end
2544 end
2545
2546
2547