0001
0002
0003 close all
0004 clear all
0005
0006 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0007
0008 nx = 101;
0009 ixc = (nx + 1)/2;
0010
0011 guess_mode = 1;
0012 opt_mode = 1;
0013
0014 omega_rf = [4.0]*2*pi*1e9;
0015 Npar = 2;
0016
0017 ne = 1e20;
0018 Bfield = 5.0;
0019 Te = mc2*0.01;
0020
0021 zZi = [1,1,1];
0022 zmi = [1,2,3];
0023 zfi = [0,1,0];
0024
0025 ns = 1 + length(zZi);
0026
0027 ne_min = 0;ne_max = 0.02e20;xne = linspace(ne_min,ne_max,nx);xdisp = xne;xlab = 'n_e';xlim = [ne_min,ne_max];
0028 xBfield = Bfield*ones(1,nx);
0029 xTe = Te*ones(1,nx);
0030
0031 xTi = xTe;
0032
0033 xNpar = Npar*ones(1,nx);
0034
0035 zxni = zfi(:)*xne;
0036 zxTi = ones(ns-1,1)*xTi;
0037
0038 sxn = [xne;zxni];
0039 sxT = [xTe;zxTi];
0040
0041 sm = [me;mp*zmi(:)];
0042 sq = qe*[-1;zZi(:)];
0043
0044 sxq = repmat(sq,[1,nx]);
0045 sxm = repmat(sm,[1,nx]);
0046
0047 sxmc2 = mc2*sxm/me;
0048
0049 sxwp = sqrt(sxn.*sxq.^2/e0./sxm);
0050
0051 sxBfield = ones(ns,1)*xBfield;
0052 sxwc = sxq.*sxBfield./sxm;
0053
0054 sxa = (sxwp/omega_rf).^2;
0055 sxy = sxwc/omega_rf;
0056
0057 sxbeta = sqrt(sxT./sxmc2);
0058
0059 xNperpp_c = zeros(1,nx);
0060 xNperpm_c = zeros(1,nx);
0061 xep_c_xyz = zeros(3,nx);
0062 xem_c_xyz = zeros(3,nx);
0063 xep_c_pmz = zeros(3,nx);
0064 xem_c_pmz = zeros(3,nx);
0065 xep_c_pyk = zeros(3,nx);
0066 xem_c_pyk = zeros(3,nx);
0067 xphip_c_xyz = zeros(3,nx);
0068 xphim_c_xyz = zeros(3,nx);
0069 xphip_c_pmz = zeros(3,nx);
0070 xphim_c_pmz = zeros(3,nx);
0071 xphip_c_pyk = zeros(3,nx);
0072 xphim_c_pyk = zeros(3,nx);
0073
0074 if opt_mode >= 0,
0075
0076 xNperpp_k_ft = zeros(1,nx);
0077 xep_k_ft_xyz = zeros(3,nx);
0078 xep_k_ft_pmz = zeros(3,nx);
0079 xep_k_ft_pyk = zeros(3,nx);
0080
0081 xNperpp_k_ht = zeros(1,nx);
0082 xep_k_ht_xyz = zeros(3,nx);
0083 xep_k_ht_pmz = zeros(3,nx);
0084 xep_k_ht_pyk = zeros(3,nx);
0085 xalphaphi_lin_p_k = zeros(1,nx);
0086 xphiTp_k_xyz = zeros(3,nx);
0087 xphiPp_k_xyz = zeros(3,nx);
0088 xphiTp_k_pmz = zeros(3,nx);
0089 xphiPp_k_pmz = zeros(3,nx);
0090 xphiTp_k_pyk = zeros(3,nx);
0091 xphiPp_k_pyk = zeros(3,nx);
0092
0093 end
0094
0095 if opt_mode <= 0,
0096
0097 xNperpm_k_ft = zeros(1,nx);
0098 xem_k_ft_xyz = zeros(3,nx);
0099 xem_k_ft_pmz = zeros(3,nx);
0100 xem_k_ft_pyk = zeros(3,nx);
0101
0102 xNperpm_k_ht = zeros(1,nx);
0103 xem_k_ht_xyz = zeros(3,nx);
0104 xem_k_ht_pmz = zeros(3,nx);
0105 xem_k_ht_pyk = zeros(3,nx);
0106 xalphaphi_lin_m_k = zeros(1,nx);
0107 xphiTm_k_xyz = zeros(3,nx);
0108 xphiPm_k_xyz = zeros(3,nx);
0109 xphiTm_k_pmz = zeros(3,nx);
0110 xphiPm_k_pmz = zeros(3,nx);
0111 xphiTm_k_pyk = zeros(3,nx);
0112 xphiPm_k_pyk = zeros(3,nx);
0113
0114 end
0115
0116 for ix = 1:nx,
0117
0118 [xNperpp_c(ix),xNperpm_c(ix),Kxyz,xep_c_xyz(:,ix),xem_c_xyz(:,ix),xep_c_pmz(:,ix),xem_c_pmz(:,ix),xep_c_pyk(:,ix),xem_c_pyk(:,ix),...
0119 xphip_c_xyz(:,ix),xphim_c_xyz(:,ix),xphip_c_pmz(:,ix),xphim_c_pmz(:,ix),xphip_c_pyk(:,ix),xphim_c_pyk(:,ix)]...
0120 = colddisp_dke_jd(xNpar(ix),omega_rf,xne(ix),zxni(:,ix),zZi,zmi,xBfield(ix));
0121
0122 end
0123
0124 wb = waitbar(0,'dispersion calculations');
0125
0126
0127
0128 for ix = ixc:nx,
0129
0130
0131
0132
0133 if opt_mode >= 0,
0134
0135 if guess_mode == 0 | ix == ixc | isnan(xNperpp_k_ft(ix-1)),
0136 xNperp_guess = xNperpp_c(ix);
0137 else
0138 xNperp_guess = xNperpp_k_ft(ix-1);
0139 end
0140 [xNperpp_k_ft(ix),xep_k_ft_xyz(:,ix),xep_k_ft_pmz(:,ix),xep_k_ft_pyk(:,ix)]...
0141 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,0);
0142
0143 if guess_mode == 0 | ix == ixc | isnan(xNperpp_k_ht(ix-1)),
0144 xNperp_guess = xNperpp_c(ix);
0145 else
0146 xNperp_guess = xNperpp_k_ht(ix-1);
0147 end
0148 [xNperpp_k_ht(ix),xep_k_ht_xyz(:,ix),xep_k_ht_pmz(:,ix),xep_k_ht_pyk(:,ix),D,X,X_s,X_sn,...
0149 xalphaphi_lin_p_k(ix),xphiTp_k_xyz(:,ix),xphiPp_k_xyz(:,ix),xphiTp_k_pmz(:,ix),...
0150 xphiPp_k_pmz(:,ix),xphiTp_k_pyk(:,ix),xphiPp_k_pyk(:,ix)]...
0151 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,1);
0152
0153 end
0154
0155
0156
0157 if opt_mode <= 0,
0158
0159 if guess_mode == 0 | ix == ixc | isnan(xNperpm_k_ft(ix-1)),
0160 xNperp_guess = xNperpm_c(ix);
0161 else
0162 xNperp_guess = xNperpm_k_ft(ix-1);
0163 end
0164 [xNperpm_k_ft(ix),xem_k_ft_xyz(:,ix),xem_k_ft_pmz(:,ix),xem_k_ft_pyk(:,ix)]...
0165 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,0);
0166
0167 if guess_mode == 0 | ix == ixc | isnan(xNperpm_k_ht(ix-1)),
0168 xNperp_guess = xNperpm_c(ix);
0169 else
0170 xNperp_guess = xNperpm_k_ht(ix-1);
0171 end
0172 [xNperpm_k_ht(ix),xem_k_ht_xyz(:,ix),xem_k_ht_pmz(:,ix),xem_k_ht_pyk(:,ix),D,X,X_s,X_sn,...
0173 xalphaphi_lin_m_k(ix),xphiTm_k_xyz(:,ix),xphiPm_k_xyz(:,ix),xphiTm_k_pmz(:,ix),...
0174 xphiPm_k_pmz(:,ix),xphiTm_k_pyk(:,ix),xphiPm_k_pyk(:,ix)]...
0175 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,1);
0176 end
0177
0178 waitbar((ix-ixc)/nx,wb);
0179
0180 end
0181
0182
0183
0184 for ix = ixc:-1:1,
0185
0186
0187
0188 if opt_mode >= 0,
0189
0190 if guess_mode == 0 | ix == ixc | isnan(xNperpp_k_ft(ix+1)),
0191 xNperp_guess = xNperpp_c(ix);
0192 else
0193 xNperp_guess = xNperpp_k_ft(ix+1);
0194 end
0195 [xNperpp_k_ft(ix),xep_k_ft_xyz(:,ix),xep_k_ft_pmz(:,ix),xep_k_ft_pyk(:,ix)]...
0196 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,0);
0197
0198 if guess_mode == 0 | ix == ixc | isnan(xNperpp_k_ht(ix+1)),
0199 xNperp_guess = xNperpp_c(ix);
0200 else
0201 xNperp_guess = xNperpp_k_ht(ix+1);
0202 end
0203 [xNperpp_k_ht(ix),xep_k_ht_xyz(:,ix),xep_k_ht_pmz(:,ix),xep_k_ht_pyk(:,ix),D,X,X_s,X_sn,...
0204 xalphaphi_lin_p_k(ix),xphiTp_k_xyz(:,ix),xphiPp_k_xyz(:,ix),xphiTp_k_pmz(:,ix),...
0205 xphiPp_k_pmz(:,ix),xphiTp_k_pyk(:,ix),xphiPp_k_pyk(:,ix)]...
0206 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,1);
0207
0208 end
0209
0210
0211
0212
0213 if opt_mode <= 0,
0214
0215 if guess_mode == 0 | ix == ixc | isnan(xNperpm_k_ft(ix+1)),
0216 xNperp_guess = xNperpm_c(ix);
0217 else
0218 xNperp_guess = xNperpm_k_ft(ix+1);
0219 end
0220 [xNperpm_k_ft(ix),xem_k_ft_xyz(:,ix),xem_k_ft_pmz(:,ix),xem_k_ft_pyk(:,ix)]...
0221 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,0);
0222
0223 if guess_mode == 0 | ix == ixc | isnan(xNperpm_k_ht(ix+1)),
0224 xNperp_guess = xNperpm_c(ix);
0225 else
0226 xNperp_guess = xNperpm_k_ht(ix+1);
0227 end
0228 [xNperpm_k_ht(ix),xem_k_ht_xyz(:,ix),xem_k_ht_pmz(:,ix),xem_k_ht_pyk(:,ix),D,X,X_s,X_sn,...
0229 xalphaphi_lin_m_k(ix),xphiTm_k_xyz(:,ix),xphiPm_k_xyz(:,ix),xphiTm_k_pmz(:,ix),...
0230 xphiPm_k_pmz(:,ix),xphiTm_k_pyk(:,ix),xphiPm_k_pyk(:,ix)]...
0231 = kineticdisp_jd(sxa(:,ix),sxy(:,ix),sxbeta(:,ix),xNpar(ix),xNperp_guess,1);
0232
0233 end
0234
0235 waitbar((nx - ix + 1)/nx,wb);
0236
0237 end
0238
0239
0240 close(wb)
0241 clear Kxyz D X X_s X_sn xNperp_guess
0242
0243
0244
0245 style1 = '-';
0246 style2 = '--';
0247 marker1 = 'none';
0248 marker2 = '.';
0249 marker3 = 'o';
0250 color0 = NaN;
0251 color1 = 'r';
0252 color2 = 'b';
0253 color3 = 'm';
0254 color4 = 'g';
0255 color5 = 'k';
0256 width1 = 2;
0257 width2 = 0.5;
0258 siz = 20;
0259 red = 0.9;
0260 lspace = 0.7;
0261 bspace = 0.7;
0262
0263 if opt_mode >= 0,
0264
0265 xphip_k_xyz = xphiTp_k_xyz + xphiPp_k_xyz;
0266 xNperpp_k_ht_im = xalphaphi_lin_p_k./xphip_k_xyz(1,:)/2;
0267
0268 figure(1),clf
0269
0270 ylim = NaN;
0271 ylab = 'Re(N_{\perp})';
0272
0273 tit = 'p wave';
0274 leg = ['Cold ';...
0275 'Kinetic FT';...
0276 'Kinetic HT'];
0277
0278 graph1D_jd(xdisp,[real(xNperpp_c);real(xNperpp_k_ft);real(xNperpp_k_ht)],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0279 ylim = get(gca,'ylim');
0280 ylim(1) = 0;
0281 set(gca,'ylim',ylim);
0282
0283 figure(2),clf
0284
0285 ylim = [-ylim(2),0];
0286 ylab = 'Im(N_{\perp})';
0287
0288 graph1D_jd(xdisp,[imag(xNperpp_c);imag(xNperpp_k_ft);real(xNperpp_k_ht_im)],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0289
0290 figure(3),clf
0291
0292 ylim = [0,1];
0293 ylab = '|e_x|';
0294
0295 graph1D_jd(xdisp,[abs(xep_c_xyz(1,:));abs(xep_k_ft_xyz(1,:));abs(xep_k_ht_xyz(1,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0296
0297 figure(4),clf
0298
0299 ylab = '|e_y|';
0300
0301 graph1D_jd(xdisp,[abs(xep_c_xyz(2,:));abs(xep_k_ft_xyz(2,:));abs(xep_k_ht_xyz(2,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0302
0303 figure(5),clf
0304
0305 ylab = '|e_z|';
0306
0307 graph1D_jd(xdisp,[abs(xep_c_xyz(3,:));abs(xep_k_ft_xyz(3,:));abs(xep_k_ht_xyz(3,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0308
0309 figure(6),clf
0310
0311 ylim = NaN;
0312 ylab = 'S_x';
0313
0314 leg = ['Cold ';...
0315 'Kinetic P';...
0316 'Kinetic T'];
0317
0318 graph1D_jd(xdisp,[real(xphip_c_xyz(1,:));real(xphiPp_k_xyz(1,:));real(xphiTp_k_xyz(1,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0319 ylim = get(gca,'ylim');
0320 ylim(2) = 0;
0321 set(gca,'ylim',ylim);
0322
0323 figure(7),clf
0324
0325 ylim = NaN;
0326 ylab = 'S_y';
0327
0328 graph1D_jd(xdisp,[real(xphip_c_xyz(2,:));real(xphiPp_k_xyz(2,:));real(xphiTp_k_xyz(2,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0329
0330
0331
0332
0333 figure(8),clf
0334
0335 ylim = NaN;
0336 ylab = 'S_z';
0337
0338 graph1D_jd(xdisp,[real(xphip_c_xyz(3,:));real(xphiPp_k_xyz(3,:));real(xphiTp_k_xyz(3,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0339 ylim = get(gca,'ylim');
0340 ylim(1) = 0;
0341 set(gca,'ylim',ylim);
0342
0343 end
0344
0345 if opt_mode <= 0,
0346
0347 xphim_k_xyz = xphiTm_k_xyz + xphiPm_k_xyz;
0348 xNperpm_k_ht_im = xalphaphi_lin_m_k./xphim_k_xyz(1,:)/2;
0349
0350 figure(11),clf
0351
0352 ylim = NaN;
0353 ylab = 'Re(N_{\perp})';
0354
0355 tit = 'm wave';
0356 leg = ['Cold ';...
0357 'Kinetic FT';...
0358 'Kinetic HT'];
0359
0360 graph1D_jd(xdisp,[real(xNperpm_c);real(xNperpm_k_ft);real(xNperpm_k_ht)],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0361 ylim = get(gca,'ylim');
0362 ylim(1) = 0;
0363 set(gca,'ylim',ylim);
0364
0365 figure(12),clf
0366
0367 ylab = 'Im(N_{\perp})';
0368
0369 graph1D_jd(xdisp,[imag(xNperpm_c);imag(xNperpm_k_ft);real(xNperpm_k_ht_im)],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0370
0371 figure(13),clf
0372
0373 ylim = [0,1];
0374 ylab = '|e_x|';
0375
0376 graph1D_jd(xdisp,[abs(xem_c_xyz(1,:));abs(xem_k_ft_xyz(1,:));abs(xem_k_ht_xyz(1,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0377
0378 figure(14),clf
0379
0380 ylab = '|e_y|';
0381
0382 graph1D_jd(xdisp,[abs(xem_c_xyz(2,:));abs(xem_k_ft_xyz(2,:));abs(xem_k_ht_xyz(2,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0383
0384 figure(15),clf
0385
0386 ylab = '|e_z|';
0387
0388 graph1D_jd(xdisp,[abs(xem_c_xyz(3,:));abs(xem_k_ft_xyz(3,:));abs(xem_k_ht_xyz(3,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0389
0390 figure(16),clf
0391
0392 ylim = NaN;
0393 ylab = 'S_x';
0394
0395 leg = ['Cold ';...
0396 'Kinetic P';...
0397 'Kinetic T'];
0398
0399 graph1D_jd(xdisp,[real(xphim_c_xyz(1,:));real(xphiPm_k_xyz(1,:));real(xphiTm_k_xyz(1,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0400 ylim = get(gca,'ylim');
0401 ylim(1) = 0;
0402 set(gca,'ylim',ylim);
0403
0404 figure(17),clf
0405
0406 ylim = NaN;
0407 ylab = 'S_y';
0408
0409 graph1D_jd(xdisp,[real(xphim_c_xyz(2,:));real(xphiPm_k_xyz(2,:));real(xphiTm_k_xyz(2,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0410
0411
0412
0413
0414 figure(18),clf
0415
0416 ylim = NaN;
0417 ylab = 'S_z';
0418
0419 graph1D_jd(xdisp,[real(xphim_c_xyz(3,:));real(xphiPm_k_xyz(3,:));real(xphiTm_k_xyz(3,:))],0,0,xlab,ylab,tit,leg,xlim,ylim,style1,marker1,color0,width1,siz,gca,red,lspace,bspace);
0420 ylim = get(gca,'ylim');
0421 ylim(1) = 0;
0422 set(gca,'ylim',ylim);
0423
0424 end
0425
0426