0001 function proc_wave_luke_jd(filename,opt_disp,nr_dep)
0002
0003
0004
0005 if nargin < 3,
0006 nr_dep = 20;
0007 end
0008
0009 if nargin < 2 | isempty(opt_disp),
0010 opt_disp.equil = 0;
0011 opt_disp.si = 1;
0012 end
0013
0014 load([filename,'.mat'],'equil','waves','dke_out','rho_S','ZP0','Zcurr','mksa','equilDKE','radialDKE','momentumDKE')
0015
0016 pn = momentumDKE.pn;
0017 mhu = momentumDKE.mhu;
0018 xrho = equilDKE.xrho;
0019 rdke = radialDKE.r_dke;
0020 xdrho = diff(radialDKE.xrho_S_dke);
0021
0022 clear gridDKE radialDKE
0023
0024 xys = dke_out.xys;
0025 xyP0_2piRp_mod = dke_out.xyP0_2piRp_mod;
0026 xyprop_dke = dke_out.xyprop_dke;
0027 yb = dke_out.yb;
0028 waveparam = dke_out.waveparam;
0029
0030 ZXXD_rf = dke_out.ZXXD_rf;
0031 ZXXD_rf_tp = dke_out.ZXXD_rf_tp;
0032 ZXXF_rf_tp = dke_out.ZXXF_rf_tp;
0033
0034 XXfinit = dke_out.XXfinit;
0035 XXf0 = dke_out.XXf0;
0036
0037 clear dke_out
0038
0039 xmhu0T2 = equilDKE.xmhu0T2;
0040 xdV_2piRp_dke = equilDKE.xdV_2piRp_dke;
0041
0042 clear equilDKE
0043
0044 xyPabs_rf_dke_fsav_mksa = ZP0.xy_rf_fsav*mksa.P_ref*1e6;
0045 xJ_rf_dke_fsav_mksa = Zcurr.x_0_fsav*mksa.j_ref*1e6;
0046 betath_ref = mksa.betath_ref;
0047
0048 clear mksa ZP0 Zcurr
0049
0050 if opt_disp.equil == 1,
0051 equilibrium_jd(equil,NaN,2);
0052 end
0053
0054 [nr,ny] = size(xys);
0055
0056 xydV_2piRp_dke = xdV_2piRp_dke'*ones(1,ny);
0057 xydPrf_dke = xyPabs_rf_dke_fsav_mksa.*xydV_2piRp_dke*2*pi*equil.Rp;
0058
0059 xdPrf_dke = sum(xydPrf_dke,2).';
0060
0061 xdPdrho = xdPrf_dke./xdrho;
0062
0063 zs = raypath_jd(xys,xyprop_dke,yb,1:nr);
0064 zP0_2piRp_mod = raypath_jd(xyP0_2piRp_mod,xyprop_dke,yb,1:nr);
0065
0066 rdPdrho_lin = 0;
0067 P0tot_2piRp = 0;
0068 ibb = 0;
0069
0070 for iw = 1:length(waves),
0071
0072
0073
0074 for ib = 1:length(waves{iw}.rays),
0075
0076 ibb = ibb + 1;
0077
0078 ray = waves{iw}.rays{ib};
0079 P0tot_2piRp = P0tot_2piRp + ray.P0_2piRp;
0080
0081 ray.srho = psi2rho_jd(equil,ray.spsin);
0082 rdPdrho_lin = rdPdrho_lin + 2*pi*equil.Rp*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin,nr_dep);
0083
0084 ray.zs = zs{ibb};
0085 ray.zP0_2piRp_mod = zP0_2piRp_mod{ibb};
0086
0087 zs_fit = ray.zs;
0088 zP0_2piRp_mod_fit = ray.zP0_2piRp_mod;
0089
0090 irep = find(diff(zs_fit) == 0);
0091
0092 zs_fit(irep+1) = [];
0093 zP0_2piRp_mod_fit(irep) = (zP0_2piRp_mod_fit(irep) + zP0_2piRp_mod_fit(irep+1))/2;
0094 zP0_2piRp_mod_fit(irep+1) = [];
0095
0096 if zs_fit(1) > ray.ss(1),
0097 zs_fit = [ray.ss(1),zs_fit];
0098 zP0_2piRp_mod_fit = [ray.zP0_2piRp_mod(1),zP0_2piRp_mod_fit];
0099 end
0100
0101 if zs_fit(end) < ray.ss(end),
0102 zs_fit = [zs_fit,ray.ss(end)];
0103 zP0_2piRp_mod_fit = [zP0_2piRp_mod_fit,ray.zP0_2piRp_mod(end)];
0104 end
0105
0106 ray.sP_2piRp = interp1(zs_fit,zP0_2piRp_mod_fit,ray.ss,'pchip');
0107
0108 sP_2piRp_norm = ray.sP_2piRp/ray.P0_2piRp;
0109 delta = (1 - erf(1))/2;
0110
0111 ray.isabsmin = max(find(sP_2piRp_norm > (1 - delta)));
0112 ray.isabsmax = max(find(sP_2piRp_norm > delta));
0113 ray.ismax = max(find(diff(ray.sP_2piRp) == min(diff(ray.sP_2piRp))));
0114
0115 Pin.dke(iw,ib) = ray.zP0_2piRp_mod(1)*2*pi*equil.Rp;
0116 Pout.dke(iw,ib) = ray.zP0_2piRp_mod(end)*2*pi*equil.Rp;
0117
0118 Pin.alpha(iw,ib) = ray.sP_2piRp_lin(1)*2*pi*equil.Rp;
0119 Pout.alpha(iw,ib) = ray.sP_2piRp_lin(end)*2*pi*equil.Rp;
0120
0121 ray.absrate_dke = 1 - Pout.dke(iw,ib)/Pin.dke(iw,ib);
0122 ray.absrate_alpha = 1 - Pout.alpha(iw,ib)/Pin.alpha(iw,ib);
0123
0124 waves{iw}.rays{ib} = ray;
0125
0126 end
0127
0128 absrate.dke(iw) = 1 - sum(Pout.dke(iw,:))./sum(Pin.dke(iw,:));
0129 absrate.alpha(iw) = 1 - sum(Pout.alpha(iw,:))./sum(Pin.alpha(iw,:));
0130
0131 end
0132
0133 rho_dep_S = linspace(0,1,nr_dep+1);
0134 rho_dep = (rho_dep_S(1:end-1) + rho_dep_S(2:end))/2;
0135
0136 close all
0137
0138 color_list = {[0,0.75,0],'m'};
0139
0140 xlab = '';
0141 ylab = '';
0142 tit = '';
0143 leg = NaN;
0144 xlim = get(gca,'xlim');
0145 ylim = get(gca,'ylim');
0146 style = '-';
0147 style2 = '--';
0148 width = 0.5;
0149 width2 = 2;
0150 marker = 'none';
0151 marker2 = '+';
0152 siz = 20;
0153 red = 0.9;
0154 lspace = 0.7;
0155 bspace = 0.7;
0156
0157 color0 = NaN;
0158 color1 = 'r';
0159 color2 = 'b';
0160 color3 = 'g';
0161 color5 = 'k';
0162
0163 figure(3);clf
0164
0165 xlab = 'r/a';
0166 ylab = 'dP/d\rho (norm.)';
0167 xlim = NaN;
0168 ylim = NaN;
0169 leg = ['LUKE ';...
0170 'linear'];
0171
0172 graph1D_jd(xrho,xdPdrho/P0tot_2piRp,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker2,color1,width2,siz,gca,red,lspace,bspace);
0173
0174 graph1D_jd(rho_dep,rdPdrho_lin/P0tot_2piRp,0,0,'','','',leg,xlim,ylim,style,marker,color2,width2,siz);
0175
0176
0177 figure(4);clf
0178
0179 xlab = 'r/a';
0180 ylab = 'J (MA/m^2)';
0181 xlim = NaN;
0182 ylim = NaN;
0183
0184 graph1D_jd(xrho,opt_disp.si*xJ_rf_dke_fsav_mksa,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker2,color1,width2,siz,gca,red,lspace,bspace);
0185
0186
0187
0188
0189 xlab = 's';
0190 ylab = 'P(s)';
0191 leg = ['linear';...
0192 'LUKE '];
0193
0194 leg2 = ['\theta (\times \pi) ';...
0195 '6.5/T_e^{1/2} ';...
0196 'N_{||} ';...
0197 'N_{\perp} (\times 10)'];
0198
0199 delete([filename,'.txt'])
0200 diary([filename,'.txt'])
0201
0202 for iw = 1:length(waves),
0203
0204 disp(['======> Wave # ',num2str(iw)])
0205 disp([' '])
0206
0207 disp(['Fraction of power absorbed in linear theory : ',num2str(100*absrate.alpha(iw)),' %']);
0208 disp(['Fraction of power absorbed calculated by luke : ',num2str(100*absrate.dke(iw)),' %']);
0209 disp([' '])
0210
0211 for ib = 1:length(waves{iw}.rays),
0212
0213 figure(100*iw + ib)
0214
0215 ray = waves{iw}.rays{ib};
0216
0217
0218 graph1D_jd(ray.ss,ray.sP_2piRp_lin/ray.P0_2piRp,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width2,siz,gca,red,lspace,bspace);
0219 graph1D_jd(ray.zs,ray.zP0_2piRp_mod/ray.P0_2piRp,0,0,'','','',leg,xlim,ylim,style,marker2,color2,width2);
0220
0221 figure(100*iw + 10 + ib)
0222
0223 locylab = ' ';
0224
0225 absmask = ray.isabsmin:ray.isabsmax;
0226
0227 graph1D_jd(ray.ss,[ray.stheta/pi;6.5./sqrt(ray.sTe);ray.sNpar;ray.sNperp/10],0,0,xlab,locylab,tit,leg2,xlim,ylim,style,marker,color0,width2,siz,gca,red,lspace,bspace);
0228 graph1D_jd(ray.ss(absmask),[ray.stheta(absmask)/pi;6.5./sqrt(ray.sTe(absmask));ray.sNpar(absmask);ray.sNperp(absmask)/10],0,0,'','','',NaN,xlim,ylim,style2,marker,color5,width2,siz);
0229
0230 disp(['-----> Absorption of ray # ',num2str(ib)])
0231 disp([' '])
0232
0233 disp(['Fraction of power absorbed in linear theory : ',num2str(100*ray.absrate_alpha),' %']);
0234 disp(['Fraction of power absorbed calculated by luke : ',num2str(100*ray.absrate_dke),' %']);
0235 disp([' '])
0236
0237 disp(['At the peak of power deposition profile'])
0238 disp(['--> Npar : ',num2str(ray.sNpar(ray.ismax)),'']);
0239 disp(['--> Nperp: ',num2str(ray.sNperp(ray.ismax)),'']);
0240 disp(['--> rho : ',num2str(ray.srho(ray.ismax)),'']);
0241 disp(['--> theta: ',num2str(ray.stheta(ray.ismax)/pi),' x pi']);
0242
0243 disp([' '])
0244
0245 end
0246
0247 end
0248
0249 diary off
0250
0251 ir_display = input('radial grid index? ');
0252
0253 while ~isempty(ir_display),
0254
0255 load(filename,'dkeparam','dkedisplay')
0256
0257
0258
0259
0260 [logXfinit_cyl,ppar_cyl,dppar_cyl,pperp_cyl,dpperp_cyl] = s2c_dke_yp(log(XXfinit(:,:,ir_display)),pn,mhu,dkedisplay.dp_cyl);
0261 Xfinit_cyl = exp(logXfinit_cyl);
0262
0263 [logXf0_cyl,ppar_cyl,dppar_cyl,pperp_cyl,dpperp_cyl] = s2c_dke_yp(log(XXf0(:,:,ir_display)),pn,mhu,dkedisplay.dp_cyl);
0264 Xf0_cyl = exp(logXf0_cyl);
0265
0266
0267 XDrf = ZXXD_rf.pp_ij(:,:,ir_display);
0268
0269 [XDrf_cyl,ppar_cyl,dppar_cyl,pperp_cyl,dpperp_cyl] = s2c_dke_yp(XDrf,pn,mhu,dkedisplay.dp_cyl);
0270
0271 Xppar_cyl = ones(length(pperp_cyl),1)*ppar_cyl(:)';
0272 Xpperp_cyl = pperp_cyl(:)*ones(1,length(ppar_cyl));
0273 Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0274
0275
0276
0277
0278 Xf0_cyl = Xf0_cyl.*(Xf0_cyl>0);
0279 Xf0_cyl = Xf0_cyl.*(Xp_cyl<max(ppar_cyl));
0280 Xf0_cyl(isnan(Xf0_cyl)) = 0;
0281
0282 Xfinit_cyl = Xfinit_cyl.*(Xfinit_cyl>0);
0283 Xfinit_cyl = Xfinit_cyl.*(Xp_cyl<max(ppar_cyl));
0284 Xfinit_cyl(isnan(Xfinit_cyl)) = 0;
0285
0286 XDrf_cyl = XDrf_cyl.*(XDrf_cyl>0);
0287 XDrf_cyl = XDrf_cyl.*(Xp_cyl<max(ppar_cyl));
0288 XDrf_cyl(isnan(XDrf_cyl))=0;
0289
0290
0291
0292 figure(10),clf
0293
0294 x1 = ppar_cyl;
0295 xlab1 = 'p_{||}/p_{Te}';
0296 xlim1 = [-dkeparam.pnmax_S,dkeparam.pnmax_S];
0297 x2 = pperp_cyl;
0298 xlab2 = 'p_{\perp}/p_{Te}';
0299 xlim2 = [0,dkeparam.pnmax_S];
0300 leg = ['f_{init}';'f_0 '];
0301
0302 x9 = x1;
0303 y9 = sqrt(dkeparam.pnmax_S^2 - x9.^2);
0304 y1 = Xfinit_cyl;
0305 y2 = Xf0_cyl;
0306
0307 graph2D_jd(x1,x2,y1','','','',xlim1,xlim2,0,betath_ref,[0:0.5:dkeparam.pnmax_S],0,style,color2,width,siz,red,max(max(y1)));
0308 graph2D_jd(x1,x2,y2','','','',xlim1,xlim2,0,betath_ref,[0:0.5:dkeparam.pnmax_S],0,style,color1,width2,siz,1,max(max(y1)));
0309 graph1D_jd(x9,y9,0,0,xlab1,xlab2,'2D Distribution function f_0',NaN,xlim1,xlim2,style,marker,color5,width2,siz);
0310 axis('equal');axis([xlim1,xlim2])
0311 bounce_ripple_jd(dkeparam.bounce_mode,dkeparam.pnmax_S,xmhu0T2,ir_display,[]);
0312
0313 figure(11),clf
0314
0315 y1 = XDrf_cyl;
0316 cont = 10;[0:ceil(max(max(XDrf_cyl)))];
0317
0318 graph2D_jd(x1,x2,y1','','','',xlim1,xlim2,0,NaN,cont,0,style,NaN,width2,siz,red);
0319 graph1D_jd(x9,y9,0,0,xlab1,xlab2,'LH dif. coef. [\nu^{ref}(p_{Te}^{ref})^2]',NaN,xlim1,xlim2,style,marker,color5,width2,siz);
0320 axis('equal');axis([xlim1,xlim2])
0321 bounce_ripple_jd(dkeparam.bounce_mode,dkeparam.pnmax_S,xmhu0T2,ir_display,[]);
0322 colorbar
0323
0324 figure(10),
0325 graph2D_jd(x1,x2,y1',xlab1,xlab2,'2D Distribution function f_0',xlim1,xlim2,0,NaN,cont,0,style,color3,width2,siz);
0326
0327 ir_display = input('radial grid index? ');
0328 end
0329