graph_wave_jd

 This function plots the characteristics of wave data

 by J. Decker 02/04/03

0001 function graph_wave_jd(equil,waves,nr_dep,opt,rho_disp,ymask,xvar,yabs,colors,styles)
0002 %
0003 % This function plots the characteristics of wave data
0004 %
0005 % by J. Decker 02/04/03
0006 %
0007 opt.abs_step = 0.01;
0008 %
0009 if nargin < 10
0010     styles = {'-','-','-','-','-','-'};
0011 end
0012 if nargin < 9
0013     colors = {'b','r','g','m','c','y'};
0014 end
0015 if nargin < 8
0016     yabs = NaN;
0017 end
0018 if nargin < 7 || isempty(xvar)
0019     xvar = 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]);
0020 end
0021 if nargin < 6
0022     ymask = NaN;
0023 end
0024 if nargin < 5
0025     rho_disp = 0:0.1:1;
0026 end
0027 if nargin < 4,
0028     opt = struct;
0029 end
0030 if nargin < 3,
0031     nr_dep = 40;
0032 end
0033 %
0034 if ~isfield(opt,'vessel_opt'),
0035     opt.vessel_opt = 'n';
0036 end
0037 if ~isfield(opt,'TeNpar'),
0038     opt.TeNpar = 1;
0039 end
0040 if ~isfield(opt,'abs_max'),
0041     opt.abs_max = input_dke_yp('Option for ray trajectories : 0 (full trajectory), 1 (cut when damped @ 99.9%), 0 <value < 1 (cut when remaining power fraction is below value) ',1,[0;1],'',[1,1]);
0042 end
0043 if ~isfield(opt,'mode'),
0044     opt.mode = input_dke_yp('Option for showing ray polarization (n/y)','n',{'n','y'});
0045 end
0046 if ~isfield(opt,'print'),
0047     opt.print = input_dke_yp('Option for output figures : -1 (nothing), 0 (print), 1 (print and save), 2 (save)',-1,-1:2,'',[1,1]);
0048 end
0049 %
0050 if ~isfield(waves{1},'colldamp'),
0051     if isfield(waves{1}.rayparam,'colldamp'),
0052         waves{1}.colldamp = waves{1}.rayparam.colldamp;%for display options
0053     end
0054 end
0055 %
0056 nw = length(waves);
0057 %
0058 if ~iscell(ymask),
0059     ymask = '';
0060     for iw = 1:nw,
0061         ymask{iw} = 1:length(waves{iw}.rays);
0062     end
0063 end
0064 %
0065 nymask = zeros(1,nw);
0066 %
0067 for iw = 1:nw,
0068     nymask(iw) = length(ymask{iw});
0069     %
0070     if opt.TeNpar && waves{iw}.omega_rf/(2*pi*1e9) < 10,%isfield(waves{iw},'info') && isfield(waves{iw}.info,'ftype') && strcmp(waves{iw}.info.ftype,'LF'),
0071         %
0072         waves{iw}.opt_TeNpar = 1;
0073         %
0074     else
0075         %
0076         waves{iw}.opt_TeNpar = 0;
0077         %
0078     end        
0079 end
0080 %
0081 if ~isfield(opt,'luke'),
0082     opt.luke = false;
0083 end
0084 %
0085 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0086 %
0087 method = 'cubic';
0088 %format
0089 %n_ec = 2;% to be automatized
0090 %N = 5;
0091 %vparn = 6;
0092 %
0093 % Equilibrium
0094 %
0095 Rp = equil.Rp;
0096 Zp = equil.Zp;
0097 ptx = equil.ptx;
0098 pty = equil.pty;
0099 %ptBx = equil.ptBx;
0100 %ptBy = equil.ptBy;
0101 %ptBPHI = equil.ptBPHI;
0102 theta = equil.theta;
0103 psi_apRp = equil.psi_apRp;
0104 %
0105 %[np,nt] = size(ptx);
0106 psin = psi_apRp/psi_apRp(end);
0107 %
0108 ap =  ptx(end,1);%Plasma minor radius on the LFS midplane
0109 %
0110 prho_g = ptx(:,1)/ap;
0111 %
0112 equilDKE = equilibrium_jd(equil);
0113 %
0114 if xvar == 'g',%rho defined as r/a on LFS midplane
0115     prho = prho_g;
0116     plab = 'r/a';
0117 elseif xvar == 'p',%rho defined as square root of poloidal flux
0118     prho = sqrt(psin);
0119     plab = '\rho_{\psi_P}';
0120 elseif xvar == 't',%rho defined as square root of toroidal flux
0121     prho = equilDKE.xrhoT;
0122     plab = '\rho_{\psi_T}';
0123 elseif xvar == 'v',%rho defined as square root of toroidal flux
0124     prho = equilDKE.xrhoV;
0125     plab = '\rho_V';
0126 end
0127 %
0128 %ptB = sqrt(ptBx.^2 + ptBy.^2 + ptBPHI.^2);% total magnetic field
0129 %ptr = sqrt(ptx.^2 + pty.^2); %minor radius on the (psi,stheta) grid
0130 %
0131 ptx_disp = interp1(prho,ptx,rho_disp,method);  
0132 pty_disp = interp1(prho,pty,rho_disp,method); 
0133 %
0134 xmin = min(min(ptx));
0135 xmax = max(max(ptx));
0136 ymin = min(min(pty));
0137 ymax = max(max(pty));
0138 %
0139 PHI_disp = linspace(0,2*pi,361);
0140 X_disp = [(Rp + xmin)*cos(PHI_disp);(Rp + xmax)*cos(PHI_disp)];
0141 Y_disp = [(Rp + xmin)*sin(PHI_disp);(Rp + xmax)*sin(PHI_disp)];
0142 %
0143 Xmin = min(min(X_disp));
0144 Xmax = max(max(X_disp));
0145 Ymin = min(min(Y_disp));
0146 Ymax = max(max(Y_disp));
0147 %
0148 % Wave
0149 %
0150 iy = 0;
0151 rays = '';
0152 legpol = {'SW','FW'};
0153 %
0154 for iw = 1:nw,
0155     ibb = 0;
0156     for ib = ymask{iw},
0157         ibb = ibb + 1;
0158         iy = iy + 1; 
0159         rays{iy} = waves{iw}.rays{ib};
0160         %
0161         rays{iy}.omega_rf = waves{iw}.omega_rf;
0162         rays{iy}.color = colors{mod(iw-1,length(colors))+1};
0163         rays{iy}.style = styles{mod(iw-1,length(styles))+1};
0164         rays{iy}.sX = (Rp + rays{iy}.sx).*cos(-rays{iy}.sphi);
0165         rays{iy}.sY = (Rp + rays{iy}.sx).*sin(-rays{iy}.sphi);
0166         %
0167         rays{iy}.ns = length(rays{iy}.ss);
0168         %
0169         if iscell(yabs),
0170             rays{iy}.mask_abs_lin = 1:rays{iy}.ns >= yabs{iw}(1,ibb) & 1:rays{iy}.ns <= yabs{iw}(2,ibb);  
0171         else
0172             rays{iy}.mask_abs_lin = false(1,rays{iy}.ns);
0173         end
0174         rays{iy}.mask_abs = rays{iy}.mask_abs_lin;
0175         %
0176         if 1,%~isfield(rays{iy},'rayinits') || ~isfield(rays{iy}.rayinits{1},'type') || ~strcmp(rays{iy}.rayinits{1}.type,'EC'),
0177             rays{iy}.refs = '';
0178             rays{iy}.mask_pass = {true(1,rays{iy}.ns)};
0179         else
0180             %
0181             legpol = {'O','X'};
0182             rays{iy}.mask_pass = {};
0183             is0 = 1;
0184             %
0185             for iref = 1:length(rays{iy}.rayinits),% for each reflection
0186                 %
0187                 rayinit = rays{iy}.rayinits{iref};
0188                 %
0189                 launch = rayinit.launch;
0190                 %
0191                 rays{iy}.refs{iref}.x = launch.yR_L - Rp;
0192                 rays{iy}.refs{iref}.y = launch.yZ_L - Zp;
0193                 rays{iy}.refs{iref}.phi = launch.yphi_L;
0194                 %
0195                 while isfield(launch,'launch_back') && ~isempty(launch.launch_back),
0196                     %
0197                     launch = launch.launch_back;
0198                     %
0199                     rays{iy}.refs{iref}.x = [launch.yR_L - Rp,rays{iy}.refs{iref}.x];
0200                     rays{iy}.refs{iref}.y = [launch.yZ_L - Zp,rays{iy}.refs{iref}.y];
0201                     rays{iy}.refs{iref}.phi = [launch.yphi_L,rays{iy}.refs{iref}.phi];
0202                     %
0203                 end
0204                 %
0205                 is = find(rays{iy}.ss == rayinit.ss0);
0206                 %
0207                 if is > 1 && is < length(rays{iy}.ss),
0208                     %
0209                     rays{iy}.refs{iref}.x = [rays{iy}.sx(is),rays{iy}.refs{iref}.x,rays{iy}.sx(is+1)];
0210                     rays{iy}.refs{iref}.y = [rays{iy}.sy(is),rays{iy}.refs{iref}.y,rays{iy}.sy(is+1)];
0211                     rays{iy}.refs{iref}.phi = [rays{iy}.sphi(is),rays{iy}.refs{iref}.phi,rays{iy}.sphi(is+1)];
0212                     %
0213                     rays{iy}.mask_pass = [rays{iy}.mask_pass,1:rays{iy}.ns >= is0 & 1:rays{iy}.ns <= is];
0214                     is0 = is + 1;
0215                     %
0216                 elseif is == 1,
0217                     %
0218                     rays{iy}.refs{iref}.x = [rays{iy}.refs{iref}.x,rays{iy}.sx(is)];
0219                     rays{iy}.refs{iref}.y = [rays{iy}.refs{iref}.y,rays{iy}.sy(is)];
0220                     rays{iy}.refs{iref}.phi = [rays{iy}.refs{iref}.phi,rays{iy}.sphi(is)];
0221                     %
0222                 else
0223                     %
0224                     rays{iy}.refs{iref}.x = [rays{iy}.sx(is),rays{iy}.refs{iref}.x];
0225                     rays{iy}.refs{iref}.y = [rays{iy}.sy(is),rays{iy}.refs{iref}.y];
0226                     rays{iy}.refs{iref}.phi = [rays{iy}.sphi(is),rays{iy}.refs{iref}.phi];
0227                     %
0228                     rays{iy}.mask_pass = [rays{iy}.mask_pass,1:rays{iy}.ns >= is0 & 1:rays{iy}.ns <= is];
0229                     %
0230                 end
0231                 %
0232                 rays{iy}.refs{iref}.R = rays{iy}.refs{iref}.x + Rp;
0233                 rays{iy}.refs{iref}.Z = rays{iy}.refs{iref}.y + Zp;
0234                 rays{iy}.refs{iref}.X = rays{iy}.refs{iref}.R.*cos(-rays{iy}.refs{iref}.phi);
0235                 rays{iy}.refs{iref}.Y = rays{iy}.refs{iref}.R.*sin(-rays{iy}.refs{iref}.phi);
0236                 %
0237                 xmin = min([xmin,rays{iy}.refs{iref}.x]);    
0238                 xmax = max([xmax,rays{iy}.refs{iref}.x]);    
0239                 ymin = min([ymin,rays{iy}.refs{iref}.y]);    
0240                 ymax = max([ymax,rays{iy}.refs{iref}.y]);    
0241                 %
0242                 Xmin = min([Xmin,rays{iy}.refs{iref}.X]);    
0243                 Xmax = max([Xmax,rays{iy}.refs{iref}.X]);    
0244                 Ymin = min([Ymin,rays{iy}.refs{iref}.Y]);    
0245                 Ymax = max([Ymax,rays{iy}.refs{iref}.Y]);    
0246                 %
0247             end
0248             %
0249         end
0250         %
0251     end
0252 end
0253 ny = iy;
0254 %
0255 % Damping
0256 %
0257 yrho_abs_lin = NaN(1,ny);
0258 yPtot_abs_lin = NaN(1,ny);
0259 yNpar_abs_lin = NaN(1,ny);
0260 %
0261 yrho_abs_lin_coll = NaN(1,ny);
0262 yPtot_abs_lin_coll = NaN(1,ny);
0263 yNpar_abs_lin_coll = NaN(1,ny);
0264 %
0265 yrho_abs_coll_luke = NaN(1,ny);
0266 yPtot_abs_coll_luke = NaN(1,ny);
0267 yNpar_abs_coll_luke = NaN(1,ny);
0268 %
0269 yrho_abs = NaN(1,ny);
0270 yPtot_abs = NaN(1,ny);
0271 yNpar_abs = NaN(1,ny);
0272 %
0273 dPdrho_lin = 0;
0274 dPdrho_lin_coll = 0;
0275 dPdrho_coll_luke = 0;
0276 dPdrho = 0;
0277 %
0278 if opt.abs_max == 1,
0279     opt.abs_max = 1e-3;% default value
0280 end
0281 %
0282 if ~isfield(opt,'disp_rays'),
0283     opt.disp_rays = input_dke_yp('Option for showing ray absorption (n/y)','n',{'n','y'});
0284 end
0285 %
0286 for iy = 1:ny,
0287     %
0288     ns = length(rays{iy}.ss);
0289     %
0290     rays{iy}.is_abs_max_lin = ns;
0291     rays{iy}.is_abs_max_lin_coll = ns;
0292     rays{iy}.is_abs_max_coll_luke = ns;
0293     rays{iy}.is_abs_lin = NaN;
0294     rays{iy}.dP_lin_normdrho = NaN;
0295     rays{iy}.dP_lin_coll_normdrho = NaN;
0296     rays{iy}.dP_coll_luke_normdrho = NaN;
0297     yrho_abs_lin(iy) = NaN;
0298     yNpar_abs_lin(iy) = NaN;
0299     yPtot_abs_lin(iy) = NaN;
0300     yrho_abs_lin_coll(iy) = NaN;
0301     yNpar_abs_lin_coll(iy) = NaN;
0302     yPtot_abs_lin_coll(iy) = NaN;
0303     yrho_abs_coll_luke(iy) = NaN;
0304     yNpar_abs_coll_luke(iy) = NaN;
0305     yPtot_abs_coll_luke(iy) = NaN;
0306     %
0307     rhog_dep = NaN;
0308     %
0309     if opt.disp_rays == 'y',
0310         figure(1000+iy),clf;%ray figures
0311     else
0312         if ishandle(figure(1000+iy)),
0313             close(figure(1000+iy));
0314         end
0315     end
0316     %
0317     if isfield(rays{iy},'sP_2piRp_lin'),%linear resonant Landau absorption (RLA)
0318         %
0319         sP_2piRp_lin_norm = rays{iy}.sP_2piRp_lin/rays{iy}.sP_2piRp_lin(1);
0320         sdP_2piRp_ds_lin_norm = rays{iy}.sdP_2piRp_ds_lin/max(rays{iy}.sdP_2piRp_ds_lin+eps).*(1 - sP_2piRp_lin_norm(end)/sP_2piRp_lin_norm(1));
0321         if ~iscell(yabs),
0322             rays{iy}.mask_abs_lin = sdP_2piRp_ds_lin_norm > opt.abs_step;
0323         end
0324         if opt.abs_max > 0,
0325             rays{iy}.is_abs_max_lin = min([ns,find(sP_2piRp_lin_norm < opt.abs_max)]);        
0326         else
0327             rays{iy}.is_abs_max_lin = ns;
0328         end
0329         %
0330         if isfield(opt,'mode') && opt.mode == 'y',
0331 %             rays{iy}.mask_p_lin = rays{iy}.spmode == 1;
0332 %             rays{iy}.mask_m_lin = rays{iy}.spmode == 0;
0333 %             %
0334 %             rays{iy}.mask_p_lin(rays{iy}.is_abs_max_lin+1:ns) = false;
0335 %             rays{iy}.mask_m_lin(rays{iy}.is_abs_max_lin+1:ns) = false;
0336 %             %
0337 %             rays{iy}.mask_p_lin(1) = true;%to avoid legend problems
0338 %             rays{iy}.mask_m_lin(1) = true;
0339             %
0340             spmode = rays{iy}.spmode(1:rays{iy}.is_abs_max_lin);
0341             %
0342             rays{iy}.mask_p_lin{1} = 1:rays{iy}.ns == 1;%to avoid legend problems
0343             rays{iy}.mask_m_lin{1} = 1:rays{iy}.ns == 1;%to avoid legend problems
0344             %
0345             imc_list = 1+[0,find(spmode(2:end) ~= spmode(1:end-1)),rays{iy}.is_abs_max_lin];
0346             %
0347             pmode = spmode(1);
0348             %
0349             for imc = 1:length(imc_list)-1,
0350                 %
0351                 mask_mc = 1:rays{iy}.ns >= imc_list(imc) & 1:rays{iy}.ns < imc_list(imc+1);
0352                 %
0353                 if pmode == 1,
0354                     rays{iy}.mask_p_lin = [rays{iy}.mask_p_lin,mask_mc];
0355                 else
0356                     rays{iy}.mask_m_lin = [rays{iy}.mask_m_lin,mask_mc];
0357                 end
0358                 %
0359                 pmode = ~pmode;
0360                 %
0361             end
0362             %
0363         else
0364             %
0365             rays{iy}.mask_p_lin = NaN;
0366             rays{iy}.mask_m_lin = NaN;               
0367             %
0368         end
0369         %
0370         rays{iy}.is_abs_lin = min(find(sdP_2piRp_ds_lin_norm == max(sdP_2piRp_ds_lin_norm))); 
0371         %
0372         yrho_abs_lin(iy) = rays{iy}.srho(rays{iy}.is_abs_lin);
0373         yNpar_abs_lin(iy) = real(rays{iy}.sNpar(rays{iy}.is_abs_lin));
0374         yPtot_abs_lin(iy) = 1 - sP_2piRp_lin_norm(end);
0375         %
0376         [rays{iy}.dP_lin_normdrho,rhog_dep] = radialdampingprofile_jd(rays{iy}.srho,rays{iy}.sP_2piRp_lin/rays{iy}.sP_2piRp_lin(1),nr_dep);
0377         %
0378         dPdrho_lin = dPdrho_lin + rays{iy}.dP_lin_normdrho*rays{iy}.P0_2piRp*2*pi*Rp;
0379         %
0380         if opt.disp_rays == 'y',
0381             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0382                 graph1D_jd(rays{iy}.ss,rays{iy}.sP_2piRp_lin*2*pi*equil.Rp,0,0,'ray length (m)','Power (W)',['ray #',int2str(iy),', N_{||0} = ',num2str(real(rays{iy}.sNpar(1)))],NaN,'','','--','none','r',2);
0383             else
0384                 graph1D_jd(rays{iy}.ss,rays{iy}.sP_2piRp_lin*2*pi*equil.Rp,0,0,'ray length (m)','Power (W)',['ray #',int2str(iy),', N_{||0} = ',num2str(real(rays{iy}.sNpar(1)))],NaN,'','','-','none','r',2);
0385             end
0386             %
0387             legend('Linear absorption (RLA)','Location','SouthWest')
0388         end
0389         %
0390     end
0391     %
0392     if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,%linear non-resonant collisional absorption (NRCA) + linear resonant Landau absorption (RLA)
0393         %
0394         sP_2piRp_lin_coll_norm = rays{iy}.sP_2piRp_lin_coll/rays{iy}.sP_2piRp_lin_coll(1);
0395         sdP_2piRp_ds_lin_coll_norm = rays{iy}.sdP_2piRp_ds_lin_coll/max(rays{iy}.sdP_2piRp_ds_lin_coll+eps).*(1 - sP_2piRp_lin_coll_norm(end)/sP_2piRp_lin_coll_norm(1));
0396         if ~iscell(yabs),
0397             rays{iy}.mask_abs_lin_coll = sdP_2piRp_ds_lin_coll_norm > opt.abs_step;
0398         end
0399         if opt.abs_max > 0,
0400             rays{iy}.is_abs_max_lin_coll = min([ns,find(sP_2piRp_lin_coll_norm < opt.abs_max)]);        
0401         else
0402             rays{iy}.is_abs_max_lin_coll = ns;
0403         end
0404         %
0405         if isfield(opt,'mode') && opt.mode == 'y',
0406 %             rays{iy}.mask_p = rays{iy}.spmode == 1;
0407 %             rays{iy}.mask_m = rays{iy}.spmode == 0;
0408 %             %
0409 %             rays{iy}.mask_p(rays{iy}.is_abs_max+1:ns) = false;
0410 %             rays{iy}.mask_m(rays{iy}.is_abs_max+1:ns) = false;
0411 %             %
0412 %             rays{iy}.mask_p(1) = true;%to avoid legend problems
0413 %             rays{iy}.mask_m(1) = true;
0414             %
0415             spmode = rays{iy}.spmode(1:rays{iy}.is_abs_max_lin_coll);
0416             %
0417             rays{iy}.mask_p_lin_coll{1} = 1;%to avoid legend problems
0418             rays{iy}.mask_m_lin_coll{1} = 1;%to avoid legend problems
0419             %
0420             imc_list = 1+[0,find(spmode(2:end) ~= spmode(1:end-1)),rays{iy}.is_abs_max_lin_coll];
0421             %
0422             pmode = spmode(1);
0423             %
0424             for imc = 1:length(imc_list)-1,
0425                 %
0426                 mask_mc = imc_list(imc):(imc_list(imc+1) - 1);
0427                 %
0428                 if pmode == 1,
0429                     rays{iy}.mask_p_lin_coll = [rays{iy}.mask_p_lin_coll,mask_mc];
0430                 else
0431                     rays{iy}.mask_m_lin_coll = [rays{iy}.mask_m_lin_coll,mask_mc];
0432                 end
0433                 %
0434                 pmode = ~pmode;
0435                 %
0436             end
0437             %
0438         else
0439             %
0440             rays{iy}.mask_p_lin_coll = NaN;
0441             rays{iy}.mask_m_lin_coll = NaN;               
0442             %
0443         end
0444         %
0445         rays{iy}.is_abs_lin_coll = min(find(sdP_2piRp_ds_lin_coll_norm == max(sdP_2piRp_ds_lin_coll_norm))); 
0446         %
0447         yrho_abs_lin_coll(iy) = rays{iy}.srho(rays{iy}.is_abs_lin_coll);
0448         yNpar_abs_lin_coll(iy) = real(rays{iy}.sNpar(rays{iy}.is_abs_lin_coll));
0449         yPtot_abs_lin_coll(iy) = 1 - sP_2piRp_lin_coll_norm(end);
0450         %
0451         [rays{iy}.dP_lin_coll_normdrho,rhog_dep] = radialdampingprofile_jd(rays{iy}.srho,rays{iy}.sP_2piRp_lin_coll/rays{iy}.sP_2piRp_lin_coll(1),nr_dep);
0452         %
0453         dPdrho_lin_coll = dPdrho_lin_coll + rays{iy}.dP_lin_coll_normdrho*rays{iy}.P0_2piRp*2*pi*Rp;
0454         %
0455         if opt.disp_rays == 'y',
0456             hold on,graph1D_jd(rays{iy}.ss,rays{iy}.sP_2piRp_lin_coll*2*pi*equil.Rp,0,0,'ray length (m)','Power (W)',['ray #',int2str(iy),', N_{||0} = ',num2str(real(rays{iy}.sNpar(1)))],NaN,'','','-','none','r',2);
0457             %
0458             legend('Linear absorption (RLA)','Linear absorption (RLA + NRCA)','Location','SouthWest');
0459         end
0460         %
0461     else
0462         %
0463         rays{iy}.is_abs_max_lin_coll = ns;
0464         rays{iy}.is_abs_lin_coll = NaN;
0465         rays{iy}.dP_normdrho_lin_coll = NaN;
0466         yrho_abs_lin_coll(iy) = NaN;
0467         yNpar_abs_lin_coll(iy) = NaN;
0468         yPtot_abs_lin_coll(iy) = NaN;
0469         %
0470     end    
0471     %
0472     if opt.luke && isfield(rays{iy},'sP_2piRp'),%LUKE absorption (without non-resonant collisional absorption if any)
0473         sP_2piRp_norm = rays{iy}.sP_2piRp{end}/rays{iy}.sP_2piRp{end}(1);
0474         sdP_2piRp_ds_norm = rays{iy}.sdP_2piRp_ds{end}/max(rays{iy}.sdP_2piRp_ds{end}+eps).*(1 - sP_2piRp_norm(end)/sP_2piRp_norm(1));
0475         if ~iscell(yabs),
0476             rays{iy}.mask_abs = sdP_2piRp_ds_norm > opt.abs_step;
0477         end
0478         if opt.abs_max > 0,
0479             rays{iy}.is_abs_max = min([ns,find(sP_2piRp_norm < opt.abs_max)]);        
0480         else
0481             rays{iy}.is_abs_max = ns;
0482         end
0483         %
0484         if isfield(opt,'mode') && opt.mode == 'y',
0485 %             rays{iy}.mask_p = rays{iy}.spmode == 1;
0486 %             rays{iy}.mask_m = rays{iy}.spmode == 0;
0487 %             %
0488 %             rays{iy}.mask_p(rays{iy}.is_abs_max+1:ns) = false;
0489 %             rays{iy}.mask_m(rays{iy}.is_abs_max+1:ns) = false;
0490 %             %
0491 %             rays{iy}.mask_p(1) = true;%to avoid legend problems
0492 %             rays{iy}.mask_m(1) = true;
0493             %
0494             spmode = rays{iy}.spmode(1:rays{iy}.is_abs_max);
0495             %
0496             rays{iy}.mask_p{1} = 1;%to avoid legend problems
0497             rays{iy}.mask_m{1} = 1;%to avoid legend problems
0498             %
0499             imc_list = 1+[0,find(spmode(2:end) ~= spmode(1:end-1)),rays{iy}.is_abs_max];
0500             %
0501             pmode = spmode(1);
0502             %
0503             for imc = 1:length(imc_list)-1,
0504                 %
0505                 mask_mc = imc_list(imc):(imc_list(imc+1) - 1);
0506                 %
0507                 if pmode == 1,
0508                     rays{iy}.mask_p = [rays{iy}.mask_p,mask_mc];
0509                 else
0510                     rays{iy}.mask_m = [rays{iy}.mask_m,mask_mc];
0511                 end
0512                 %
0513                 pmode = ~pmode;
0514                 %
0515             end
0516             %
0517         else
0518             %
0519             rays{iy}.mask_p = NaN;
0520             rays{iy}.mask_m = NaN;               
0521             %
0522         end
0523         %
0524         rays{iy}.is_abs = min(find(sdP_2piRp_ds_norm == max(sdP_2piRp_ds_norm))); 
0525         %
0526         yrho_abs(iy) = rays{iy}.srho(rays{iy}.is_abs);
0527         yNpar_abs(iy) = real(rays{iy}.sNpar(rays{iy}.is_abs));
0528         yPtot_abs(iy) = 1 - sP_2piRp_norm(end);
0529         %
0530         rays{iy}.dP_normdrho = radialdampingprofile_jd(rays{iy}.srho,rays{iy}.sP_2piRp{end}/rays{iy}.sP_2piRp{end}(1),nr_dep);
0531         %
0532         dPdrho = dPdrho + rays{iy}.dP_normdrho*rays{iy}.P0_2piRp*2*pi*Rp;
0533         %
0534         if opt.disp_rays == 'y',
0535             if ~(isfield(waves{1},'colldamp') && waves{1}.colldamp == 1),
0536                 hold on,graph1D_jd(rays{iy}.ss,rays{iy}.sP_2piRp{1}*2*pi*equil.Rp,0,0,'ray length (m)','Power (W)',['ray #',int2str(iy),', N_{||0} = ',num2str(real(rays{iy}.sNpar(1)))],NaN,'','','--','none','k',2);
0537                 legend('Linear absorption (RLA)','LUKE quasilinear absorption (RLA)','Location','SouthWest');
0538             end
0539         end
0540         %
0541     else
0542         %
0543         rays{iy}.is_abs_max = ns;
0544         rays{iy}.is_abs = NaN;
0545         rays{iy}.dP_normdrho = NaN;
0546         yrho_abs(iy) = NaN;
0547         yNpar_abs(iy) = NaN;
0548         yPtot_abs(iy) = NaN;
0549         %
0550     end    
0551     %
0552     if opt.luke && isfield(waves{1},'colldamp') && waves{1}.colldamp == 1 && isfield(rays{iy},'sP_2piRp_coll_luke'),%LUKE absorption (with non-resonant collisional absorption if any)
0553         sP_2piRp_coll_luke_norm = rays{iy}.sP_2piRp_coll_luke{end}/rays{iy}.sP_2piRp_coll_luke{end}(1);
0554         sdP_2piRp_coll_luke_ds_norm = rays{iy}.sdP_2piRp_coll_luke_ds{end}/max(rays{iy}.sdP_2piRp_coll_luke_ds{end}+eps).*(1 - sP_2piRp_coll_luke_norm(end)/sP_2piRp_coll_luke_norm(1));
0555         if ~iscell(yabs),
0556             rays{iy}.mask_abs_coll_luke = sdP_2piRp_coll_luke_ds_norm > opt.abs_step;
0557         end
0558         if opt.abs_max > 0,
0559             rays{iy}.is_abs_max_coll_luke = min([ns,find(sP_2piRp_coll_luke_norm < opt.abs_max)]);        
0560         else
0561             rays{iy}.is_abs_max_coll_luke = ns;
0562         end
0563         %
0564         if isfield(opt,'mode') && opt.mode == 'y',
0565 %             rays{iy}.mask_p = rays{iy}.spmode == 1;
0566 %             rays{iy}.mask_m = rays{iy}.spmode == 0;
0567 %             %
0568 %             rays{iy}.mask_p(rays{iy}.is_abs_max+1:ns) = false;
0569 %             rays{iy}.mask_m(rays{iy}.is_abs_max+1:ns) = false;
0570 %             %
0571 %             rays{iy}.mask_p(1) = true;%to avoid legend problems
0572 %             rays{iy}.mask_m(1) = true;
0573             %
0574             spmode = rays{iy}.spmode(1:rays{iy}.is_abs_max_coll_luke);
0575             %
0576             rays{iy}.mask_p_coll_luke{1} = 1;%to avoid legend problems
0577             rays{iy}.mask_m_coll_luke{1} = 1;%to avoid legend problems
0578             %
0579             imc_list = 1+[0,find(spmode(2:end) ~= spmode(1:end-1)),rays{iy}.is_abs_max_coll_luke];
0580             %
0581             pmode = spmode(1);
0582             %
0583             for imc = 1:length(imc_list)-1,
0584                 %
0585                 mask_mc = imc_list(imc):(imc_list(imc+1) - 1);
0586                 %
0587                 if pmode == 1,
0588                     rays{iy}.mask_p_coll_luke = [rays{iy}.mask_p_coll_luke,mask_mc];
0589                 else
0590                     rays{iy}.mask_m_coll_luke = [rays{iy}.mask_m_coll_luke,mask_mc];
0591                 end
0592                 %
0593                 pmode = ~pmode;
0594                 %
0595             end
0596             %
0597         else
0598             %
0599             rays{iy}.mask_p_coll_luke = NaN;
0600             rays{iy}.mask_m_coll_luke = NaN;               
0601             %
0602         end
0603         %
0604         rays{iy}.is_abs_coll_luke = min(find(sdP_2piRp_coll_luke_ds_norm == max(sdP_2piRp_coll_luke_ds_norm))); 
0605         %
0606         yrho_abs_coll_luke(iy) = rays{iy}.srho(rays{iy}.is_abs_coll_luke);
0607         yNpar_abs_coll_luke(iy) = real(rays{iy}.sNpar(rays{iy}.is_abs_coll_luke));
0608         yPtot_abs_coll_luke(iy) = 1 - sP_2piRp_coll_luke_norm(end);
0609         %
0610         [rays{iy}.dP_coll_luke_normdrho,rhog_dep] = radialdampingprofile_jd(rays{iy}.srho,rays{iy}.sP_2piRp_coll_luke{end}/rays{iy}.sP_2piRp_coll_luke{end}(1),nr_dep);
0611         %
0612         rays{iy}.dP_coll_luke_normdrho = radialdampingprofile_jd(rays{iy}.srho,rays{iy}.sP_2piRp_coll_luke{end}/rays{iy}.sP_2piRp_coll_luke{end}(1),nr_dep);
0613         %
0614         dPdrho_coll_luke = dPdrho_coll_luke + rays{iy}.dP_coll_luke_normdrho*rays{iy}.P0_2piRp*2*pi*Rp;
0615         %
0616         if opt.disp_rays == 'y',
0617             hold on,graph1D_jd(rays{iy}.ss,rays{iy}.sP_2piRp_coll_luke{1}*2*pi*equil.Rp,0,0,'ray length (m)','Power (W)',['ray #',int2str(iy),', N_{||0} = ',num2str(real(rays{iy}.sNpar(1)))],NaN,'','','-','none','k',2);
0618             %
0619             legend('Linear absorption (RLA)','Linear absorption (RLA + NRCA)','LUKE quasilinear absorption (RLA + NRCA)','Location','SouthWest');
0620         end
0621         %
0622     else
0623         %
0624         rays{iy}.sP_2piRp_coll_luke = NaN;
0625         rays{iy}.is_abs_max_coll_luke = ns;
0626         rays{iy}.is_abs_coll_luke = NaN;
0627         rays{iy}.dP_normdrho_coll_luke = NaN;
0628         rays{iy}.mask_abs_coll_luke = NaN;
0629         yrho_abs_coll_luke(iy) = NaN;
0630         yNpar_abs_coll_luke(iy) = NaN;
0631         yPtot_abs_coll_luke(iy) = NaN;
0632         
0633         %
0634     end    
0635     %
0636     if opt.disp_rays == 'y',
0637         if opt.luke,
0638             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0639                 disp(['Ray # ',num2str(iy),': fraction of power absorbed (LUKE: quasilinear, RLA + NRCA): ',num2str(round(yPtot_abs_coll_luke(iy)*100)),' %'])
0640             else     
0641                 disp(['Ray # ',num2str(iy),': fraction of power absorbed (LUKE: quasilinear, RLA): ',num2str(round(yPtot_abs(iy)*100)),' %']) 
0642             end
0643         end
0644         %
0645         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0646             disp(['Ray # ',num2str(iy),': fraction of power absorbed (C3PO: linear, RLA + NRCA): ',num2str(round(yPtot_abs_lin_coll(iy)*100)),' %'])  
0647         else
0648             disp(['Ray # ',num2str(iy),': fraction of power absorbed (C3PO: linear, RLA): ',num2str(round(yPtot_abs_lin(iy)*100)),' %'])  
0649         end
0650     end
0651     %
0652 end
0653 %
0654 if sum(dPdrho) == 0,
0655     opt.luke = false;
0656 end
0657 %
0658 rho_dep = interp1(prho_g,prho,rhog_dep);
0659 %
0660 drhog_dep = diff([0,(rhog_dep(1:end-1) + rhog_dep(2:end))/2,1]);
0661 drho_dep = diff([0,(rho_dep(1:end-1) + rho_dep(2:end))/2,1]);
0662 %
0663 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0664     dPdrho_lin = dPdrho_lin.*drhog_dep./drho_dep;
0665     dPdrho_lin_coll = dPdrho_lin_coll.*drhog_dep./drho_dep;
0666     dPdrho_coll_luke = dPdrho_coll_luke.*drhog_dep./drho_dep;
0667 else
0668     dPdrho_lin = dPdrho_lin.*drhog_dep./drho_dep;
0669     dPdrho_lin_coll = zeros(size(dPdrho_lin));
0670     dPdrho_coll_luke = zeros(size(dPdrho_lin));
0671 end
0672 %
0673 dPdrho = dPdrho.*drhog_dep./drho_dep;    
0674 %
0675 % %
0676 % % Resonance
0677 % %
0678 % if opt(2) >= 1 | opt(3) >= 1 | opt(4) >= 1,
0679 %     if isempty(sypsin) | isempty(sytheta),
0680 %         for ib = 1:ny
0681 %         %
0682 %         % ray position in (r,theta) coordinates
0683 %         %
0684 %             sR = syR(:,ib) + eps*(syR(:,ib) == Rp); %to avoid singularities
0685 %             sZ = syZ(:,ib) + eps*(syZ(:,ib) == Zp); %to avoid singularities
0686 %             sr = sqrt((sR - Rp).^2 + (sZ - Zp).^2); %minor radius
0687 %             stheta = atan((sZ - Zp)./(sR - Rp)) + pi - pi*(sR > Rp).*sign(sZ - Zp); %poloidal angle [0,2*pi]
0688 %         %
0689 %         % interpolation for r at ray poloidal locations
0690 %         %
0691 %             psr = interp1(theta,r',stheta,method)';%r(psi,st)
0692 %         %
0693 %         % interpolation for psi at ray position (along cst theta lines)
0694 %         %
0695 %         %spsin = interp1(pstr,ppsin,sr,method);%psi(s)
0696 %             method2 = 'linear';
0697 %             for is = 1:ns
0698 %                 spsin(is,1) = interp1(psr(:,is),ppsin,sr(is),method2);%psi(s)
0699 %             end
0700 %         %
0701 %             imask_1 = min(find(~isnan(spsin)));
0702 %             imask_2 = min(find(isnan(spsin(imask_1:ns))));
0703 %             mask = [1:(imask_1 - 1),(imask_1 + imask_2 - 1):ns];
0704 %         %
0705 %             spsin(mask) = 1;%part of ray outside the last flux-surface
0706 %         %
0707 %             sypsin(1:ns,ib) = spsin;
0708 %             sytheta(1:ns,ib) = stheta;
0709 %         end
0710 %     end
0711 %     for ib = 1:ny
0712 %     %
0713 %         spsin = sypsin(:,ib);
0714 %         stheta = sytheta(:,ib);
0715 %     %
0716 %         syB(:,ib) = interp2(theta,ppsin,B,stheta,spsin,method);
0717 %         if opt(4) >= 1,
0718 %             syBR(:,ib) = interp2(theta,ppsin,BR,stheta,spsin,method);
0719 %             syBZ(:,ib) = interp2(theta,ppsin,BZ,stheta,spsin,method);
0720 %             syBPHI(:,ib) = interp2(theta,ppsin,BPHI,stheta,spsin,method);
0721 %         end
0722 %     %
0723 %         syne(:,ib) = interp1(ppsin,equil.pne,spsin,method);
0724 %         siyzni(:,:,ib) = interp1(ppsin,equil.pzni',spsin,method);
0725 %         syTe(:,ib) = interp1(ppsin,equil.pTe,spsin,method);
0726 %         siyzTi(:,:,ib) = interp1(ppsin,equil.pzTi',spsin,method);
0727 %     %
0728 %     end
0729 % %
0730 %   sybetath = sqrt(syTe/mc2);
0731 % %
0732 %   omega_ce = qe*B/me;
0733 %   yres = omega_rf./omega_ce;
0734 % %
0735 %   syomega_ce = qe*syB/me;
0736 %   syres = omega_rf./syomega_ce;
0737 %   is_syres = find(abs(diff(floor(syres))) == 1);
0738 %   if exist('wave_s','var') == 1,
0739 % %
0740 %     sybetath_s = sqrt(syTe_s/mc2);
0741 % %
0742 %     syomega_ce_s = qe*syB_s/me;
0743 %     syres_s = omega_rf./syomega_ce_s;
0744 % %
0745 %   end
0746 % end
0747 % %
0748 % % Absorption coefficient
0749 % %
0750 % if opt(3) >= 1,
0751 % %
0752 %     if isempty(syepol_pmz) | isempty(syphi_pmz),
0753 %         for ib = 1:ny
0754 %             for is = 1:ns
0755 %               if ~isnan(syNpar(is,ib))
0756 %                 [sNperpp,sNperpm,sKxyz,sep_xyz,sem_xyz,sep_pmz,sem_pmz,sep_pyk,sem_pyk,sphip_xyz,sphim_xyz,sphip_pmz,sphim_pmz] = ...
0757 %                     colddisp_dke_jd(syNpar(is,ib),omega_rf,syne(is,ib),siyzni(is,:,ib),equil.zZi,equil.zmi,syB(is,ib));
0758 %                 maskp = (abs(sNperpp - syNperp(is,ib)) <= abs(sNperpm - syNperp(is,ib)));
0759 %                 syNperp_cold(is,ib) = sNperpp.*maskp + sNperpm.*~maskp;
0760 %                 syepol_pmz(is,1:3,ib) = sep_pmz'.*maskp + sem_pmz'.*~maskp;
0761 %                 syphi_pmz(is,1:3,ib) = sphip_pmz'.*maskp + sphim_pmz'.*~maskp;
0762 %               else
0763 %                 syNperp_cold(is,ib) = NaN;
0764 %                 syepol_pmz(is,1:3,ib) = NaN;
0765 %                 syphi_pmz(is,1:3,ib) = NaN;
0766 %               end
0767 %             end
0768 %         end
0769 %     end
0770 %     %
0771 %     if id_wave(1) == 'S' & exist('syNperp_cold'),
0772 %         syNperp = syNperp_cold;
0773 %     end
0774 %     %
0775 %     sywpe = sqrt(syne*qe^2/e0/me);
0776 %     syabsphi = squeeze(sqrt(abs(syphi_pmz(:,1,:)).^2 + abs(syphi_pmz(:,2,:)).^2 + abs(syphi_pmz(:,3,:)).^2));
0777 %     %
0778 %     syalpha = 2^(1-n_ec)*n_ec/sqrt(2*pi)/fact_dke_jd(n_ec - 1).*sywpe.^2.*syres.^(2*(n_ec - 1))...
0779 %         .*sybetath.^(2*n_ec - 3)*pi/clum/omega_rf./syabsphi.*syNperp.^(2*(n_ec - 1)).*squeeze(abs(syepol_pmz(:,2,:)).^2)./abs(syNpar).*exp(-((1-n_ec./syres)./sybetath./syNpar).^2/2);
0780 %     %
0781 %     syalpha(isnan(syalpha)) = 0;
0782 %     syds = [zeros(1,ny);diff(sys)];
0783 %     syds(isnan(syds)) = 0;
0784 %     %
0785 %     syPinc = ones(ns,1)*yP0.*exp(-cumsum(syalpha.*syds));
0786 %     sydPincds = syPinc.*syalpha;
0787 %     %
0788 %     if exist('wave_s','var') == 1,
0789 %         %
0790 %         sywpe_s = sqrt(syne_s*qe^2/e0/me);
0791 %         syabsphi_s = squeeze(sqrt(abs(syphi_pmz_s(:,1,:)).^2 + abs(syphi_pmz_s(:,2,:)).^2 + abs(syphi_pmz_s(:,3,:)).^2));
0792 %         %
0793 %         syalpha_s = 2^(1-n_ec)*n_ec/sqrt(2*pi)/fact_dke_jd(n_ec - 1).*sywpe_s.^2.*syres_s.^(2*(n_ec - 1))...
0794 %             .*sybetath_s.^(2*n_ec - 3)*pi/clum/omega_rf./syabsphi_s.*syNperp_s.^(2*(n_ec - 1)).*squeeze(abs(syepol_pmz_s(:,2,:)).^2)./abs(syNpar_s).*exp(-((1-n_ec./syres_s)./sybetath_s./syNpar_s).^2/2);
0795 %         %
0796 %         syPinc_s = ones(ns_s,1)*yP0_s.*exp(-cumsum(syalpha_s.*[zeros(1,ny);diff(sys_s)]));
0797 %         sydPincds_s = syPinc_s.*syalpha_s;
0798 %         %
0799 %     end
0800 %     %
0801 % end
0802 % %
0803 % % frequencies
0804 % %
0805 % if (opt(2) >=1 | opt(3) >= 1) & ny == 1% only one ray
0806 %     %
0807 %     syrho = interp1(ppsin,prho,sypsin,method);
0808 %     %
0809 %     omega_rf_ghz = omega_rf/(2*pi*1e9);
0810 %     somega_ce_ghz = syomega_ce/(2*pi*1e9);
0811 %     %
0812 %     Nsomega_ce_ghz = somega_ce_ghz*(1:N);
0813 %     NsNpar = repmat(syNpar,[1,N]);
0814 %     NsdNpar = repmat(sydNpar,[1,N]);
0815 %     Nsybetath = repmat(sybetath,[1,N]);
0816 %     %
0817 %     Nsomega_ce_ghz_m = Nsomega_ce_ghz - (abs(NsNpar) - NsdNpar)*vparn.*Nsybetath*omega_rf_ghz;
0818 %     Nsomega_ce_ghz_p = Nsomega_ce_ghz + (abs(NsNpar) - NsdNpar)*vparn.*Nsybetath*omega_rf_ghz;
0819 %     %
0820 %     [imask_res,Nmask_res] = find(omega_rf_ghz > Nsomega_ce_ghz_m & omega_rf_ghz < Nsomega_ce_ghz_p);
0821 %     %
0822 %     for iN = 1:N
0823 %         mask_N = find(Nmask_res == iN);
0824 %         if ~isempty(mask_N),
0825 %             rhomin_res(iN) = min(syrho(imask_res(mask_N)));
0826 %             rhomax_res(iN) = max(syrho(imask_res(mask_N)));
0827 %             disp(['Harmonic: ',num2str(iN),', rho_min = ',num2str(rhomin_res(iN)),', rho_max = ',num2str(rhomax_res(iN))])
0828 %         else
0829 %             rhomin_res(iN) = NaN;
0830 %             rhomax_res(iN) = NaN;
0831 %             disp(['No resonance at harmonic: ',num2str(iN)])
0832 %         end
0833 %     end
0834 %     disp(['===> global grid: rho_min = ',num2str(min(rhomin_res)),', rho_max = ',num2str(max(rhomax_res))])
0835 % end
0836 % %
0837 % % beam size - gaussian width
0838 % %
0839 % % Note, we assume that the first ray is the central ray
0840 % %
0841 % if opt(4) >= 1,
0842 %     syispsin = NaN*zeros(ns,ny);
0843 %     %
0844 %     scsign = [-1;sign(diff(sypsin(:,1)))];
0845 %     sysign = [ones(1,ny);sign(diff(sypsin))];
0846 %     %
0847 %     for is=1:ns,
0848 %         scpsin = sypsin(is,1);%central ray
0849 %         %
0850 %         sysignmask = scsign(is) ~= sysign;
0851 %         %
0852 %         [isy,iyy] = find(abs(sypsin - scpsin) == repmat(min(abs(sypsin - scpsin) + sysignmask),[ns,1]));
0853 %         for iy = 1:length(isy),
0854 %             syispsin(is,iyy(iy)) = isy(iy);
0855 %         end
0856 %     end
0857 %     %
0858 %     sycR = repmat(syR(:,1),[1,ny]);
0859 %     sycZ = repmat(syZ(:,1),[1,ny]);
0860 %     sycphi = repmat(syphi(:,1),[1,ny]);
0861 %     %
0862 %     scBR = syBR(:,1);
0863 %     scBZ = syBZ(:,1);
0864 %     scBPHI = syBPHI(:,1);
0865 %     scB = syB(:,1);
0866 %     %
0867 %     % unit vector in the field direction
0868 %     %
0869 %     sycbR = repmat(scBR./scB,[1,ny]);
0870 %     sycbZ = repmat(scBZ./scB,[1,ny]);
0871 %     sycbPHI = repmat(scBPHI./scB,[1,ny]);
0872 %     %
0873 %     mask = isnan(syispsin);
0874 %     syispsin(mask) = 1;
0875 %     %
0876 %     syispsin = syispsin + repmat((0:ny-1)*ns,[ns,1]);
0877 %     %
0878 %     sydpar = (syR(syispsin) - sycR).*sycbR + (syZ(syispsin) - sycZ).*sycbZ + (syphi(syispsin) - sycphi).*sycR.*sycbPHI;
0879 %     %
0880 %     sydpar(mask) = NaN;
0881 %     %
0882 %     sdpar = (max(sydpar') - min(sydpar'))/2;
0883 %     sdNpar = 2*clum/omega_rf./sdpar;
0884 % end
0885 %
0886 % Display
0887 %
0888 style = '-';
0889 style2 = 'none';
0890 style3 = '--';
0891 marker = 'none';
0892 marker1 = '+';
0893 marker2 = '.';
0894 marker3 = 'o';
0895 color1 = 'r';
0896 color2 = 'b';
0897 color3 = 'm';
0898 color4 = 'g';
0899 color5 = 'k';
0900 %
0901 colors = {color1,color2,color3,color4};
0902 ncols = length(colors);
0903 %
0904 width1 = 2;
0905 width2 = 0.5+16i;
0906 %
0907 if ~isfield(opt,'font'),
0908     opt.font = 20+14i;
0909 end
0910 %
0911 if ~isfield(opt,'axs'),
0912     figure(1),clf,set(1,'name','Poloidal trajectories')
0913     ax = gca;
0914     red = 0.9;
0915     lspace = 0.7;
0916     bspace = 0.7;
0917 else
0918     ax = opt.axs(1);
0919     red = 1;
0920     lspace = NaN;
0921     bspace = NaN;
0922 end   
0923 % Add vessel and tiles if TCV
0924 if any(strfind(equil.id,'TCV')) && strcmp(opt.vessel_opt,'y')
0925     tcvview('avt'); % Plot axes, vessel and tiles
0926     tcvview('color','t',[0.75 0.75 0.75]); % Tiles in light gray
0927     title(''); % Remove title
0928     hold on;
0929     Rpdisp = Rp; % Centering of flux surfaces and rays on the magnetic axis
0930     Zpdisp = Zp;
0931     %
0932     legpol = NaN;
0933 else
0934     Rpdisp = 0;
0935     Zpdisp = 0;
0936     %
0937 end
0938 %
0939 xlim = [xmin + Rpdisp,xmax + Rpdisp];
0940 ylim = [ymin + Zpdisp,ymax + Zpdisp];
0941 xlab = 'x (m)';
0942 ylab = 'y (m)';
0943 tit = '';
0944 %
0945 for iy = 1:ny,
0946     %
0947     if opt.luke,
0948         %
0949         mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
0950         mask_abs =rays{iy}.mask_abs;
0951         %
0952         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0953             mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
0954             mask_abs =rays{iy}.mask_abs_coll_luke;
0955         end
0956         %
0957     else
0958         mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
0959         mask_abs = rays{iy}.mask_abs_lin;
0960         %
0961         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0962             mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max;
0963             mask_abs =rays{iy}.mask_abs_lin_coll;            
0964         end
0965         %
0966     end
0967     %
0968     for ipass = 1:length(rays{iy}.mask_pass),
0969         %
0970         if isfield(opt,'mode') && opt.mode == 'y',
0971             %
0972             if opt.luke,
0973                 mask_p = rays{iy}.mask_p;
0974                 mask_m = rays{iy}.mask_m;
0975                 %
0976                 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0977                      mask_p = rays{iy}.mask_p_coll_luke;
0978                      mask_m = rays{iy}.mask_m_coll_luke;
0979                 end
0980             else
0981                 mask_p = rays{iy}.mask_p_lin;
0982                 mask_m = rays{iy}.mask_m_lin;
0983                 %
0984                 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
0985                     mask_p = rays{iy}.mask_p_lin_coll;
0986                     mask_m = rays{iy}.mask_m_lin_coll;
0987                 end
0988             end
0989             %
0990             graph1D_jd(rays{iy}.sx(mask_p{1} & rays{iy}.mask_pass{ipass}) + Rpdisp,rays{iy}.sy(mask_p{1} & rays{iy}.mask_pass{ipass}) + Zpdisp,0,0,'','','',NaN,xlim,ylim,style,marker,rays{iy}.color,width1,opt.font,ax);
0991             graph1D_jd(rays{iy}.sx(mask_m{1} & rays{iy}.mask_pass{ipass}) + Rpdisp,rays{iy}.sy(mask_m{1} & rays{iy}.mask_pass{ipass}) + Zpdisp,0,0,'','','',legpol,xlim,ylim,style3,marker,rays{iy}.color,width1,opt.font,ax);
0992             for ip = 2:length(mask_p),
0993                 graph1D_jd(rays{iy}.sx(mask_p{ip} & rays{iy}.mask_pass{ipass}) + Rpdisp,rays{iy}.sy(mask_p{ip} & rays{iy}.mask_pass{ipass}) + Zpdisp,0,0,'','','',NaN,xlim,ylim,style,marker,rays{iy}.color,width1,opt.font,ax);
0994             end
0995             for im = 2:length(mask_m),
0996                 graph1D_jd(rays{iy}.sx(mask_m{im} & rays{iy}.mask_pass{ipass}) + Rpdisp,rays{iy}.sy(mask_m{im} & rays{iy}.mask_pass{ipass}) + Zpdisp,0,0,'','','',NaN,xlim,ylim,style3,marker,rays{iy}.color,width1,opt.font,ax);
0997             end
0998         else
0999             graph1D_jd(rays{iy}.sx(mask_max & rays{iy}.mask_pass{ipass}) + Rpdisp,rays{iy}.sy(mask_max & rays{iy}.mask_pass{ipass}) + Zpdisp,0,0,'','','',NaN,xlim,ylim,rays{iy}.style,marker,rays{iy}.color,width1,opt.font,ax);
1000         end  
1001     end
1002     %
1003     graph1D_jd(rays{iy}.sx(mask_max & mask_abs) + Rpdisp,rays{iy}.sy(mask_max & mask_abs) + Zpdisp,0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font,ax);
1004     %
1005     if ~isempty(rays{iy}.refs),
1006         for iref = 1:length(rays{iy}.refs),
1007             graph1D_jd(rays{iy}.refs{iref}.x + Rpdisp,rays{iy}.refs{iref}.y + Zpdisp,0,0,'','','',NaN,xlim,ylim,style,marker3,rays{iy}.color,width2,opt.font,ax);
1008         end
1009     end 
1010 end
1011 [ax] = graph1D_jd(ptx_disp' + Rpdisp,pty_disp' + Zpdisp,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color5,width2,opt.font,ax,red,lspace,bspace);
1012 if 0,%opt.res >= 1,
1013     [ax] = graph2D_jd(R + Rpdisp,Z + Zpdisp,yres','','','',xlim,ylim,0,NaN,[1:10],0,style,color5,width1,opt.font,ax);
1014 end
1015 axis(ax,'equal');
1016 axis(ax,[xlim,ylim]);
1017 if any(strfind(equil.id,'TCV')) && strcmp(opt.vessel_opt,'y')
1018    axis(ax,'auto');
1019    set(ax,'XTick',[0:0.2:1.4],'YTick',[-0.8:0.2:0.8]);
1020 end
1021 %
1022 figname = [waves{1}.equil_id,'-'];
1023 for ii = 1:length(waves)
1024     figname = [figname,waves{ii}.id,'-'];
1025 end
1026 print_jd(opt.print,[figname,'Pol_traj'],'figures_WAVE');
1027 %
1028 if isfield(opt,'pol') && opt.pol == 1,
1029     return
1030 end
1031 %
1032 if ~isfield(opt,'axs'),
1033     figure(2),clf,set(2,'name','Toroidal trajectories')
1034     ax = gca;
1035 else
1036     ax = opt.axs(2);
1037 end   
1038 %
1039 xlim = [Xmin,Xmax];
1040 ylim = [Ymin,Ymax];
1041 xlab = 'X (m)';
1042 ylab = 'Y (m)';
1043 tit = '';
1044 %
1045 for iy = 1:ny,
1046     %
1047     if opt.luke,
1048         %
1049         mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1050         mask_abs =rays{iy}.mask_abs;
1051         %
1052         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1053             mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1054             mask_abs =rays{iy}.mask_abs_coll_luke;
1055         end
1056         %
1057     else
1058         mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1059         mask_abs = rays{iy}.mask_abs_lin;
1060         %
1061         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1062             mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max;
1063             mask_abs =rays{iy}.mask_abs_lin_coll;            
1064         end
1065         %
1066     end
1067     %
1068     if isfield(opt,'mode') && opt.mode == 'y',
1069         %
1070         if opt.luke,
1071             mask_p = rays{iy}.mask_p;
1072             mask_m = rays{iy}.mask_m;
1073             %
1074             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1075                  mask_p = rays{iy}.mask_p_coll_luke;
1076                  mask_m = rays{iy}.mask_m_coll_luke;
1077             end
1078         else
1079             mask_p = rays{iy}.mask_p_lin;
1080             mask_m = rays{iy}.mask_m_lin;
1081             %
1082             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1083                 mask_p = rays{iy}.mask_p_lin_coll;
1084                 mask_m = rays{iy}.mask_m_lin_coll;
1085             end
1086         end
1087         %
1088         graph1D_jd(rays{iy}.sX(mask_p{1}),rays{iy}.sY(mask_p{1}),0,0,'','','',NaN,xlim,ylim,style,marker,rays{iy}.color,width1,opt.font,ax);
1089         graph1D_jd(rays{iy}.sX(mask_m{1}),rays{iy}.sY(mask_m{1}),0,0,'','','',legpol,xlim,ylim,style3,marker,rays{iy}.color,width1,opt.font,ax);
1090         for ip = 2:length(mask_p),
1091             graph1D_jd(rays{iy}.sX(mask_p{ip}),rays{iy}.sY(mask_p{ip}),0,0,'','','',NaN,xlim,ylim,style,marker,rays{iy}.color,width1,opt.font,ax);
1092         end
1093         for im = 2:length(mask_m),
1094             graph1D_jd(rays{iy}.sX(mask_m{im}),rays{iy}.sY(mask_m{im}),0,0,'','','',NaN,xlim,ylim,style3,marker,rays{iy}.color,width1,opt.font,ax);
1095         end
1096     else
1097         graph1D_jd(rays{iy}.sX(mask_max),rays{iy}.sY(mask_max),0,0,'','','',NaN,xlim,ylim,style,marker,rays{iy}.color,width1,opt.font,ax);
1098     end
1099     %
1100     graph1D_jd(rays{iy}.sX(mask_max & mask_abs),rays{iy}.sY(mask_max & mask_abs),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font,ax);
1101     if ~isempty(rays{iy}.refs),
1102         for iref = 1:length(rays{iy}.refs),
1103             graph1D_jd(rays{iy}.refs{iref}.X,rays{iy}.refs{iref}.Y,0,0,'','','',NaN,xlim,ylim,style,marker3,rays{iy}.color,width2,opt.font,ax);
1104         end
1105     end 
1106 end
1107 [ax] = graph1D_jd(X_disp',Y_disp',0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color5,width2,opt.font,ax,red,lspace,bspace);
1108 axis(ax,'equal');
1109 axis(ax,[xlim,ylim]);
1110 %
1111 print_jd(opt.print,[figname,'Tor_traj'],'figures_WAVE');
1112 %
1113 smax = 0;
1114 %Nmax = 0;
1115 Spmax = 0;
1116 for iy = 1:ny,
1117     smax = max(smax,rays{iy}.ss(end));
1118     %Nmax = max([Nmax,abs(real(rays{iy}.sNpar))]);
1119     Spmax = max([Spmax,abs(rays{iy}.sphi_xyz(1,:))]);
1120 end
1121 %
1122 if ~isfield(opt,'axs'),
1123     figure(3),clf,set(3,'name','Radius of deposition peak')
1124     xlim = [0,ny+1];
1125     ylim = [0,1];
1126     xlab = 'ray #';
1127     ylab = '\rho^{abs}';
1128     tit = '';
1129     %
1130     if opt.luke,
1131         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1132             graph1D_jd(1:ny,yrho_abs_coll_luke,0,0,'','','',NaN,xlim,ylim,style2,marker3,color1,width1);
1133             leg = {'LUKE (RLA + NRCA)','C3PO (RLA + NRCA)'};
1134         else
1135             graph1D_jd(1:ny,yrho_abs,0,0,'','','',NaN,xlim,ylim,style2,marker3,color1,width1);
1136             leg = {'LUKE (RLA)','C3PO (RLA)'};
1137         end
1138     else
1139         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,         
1140             leg = {'C3PO (RLA + NRCA)'};
1141         else
1142             leg = {'C3PO (RLA)'};
1143         end
1144     end
1145     %
1146     if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1147         graph1D_jd(1:ny,yrho_abs_lin_coll,0,0,xlab,ylab,tit,leg,xlim,ylim,style2,marker1,color2,width1,opt.font,gca,red,lspace,bspace);
1148     else        
1149         graph1D_jd(1:ny,yrho_abs_lin,0,0,xlab,ylab,tit,leg,xlim,ylim,style2,marker1,color2,width1,opt.font,gca,red,lspace,bspace);
1150     end
1151     %
1152     print_jd(opt.print,[figname,'Rad_dep_peak'],'figures_WAVE');
1153     %
1154     figure(4),clf,set(4,'name','n// of deposition peak')
1155     xlim = [0,ny+1];
1156     ylim = NaN;%[-Nmax,Nmax];
1157     xlab = 'ray #';
1158     ylab = 'n_{||}^{abs}';
1159     tit = '';
1160     %
1161     if opt.luke,
1162         if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1163             graph1D_jd(1:ny,real(yNpar_abs_coll_luke),0,0,'','','',NaN,xlim,ylim,style2,marker3,color1,width1);
1164         else
1165             graph1D_jd(1:ny,real(yNpar_abs),0,0,'','','',NaN,xlim,ylim,style2,marker3,color1,width1);
1166         end
1167     end
1168     %
1169     if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1170         graph1D_jd(1:ny,real(yNpar_abs_lin_coll),0,0,xlab,ylab,tit,leg,xlim,ylim,style2,marker1,color2,width1,opt.font,gca,red,lspace,bspace);
1171     else
1172         graph1D_jd(1:ny,real(yNpar_abs_lin),0,0,xlab,ylab,tit,leg,xlim,ylim,style2,marker1,color2,width1,opt.font,gca,red,lspace,bspace);
1173     end
1174     %
1175     print_jd(opt.print,[figname,'Npar_dep_peak'],'figures_WAVE');
1176 end
1177 %
1178 if ~isfield(opt,'axs'),
1179     figure(5),clf,set(5,'name','Total deposition profile')
1180     ax = gca;
1181     ylim = NaN;
1182 else
1183     ax = opt.axs(3);
1184     leg = NaN;
1185     if sum(dPdrho_lin_coll) ~= 0,
1186         ylim = [0,max(dPdrho_lin_coll/1e6)]*1.2;
1187     else
1188         ylim = [0,max(dPdrho_lin/1e6)]*1.2;
1189     end
1190 end   
1191 %
1192 xlim = [0,1];
1193 xlab = plab;
1194 ylab = 'dPd\rho (MW)';
1195 tit = '';
1196 %
1197 if opt.luke,
1198     if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1199         graph1D_jd(rho_dep,dPdrho_coll_luke/1e6,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width1,opt.font,ax);
1200     else
1201         graph1D_jd(rho_dep,dPdrho/1e6,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width1,opt.font,ax);
1202     end
1203 end
1204 %
1205 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1206     graph1D_jd(rho_dep,dPdrho_lin_coll/1e6,0,0,xlab,ylab,tit,leg,xlim,ylim,style3,marker1,color2,width1,opt.font,ax,red,lspace,bspace); 
1207 else
1208     graph1D_jd(rho_dep,dPdrho_lin/1e6,0,0,xlab,ylab,tit,leg,xlim,ylim,style3,marker1,color2,width1,opt.font,ax,red,lspace,bspace);
1209 end
1210 %
1211 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1 && opt.luke,
1212     if ~isfield(opt,'axs'),
1213         figure(6),clf,set(6,'name','Total deposition profile (quasilinear, LUKE)')
1214         ax = gca;
1215         ylim = NaN;
1216     else
1217         ax = opt.axs(3);
1218         leg = NaN;
1219         if sum(dPdrho_coll_luke) ~= 0,
1220             ylim = [0,max(dPdrho_coll_luke/1e6)]*1.2;
1221         end
1222     end   
1223     %
1224     xlim = [0,1];
1225     xlab = plab;
1226     ylab = 'dPd\rho (MW)';
1227     tit = '';
1228     leg = {'LUKE (RLA + NRCA)','LUKE (RLA)'};
1229     %
1230     graph1D_jd(rho_dep,dPdrho_coll_luke/1e6,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width1,opt.font,ax);
1231     graph1D_jd(rho_dep,dPdrho/1e6,0,0,xlab,ylab,tit,leg,xlim,ylim,style3,marker3,color2,width1,opt.font,ax,red,lspace,bspace);
1232     %
1233 else
1234     if ishandle(figure(6))
1235         close(figure(6))
1236     end
1237 end
1238 %
1239 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1240     if ~isfield(opt,'axs'),
1241         figure(7),clf,set(7,'name','Total deposition profile (linear, C3PO)')
1242         ax = gca;
1243         ylim = NaN;
1244     else
1245         ax = opt.axs(3);
1246         leg = NaN;
1247         if sum(dPdrho_lin_coll) ~= 0,
1248             ylim = [0,max(dPdrho_lin_coll/1e6)]*1.2;
1249         end
1250     end   
1251     %
1252     xlim = [0,1];
1253     xlab = plab;
1254     ylab = 'dPd\rho (MW)';
1255     tit = '';
1256     leg = {'C3PO (RLA + NRCA)','C3PO (RLA)'};
1257     %
1258     graph1D_jd(rho_dep,dPdrho_lin_coll/1e6,0,0,'','','',NaN,xlim,ylim,style,marker3,color1,width1,opt.font,ax);
1259     graph1D_jd(rho_dep,dPdrho_lin/1e6,0,0,xlab,ylab,tit,leg,xlim,ylim,style3,marker3,color2,width1,opt.font,ax,red,lspace,bspace);
1260     %
1261 else
1262     if ishandle(figure(7))
1263         close(figure(7))
1264     end
1265 end
1266 %
1267 set(ax,'xtickmode','auto');drawnow;
1268 set(ax,'ytickmode','auto');drawnow;
1269 %
1270 pause(0.1);
1271 %
1272 print_jd(opt.print,[figname,'Tot_dep_prof'],'figures_WAVE');
1273 %
1274 % Ray specific figures
1275 %
1276 disp(' ')
1277 disp('-----> Display specific ray information')
1278 disp(' ')
1279 %
1280 iw = NaN;
1281 ibb = NaN;
1282 ib = 0;
1283 iwdef = 1;
1284 %
1285 if isfield(opt,'nospec') && opt.nospec == 1,
1286     iw = 0;% no specific ray display
1287 end
1288 %
1289 while iw ~= 0 && any(ibb ~= 0)
1290     %
1291     ib = ib + 1;
1292     %
1293     iw = input_dke_yp('wave index (0 to return)',iwdef,0:nw,'',[1,1]);
1294     %
1295     if iw == 0,
1296         continue
1297     end
1298     %
1299     if (isfield(waves{iw}.rayinit,'type') && strcmp(waves{iw}.rayinit.type,'LH')) || (isfield(waves{iw}.rayinit.launch,'type') && strcmp(waves{iw}.rayinit.launch.type,'LH'))
1300         waves{iw}.rayinit.Rp = Rp;
1301         plot_spectrum_from_launch_jd(waves{iw}.rayinit,NaN,-1,[])
1302         %
1303         plot_spectrum_from_launch_jd(waves{iw}.rayinit.launch,NaN,-1,[])
1304         %
1305     end
1306     %
1307     disp(['-----------------------------------------------------------'])
1308     disp(['Initial Z (m) : ',num2str(waves{iw}.rayinit.launch.rZ0)])
1309     disp(['Initial N_par0 : ',num2str(waves{iw}.rayinit.launch.bNpar0)])
1310     disp(['-----------------------------------------------------------'])
1311     disp(['-----------------------------------------------------------'])
1312     %
1313     disp(['Index table for each ray. Missing rays have not propagated in the plasma. '])
1314     disp(['-----------------------------------------------------------'])
1315     for iray = 1:length(waves{iw}.rayinit.launch.bNpar0)
1316         disp_ray_launch = ['Npar0 = ',num2str(waves{iw}.rayinit.launch.bNpar0(iray))];
1317         for iz = 1:length(waves{iw}.rayinit.launch.rZ0)
1318             for ibb0 = 1:nymask(iw),
1319                 iy0 = sum(nymask(1:iw - 1)) + ibb0;
1320                 if abs(waves{iw}.rayinit.launch.bNpar0(iray) - real(rays{iy0}.sNpar(1))) < 0.01,
1321                     if abs(waves{iw}.rayinit.launch.rZ0(iz) - rays{iy0}.sy(1)) < 0.01, 
1322                         disp_ray_launch = [disp_ray_launch,', (#',int2str(ibb0),') ZO = ',num2str(waves{iw}.rayinit.launch.rZ0(iz))];
1323                     end
1324                 end
1325             end
1326         end
1327         disp(disp_ray_launch);
1328     end
1329     disp(['-----------------------------------------------------------'])
1330     %
1331     ibb = input_dke_yp('ray indices (0 to return)',1,0:nymask(iw));
1332     %
1333     if ibb == 0,
1334         continue
1335     end
1336     %
1337     iwdef = 0;
1338     %
1339     iy = sum(nymask(1:iw - 1)) + ibb;
1340     %
1341     rays_disp = rays(iy);
1342     nray = length(rays_disp);
1343     %
1344     f_opt = input_dke_yp('single figure with subplots (0) or multiple figures (1)',0,0:1,'',[1,1]);
1345     %
1346     pow_opt = input_dke_yp('Normalized power deposition (0) or absolute power deposition (1)',0,0:1,'',[1,1]);
1347     if opt.luke,
1348         c_opt = input_dke_yp('Also plot C3PO (1) or not (0)',1,0:1,'',[1,1]);
1349     end
1350     %
1351     figlab = ['wave #',num2str(iw),'; ray #',num2str(ibb)];
1352     %
1353     if waves{iw}.opt_TeNpar == 1,
1354         %
1355         NparsTe = input_dke_yp('Value for n//*sqrt(Te)',6.5);
1356         %
1357     end
1358     %
1359     % poloidal ray trajectory
1360     %
1361     setfig(f_opt,ib,1,figlab)
1362     %
1363     xlim = [xmin,xmax];
1364     ylim = [ymin,ymax];
1365     xlab = 'x (m)';
1366     ylab = 'y (m)';
1367     tit = 'poloidal ray trajectory';
1368     %
1369     for iray  = 1:nray,
1370         %
1371         ray = rays_disp{iray};
1372         %
1373         if opt.luke,
1374             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1375                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1376                 mask_abs = ray.mask_abs_coll_luke;
1377                 mask_p = ray.mask_p_coll_luke;
1378                 mask_m = ray.mask_m_coll_luke;
1379             else
1380                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1381                 mask_abs = ray.mask_abs;
1382                 mask_p = ray.mask_p;
1383                 mask_m = ray.mask_m;
1384             end
1385         else
1386             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1387                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max_lin;
1388                 mask_abs = ray.mask_abs_lin_coll;
1389                 mask_p = ray.mask_p_lin_coll;
1390                 mask_m = ray.mask_m_lin_coll;
1391             else
1392                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1393                 mask_abs = ray.mask_abs_lin;
1394                 mask_p = ray.mask_p_lin;
1395                 mask_m = ray.mask_m_lin;
1396             end
1397         end
1398         %
1399         if isfield(opt,'mode') && opt.mode == 'y',
1400             %
1401             graph1D_jd(ray.sx(mask_p{1}),ray.sy(mask_p{1}),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1402             graph1D_jd(ray.sx(mask_m{1}),ray.sy(mask_m{1}),0,0,'','','',legpol,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1403             for ip = 2:length(mask_p),
1404                 graph1D_jd(ray.sx(mask_p{ip}),ray.sy(mask_p{ip}),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1405             end
1406             for im = 2:length(mask_m),
1407                 graph1D_jd(ray.sx(mask_m{im}),ray.sy(mask_m{im}),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1408             end
1409         else
1410             graph1D_jd(ray.sx(mask_max),ray.sy(mask_max),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1411         end       
1412         %
1413         graph1D_jd(ray.sx(mask_max & mask_abs),ray.sy(mask_max & mask_abs),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1414         %
1415         if ~isempty(ray.refs),
1416             for iref = 1:length(ray.refs),
1417                 graph1D_jd(ray.refs{iref}.x,ray.refs{iref}.y,0,0,'','','',NaN,xlim,ylim,style,marker3,color5,width2,opt.font);
1418             end
1419         end 
1420         %
1421     end
1422     %
1423     graph1D_jd(ptx_disp',pty_disp',0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color5,width2,opt.font);
1424     %
1425     axis('equal');
1426     axis([xlim,ylim]);
1427     %
1428     % toroidal ray trajectory
1429     %
1430     setfig(f_opt,ib,2,figlab)
1431     %
1432     xlim = [Xmin,Xmax];
1433     ylim = [Ymin,Ymax];
1434     xlab = 'X (m)';
1435     ylab = 'Y (m)';
1436     tit = 'toroidal ray trajectory';
1437     %
1438     for iray  = 1:nray,
1439         %
1440         ray = rays_disp{iray};
1441         %
1442         if opt.luke,
1443             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1444                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1445                 mask_abs = ray.mask_abs_coll_luke;
1446                 mask_p = ray.mask_p_coll_luke;
1447                 mask_m = ray.mask_m_coll_luke;
1448             else
1449                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1450                 mask_abs = ray.mask_abs;
1451                 mask_p = ray.mask_p;
1452                 mask_m = ray.mask_m;
1453             end
1454         else
1455             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1456                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max_lin;
1457                 mask_abs = ray.mask_abs_lin_coll;
1458                 mask_p = ray.mask_p_lin_coll;
1459                 mask_m = ray.mask_m_lin_coll;
1460             else
1461                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1462                 mask_abs = ray.mask_abs_lin;
1463                 mask_p = ray.mask_p_lin;
1464                 mask_m = ray.mask_m_lin;
1465             end
1466                 end
1467         %
1468         if isfield(opt,'mode') && opt.mode == 'y',
1469             graph1D_jd(ray.sX(mask_p{1}),ray.sY(mask_p{1}),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1470             graph1D_jd(ray.sX(mask_m{1}),ray.sY(mask_m{1}),0,0,'','','',legpol,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1471             for ip = 2:length(mask_p),
1472                 graph1D_jd(ray.sX(mask_p{ip}),ray.sY(mask_p{ip}),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1473             end
1474             for im = 2:length(mask_m),
1475                 graph1D_jd(ray.sX(mask_m{im}),ray.sY(mask_m{im}),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1476             end
1477         else
1478             graph1D_jd(ray.sX(mask_max),ray.sY(mask_max),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1479         end 
1480         %
1481         graph1D_jd(ray.sX(mask_max & mask_abs),ray.sY(mask_max & mask_abs),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1482         %
1483         if ~isempty(ray.refs),
1484             for iref = 1:length(ray.refs),
1485                 graph1D_jd(ray.refs{iref}.X,ray.refs{iref}.Y,0,0,'','','',NaN,xlim,ylim,style,marker3,color5,width2,opt.font);
1486             end
1487         end 
1488         %
1489     end
1490     %
1491     graph1D_jd(X_disp',Y_disp',0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color5,width2,opt.font);
1492     %
1493     axis('equal');
1494     axis([xlim,ylim]);
1495     %
1496     % parallel index of refraction
1497     %
1498     setfig(f_opt,ib,3,figlab)
1499     %
1500     xlim = [0,max(ray.ss(mask_max))];
1501     ylim = [0,ceil(1.5*max(abs(real(ray.sNpar(mask_max)))))];
1502     xlab = 's (m)';
1503     ylab = 'n_{||}';
1504     tit = '';
1505     %
1506     for iray  = 1:nray,
1507         %
1508         ray = rays_disp{iray};
1509         %
1510         if opt.luke,
1511             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1512                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1513                 mask_abs = ray.mask_abs_coll_luke;
1514                 mask_p = ray.mask_p_coll_luke;
1515                 mask_m = ray.mask_m_coll_luke;
1516             else
1517                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1518                 mask_abs = ray.mask_abs;
1519                 mask_p = ray.mask_p;
1520                 mask_m = ray.mask_m;
1521             end
1522         else
1523             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1524                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max_lin;
1525                 mask_abs = ray.mask_abs_lin_coll;
1526                 mask_p = ray.mask_p_lin_coll;
1527                 mask_m = ray.mask_m_lin_coll;
1528             else
1529                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1530                 mask_abs = ray.mask_abs_lin;
1531                 mask_p = ray.mask_p_lin;
1532                 mask_m = ray.mask_m_lin;
1533             end
1534         end
1535         %
1536         if isfield(opt,'mode') && opt.mode == 'y',
1537             leg = {'|n_{||}| (SW)','|n_{||}| (FW)'};
1538             graph1D_jd(ray.ss(mask_p{1}),abs(real(ray.sNpar(mask_p{1}))),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1539             graph1D_jd(ray.ss(mask_m{1}),abs(real(ray.sNpar(mask_m{1}))),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1540         else
1541             leg = {'|n_{||}|'};
1542             graph1D_jd(ray.ss(mask_max),abs(real(ray.sNpar(mask_max))),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1543         end 
1544         %
1545         if waves{iw}.opt_TeNpar == 1,
1546             graph1D_jd(ray.ss(mask_max),NparsTe./sqrt(ray.sTe(mask_max)),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width2,opt.font);
1547             leg = [leg,[num2str(NparsTe),'/Te^{1/2}'],'n_{||acc}'];
1548         else
1549             leg = [leg,'n_{||acc}'];
1550         end
1551         %
1552         if f_opt == 1,
1553             graph1D_jd(ray.ss(mask_max),abs(imag(ray.sNpar(mask_max))),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width2,opt.font);
1554         else
1555             graph1D_jd(ray.ss(mask_max),abs(imag(ray.sNpar(mask_max))),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width2,opt.font);
1556         end        
1557         %
1558         if isfield(opt,'mode') && opt.mode == 'y',
1559             for ip = 2:length(mask_p),
1560                 graph1D_jd(ray.ss(mask_p{ip}),abs(real(ray.sNpar(mask_p{ip}))),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1561             end
1562             for im = 2:length(mask_m),
1563                 graph1D_jd(ray.ss(mask_m{im}),abs(real(ray.sNpar(mask_m{im}))),0,0,'','','',NaN,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1564             end
1565         end
1566         %
1567         if iray == 1,
1568             graph1D_jd(ray.ss(mask_max & mask_abs),abs(real(ray.sNpar(mask_max & mask_abs))),0,0,xlab,ylab,tit,leg,xlim,ylim,style2,marker2,color5,width2,opt.font);%,gca,red,lspace,bspace
1569         else
1570             graph1D_jd(ray.ss(mask_max & mask_abs),abs(real(ray.sNpar(mask_max & mask_abs))),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1571         end
1572         %
1573     end
1574     %
1575     % perpendicular index of refraction
1576     %
1577     setfig(f_opt,ib,4,figlab)
1578     %
1579     ylab = 'n_{\perp}';
1580     ylim = NaN;
1581     %
1582     for iray  = 1:nray,
1583         %
1584         ray = rays_disp{iray};
1585         %
1586         if opt.luke,
1587             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1588                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1589                 mask_abs = ray.mask_abs_coll_luke;
1590                 mask_p = ray.mask_p_coll_luke;
1591                 mask_m = ray.mask_m_coll_luke;
1592             else
1593                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1594                 mask_abs = ray.mask_abs;
1595                 mask_p = ray.mask_p;
1596                 mask_m = ray.mask_m;
1597             end
1598         else
1599             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1600                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max_lin;
1601                 mask_abs = ray.mask_abs_lin_coll;
1602                 mask_p = ray.mask_p_lin_coll;
1603                 mask_m = ray.mask_m_lin_coll;
1604             else
1605                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1606                 mask_abs = ray.mask_abs_lin;
1607                 mask_p = ray.mask_p_lin;
1608                 mask_m = ray.mask_m_lin;
1609             end
1610         end
1611         %
1612         graph1D_jd(ray.ss,real(ray.sNperp),0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1613         graph1D_jd(ray.ss(mask_abs),real(ray.sNperp(mask_abs)),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1614         %
1615     end
1616     %
1617     % perpendicular normalized energy flow
1618     %
1619     setfig(f_opt,ib,5,figlab)
1620     %
1621     ylab = 'S_{\perp}';
1622     %
1623     for iray  = 1:nray,
1624         %
1625         ray = rays_disp{iray};
1626         %
1627         if opt.luke,
1628             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1629                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_coll_luke;%1:rays{iy}.is_abs_max;
1630                 mask_abs = ray.mask_abs_coll_luke;
1631                 mask_p = ray.mask_p_coll_luke;
1632                 mask_m = ray.mask_m_coll_luke;
1633             else
1634                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max;%1:rays{iy}.is_abs_max;
1635                 mask_abs = ray.mask_abs;
1636                 mask_p = ray.mask_p;
1637                 mask_m = ray.mask_m;
1638             end
1639         else
1640             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1641                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin_coll;%1:rays{iy}.is_abs_max_lin;
1642                 mask_abs = ray.mask_abs_lin_coll;
1643                 mask_p = ray.mask_p_lin_coll;
1644                 mask_m = ray.mask_m_lin_coll;
1645             else
1646                 mask_max = 1:rays{iy}.ns <= rays{iy}.is_abs_max_lin;%1:rays{iy}.is_abs_max_lin;
1647                 mask_abs = ray.mask_abs_lin;
1648                 mask_p = ray.mask_p_lin;
1649                 mask_m = ray.mask_m_lin;
1650             end
1651         end
1652         %
1653         %sphi_abs = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
1654         sphi_abs = abs(ray.sphi_xyz(1,:));
1655         %
1656         graph1D_jd(ray.ss,sphi_abs,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1657         graph1D_jd(ray.ss(mask_abs),sphi_abs(mask_abs),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1658         %
1659     end
1660     %
1661     % power deposition
1662     %
1663     setfig(f_opt,ib,6,figlab);
1664     %
1665     xlim = [0,1];
1666     xlab = plab;
1667     if pow_opt == 0,
1668         ylab = 'dP/d\rho/P_0';
1669     else
1670         ylab = 'dP/d\rho (MW)';
1671     end
1672     %
1673     for iray  = 1:nray,
1674         %
1675         ray = rays_disp{iray};
1676         %
1677         if pow_opt == 0,
1678             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1679                 dP_drho_disp = ray.dP_coll_luke_normdrho.*drhog_dep./drho_dep;
1680             else
1681                 dP_drho_disp = ray.dP_normdrho.*drhog_dep./drho_dep;
1682             end
1683             %
1684             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1685                 dP_lin_drho_disp = ray.dP_lin_coll_normdrho.*drhog_dep./drho_dep;
1686             else
1687                 dP_lin_drho_disp = ray.dP_lin_normdrho.*drhog_dep./drho_dep;
1688             end
1689         else
1690             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1691                 dP_drho_disp = ray.dP_coll_luke_normdrho*ray.P0_2piRp*2*pi*Rp.*drhog_dep./drho_dep/1e6;
1692             else
1693                 dP_drho_disp = ray.dP_normdrho*ray.P0_2piRp*2*pi*Rp.*drhog_dep./drho_dep/1e6;
1694             end
1695             %
1696             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1697                 dP_lin_drho_disp = ray.dP_lin_coll_normdrho*ray.P0_2piRp*2*pi*Rp.*drhog_dep./drho_dep/1e6;
1698             else  
1699                 dP_lin_drho_disp = ray.dP_lin_normdrho*ray.P0_2piRp*2*pi*Rp.*drhog_dep./drho_dep/1e6;
1700             end
1701         end        
1702         %
1703         if opt.luke,
1704             graph1D_jd(rho_dep,dP_drho_disp,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1705             if c_opt == 1,
1706                 if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1707                     leg = {'LUKE (RLA + NRCA)','C3PO (RLA + NRCA)'};
1708                 else
1709                     leg = {'LUKE (RLA)','C3PO (RLA)'};
1710                 end
1711             end
1712         else
1713             if isfield(waves{1},'colldamp') && waves{1}.colldamp == 1,
1714                 leg = {'C3PO (RLA + NRCA)'};
1715             else
1716                 leg = {'C3PO (RLA)'};
1717             end
1718         end
1719         if opt.luke == 0 || c_opt == 1,
1720             graph1D_jd(rho_dep,dP_lin_drho_disp,0,0,xlab,ylab,tit,leg,xlim,ylim,style3,marker,colors{rem(iray,ncols+1)},width1,opt.font);
1721         end
1722         %
1723     end
1724     %
1725     % radial location and reflection
1726     %
1727     if f_opt && ~isempty(ray.skrho),
1728         %
1729         setfig(f_opt,ib,7,figlab)
1730         %
1731         xlim = [0,max(ray.ss(mask_max))];
1732         ylim = [-40,40];
1733         xlab = 's (m)';
1734         ylab = '';
1735         tit = '';
1736         leg = {'n_{\rho}','d\rho/ds (\times 1000)'};
1737         %
1738         snrho = ray.skrho*clum/ray.omega_rf;
1739         sdrhods = 1000*gradient_jd(ray.ss,ray.srho);
1740         %
1741         if isfield(opt,'mode') && opt.mode == 'y',
1742             graph1D_jd(ray.ss(mask_p{1}),snrho(mask_p{1}),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width2,opt.font);
1743             graph1D_jd(ray.ss(mask_p{1}),sdrhods(mask_p{1}),0,0,'','','',leg,xlim,ylim,style,marker,color2,width1,opt.font);
1744             for ip = 2:length(mask_p),
1745                 graph1D_jd(ray.ss(mask_p{ip}),snrho(mask_p{ip}),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width2,opt.font);
1746                 graph1D_jd(ray.ss(mask_p{ip}),sdrhods(mask_p{ip}),0,0,'','','',leg,xlim,ylim,style,marker,color2,width1,opt.font);
1747             end
1748             for im = 1:length(mask_m),
1749                 graph1D_jd(ray.ss(mask_m{im}),snrho(mask_m{im}),0,0,'','','',NaN,xlim,ylim,style3,marker,color1,width2,opt.font);
1750                 graph1D_jd(ray.ss(mask_m{im}),sdrhods(mask_m{im}),0,0,'','','',NaN,xlim,ylim,style3,marker,color2,width1,opt.font);
1751             end
1752         else
1753             graph1D_jd(ray.ss,snrho,0,0,'','','',NaN,xlim,ylim,style,marker,color1,width2,opt.font);
1754             graph1D_jd(ray.ss,sdrhods,0,0,'','','',leg,xlim,ylim,style,marker,color2,width1,opt.font);
1755         end
1756         %
1757         graph1D_jd(ray.ss(mask_abs),snrho(mask_abs),0,0,'','','',NaN,xlim,ylim,style2,marker2,color5,width2,opt.font);
1758         graph1D_jd(ray.ss(mask_abs),sdrhods(mask_abs),0,0,xlab,ylab,tit,NaN,xlim,ylim,style2,marker2,color5,width2,opt.font,gca,red,lspace,bspace);
1759         %
1760     end
1761     %
1762     if ~f_opt && opt.print >= 0,
1763         %
1764         figure(2000 + ib - 1);
1765         cpaxlist = copyobj(get(1000 + ib - 1,'children'),2000 + ib - 1);
1766         %
1767         set(2000 + ib - 1,'position',[0,0,600,900]),
1768         set(cpaxlist(1),'units','normalized','position',[0.65,0.1,0.3,0.2])
1769         set(cpaxlist(2),'units','normalized','position',[0.15,0.1,0.3,0.2])
1770         set(cpaxlist(3),'units','normalized','position',[0.65,0.4,0.3,0.2])
1771         set(cpaxlist(4),'units','normalized','position',[0.35,0.525,0.1,0.1])
1772         set(cpaxlist(5),'units','normalized','position',[0.15,0.4,0.3,0.2])
1773         set(cpaxlist(6),'units','normalized','position',[0.65,0.7,0.3,0.2])
1774         set(cpaxlist(7),'units','normalized','position',[0.15,0.7,0.3,0.2])    
1775         %
1776         print_jd(opt.print,[figname,'Wave_',num2str(iw),'_Ray_',num2str(ibb)],'figures_WAVE');    
1777     end
1778     %
1779     if strcmp(waves{iw}.rayinit.type,'LH')
1780         LH_acces_yp(equil,ray,mask_max,iy);
1781     end
1782 end
1783 %
1784 end
1785 %
1786 function setfig(f_opt,ib,ifig,figlab)
1787 %
1788 if f_opt == 1,
1789     %
1790     fignum = 10*(ib - 1) + 9 + ifig;
1791     figure(fignum),clf,set(fignum,'name',[figlab,'; fig ',num2str(ifig)])
1792     %
1793 else
1794     %
1795     fignum = 1000 + ib - 1;
1796     figure(fignum)
1797     if ifig == 1,
1798         clf,
1799         set(gcf,'resize','on'),
1800         set(gcf,'units','normalized'),
1801         %
1802         set(gcf,'position',[0,0,1,1]),
1803         set(fignum,'name',figlab)
1804     end    
1805     %
1806     subplot(2,3,ifig),%cla
1807     %
1808 end
1809 %
1810 end
1811 %
1812 function dydx = gradient_jd(x,y)
1813     %
1814     dydx = zeros(size(y));
1815     %
1816     dydx(2:end-1) = ...
1817         (1./(x(3:end) - x(2:end-1)) - 1./(x(3:end) - x(1:end-2))).*y(3:end) + ...
1818         (1./(x(2:end-1) - x(1:end-2)) - 1./(x(3:end) - x(2:end-1))).*y(2:end-1) + ...
1819         (1./(x(3:end) - x(1:end-2)) - 1./(x(2:end-1) - x(1:end-2))).*y(1:end-2);
1820     %
1821     dydx(1) = (y(2) - y(1))/(x(2) - x(1));
1822     dydx(end) = (y(end) - y(end-1))/(x(end) - x(end-1));
1823 end
1824 
1825 
1826 
1827 
1828 
1829 
1830 
1831 
1832 
1833 
1834 
1835 
1836 
1837 
1838 
1839 % figure(2),clf
1840 % %
1841 % xlim = [-Rmax,Rmax];
1842 % ylim = [-Rmax,Rmax];
1843 % xlab = 'X (m)';
1844 % ylab = 'Y (m)';
1845 % tit = '';
1846 % %
1847 % [ax] = graph1D_jd(syX,syY,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width2,opt.font);
1848 % if exist('launch','var') == 1,
1849 %     [ax] = graph1D_jd(yX_L,yY_L,0,0,'','','',NaN,xlim,ylim,style2,marker2,color3,width2,opt.font);
1850 %     [ax] = graph1D_jd([yX_L;yX_L + ydX_L],[yY_L;yY_L + ydY_L],0,0,'','','',NaN,xlim,ylim,style,marker,color3,width2,opt.font);
1851 % end
1852 % if exist('wave_s','var') == 1,
1853 %     [ax] = graph1D_jd(syX_s,syY_s,0,0,'','','',NaN,xlim,ylim,style,marker,color4,width2,opt.font);
1854 % end
1855 % if opt(2) >= 1,
1856 %     [ax] = graph1D_jd(syX(is_syres),syY(is_syres),0,0,'','','',NaN,xlim,ylim,style2,marker3,color5,width1,opt.font);
1857 % end
1858 % [ax] = graph1D_jd(X_disp',Y_disp',0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width1,opt.font,gca,red,lspace,bspace);
1859 % axis('equal');
1860 % axis([xlim,ylim])
1861 % %
1862 % % ------------------------
1863 % %
1864 % if xvar == 's',
1865 %     xlab = 's (m)';
1866 %     x_disp = sys;
1867 %     if exist('wave_s','var') == 1,
1868 %         x_disp_s = sys_s;
1869 %     end
1870 %     xlim = [0,smax];
1871 % elseif xvar == 'R',
1872 %     xlab = 'R (m)';
1873 %     x_disp = syR;
1874 %     if exist('wave_s','var') == 1,
1875 %         x_disp_s = syR_s;
1876 %     end
1877 %     xlim = [Rmin,Rmax];
1878 % elseif xvar == 'p',
1879 %     xlab = '(\psi/\psi_a)^{1/2}';
1880 %     x_disp = sqrt(sypsin);
1881 %     if exist('wave_s','var') == 1,
1882 %         x_disp_s = sqrt(sypsin_s);
1883 %     end
1884 %     xlim = [0,1];
1885 % elseif xvar == 'r',
1886 %     xlab = '\rho';
1887 %     x_disp = syrho;
1888 %     if exist('wave_s','var') == 1,
1889 %         x_disp_s = syrho_s;
1890 %     end
1891 %     xlim = [0,1];
1892 % end
1893 % %
1894 % figure(3),clf
1895 % %
1896 % ylim = NaN;
1897 % ylab = 'N_{||}';
1898 % tit = '';
1899 % %
1900 % if opt(1) == 1 & ~isempty(sydNpar)
1901 %     [ax] = graph1D_jd([x_disp,x_disp],[syNpar + sydNpar,syNpar - sydNpar],0,0,'','','',NaN,xlim,ylim,style,marker,color2,width2,opt.font);
1902 % end
1903 % if exist('wave_s','var') == 1,
1904 %     [ax] = graph1D_jd(x_disp_s,syNpar_s,0,0,'','','',NaN,xlim,ylim,style,marker,color4,width2,opt.font);
1905 % end
1906 % if opt(2) >= 1,
1907 %     [ax] = graph1D_jd(x_disp(is_syres),syNpar(is_syres),0,0,'','','',NaN,xlim,ylim,style2,marker3,color5,width1,opt.font);
1908 % end
1909 % [ax] = graph1D_jd(x_disp,syNpar,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width1,opt.font,gca,red,lspace,bspace);
1910 % %
1911 % figure(4),clf
1912 % %
1913 % ylim = NaN;
1914 % ylab = 'N_{\perp}';
1915 % tit = '';
1916 % %
1917 % if exist('wave_s','var') == 1,
1918 %     [ax] = graph1D_jd(x_disp_s,syNperp_s,0,0,'','','',NaN,xlim,ylim,style,marker,color4,width2,opt.font);
1919 % end
1920 % if opt(2) >= 1,
1921 %     [ax] = graph1D_jd(x_disp(is_syres),syNperp(is_syres),0,0,'','','',NaN,xlim,ylim,style2,marker3,color5,width1,opt.font);
1922 % end
1923 % [ax] = graph1D_jd(x_disp,syNperp,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width1,opt.font,gca,red,lspace,bspace);
1924 % %
1925 % if exist('launch','var') == 1,
1926 %     figure(5),clf
1927 % %
1928 %     xlim_L = [min([yphi_L,0]),max([yphi_L,0])]*1.1*180/pi;
1929 %     ylim_L = [min([yZ_L,0]),max([yZ_L,0])]*1.1;
1930 %     xlab_L = '\phi_L (deg.)';
1931 %     ylab = 'Z_L (m)';
1932 %     tit = '';
1933 % %
1934 %     [ax] = graph1D_jd(yphi_L*180/pi,yZ_L,0,0,xlab_L,ylab,tit,NaN,xlim_L,ylim_L,style2,marker2,color3,width2,opt.font,gca,red,lspace,bspace);
1935 %     for iy = 1:ny
1936 %         text(yphi_L(iy)*180/pi,yZ_L(iy),['\alpha_L = ',num2str(round(yalpha_L(iy)*180/pi)),'\newline\beta_L = ',num2str(round(ybeta_L(iy)*180/pi))],'fontsize',10);
1937 %     end
1938 % end
1939 % %
1940 % if opt(3) >= 1,
1941 %     figure(6),clf
1942 %     %
1943 %     ylim = [0,1.1];
1944 %     ylab = 'P^b/P_0^b';
1945 %     tit = '';
1946 %     %
1947 %     if exist('wave_s','var') == 1,
1948 %         [ax] = graph1D_jd(x_disp_s,syPinc_s./(ones(ns_s,1)*yP0_s),0,0,'','','',NaN,xlim,ylim,style,marker,color4,width2,opt.font);
1949 %     end
1950 %     if opt(2) >= 1,
1951 %         [ax] = graph1D_jd(x_disp(is_syres),syPinc(is_syres)./missing,0,0,'','','',NaN,xlim,ylim,style2,marker3,color5,width1,opt.font);
1952 %     end
1953 %     [ax] = graph1D_jd(x_disp,syPinc./(ones(ns,1)*yP0),0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width1,opt.font,gca,red,lspace,bspace);
1954 %     %
1955 %     figure(7),clf
1956 %     %
1957 %     ylim = NaN;
1958 %     ylab = 'dP^b/ds/P_0^b (m^{-1})';
1959 %     tit = '';
1960 %     %
1961 %     if exist('wave_s','var') == 1,
1962 %         [ax] = graph1D_jd(x_disp_s,sydPincds_s./(ones(ns_s,1)*yP0_s),0,0,'','','',NaN,xlim,ylim,style,marker,color4,width2,opt.font);
1963 %     end
1964 %     if opt(2) >= 1,
1965 %         [ax] = graph1D_jd(x_disp(is_syres),sydPincds(is_syres)./missing,0,0,'','','',NaN,xlim,ylim,style2,marker3,color5,width1,opt.font);
1966 %     end
1967 %     [ax] = graph1D_jd(x_disp,sydPincds./(ones(ns,1)*yP0),0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width1,opt.font,gca,red,lspace,bspace);
1968 %     %
1969 % end
1970 % %
1971 % if (opt(2) >=1 | opt(3) >= 1) & ny == 1% only one ray
1972 %     %
1973 %     figure(8),clf
1974 %     %
1975 %     ylim = [0,(ceil(omega_rf_ghz/10) + 1)*10];
1976 %     ylab = '\omega/(2\pi) (GHz)';
1977 %     tit = '';
1978 %     str_vparn = [num2str(vparn),'  '];
1979 %     leg = ['n\omega_{ce}                         ';...
1980 %            'n\omega_{ce} \pm ',str_vparn(1:3),'v_{Te}\cdotk_{//}'];
1981 %     %
1982 %     graph1D_jd(x_disp,Nsomega_ce_ghz(:,1),0,0,'','','',NaN,xlim,ylim,style,marker,color2,width1,opt.font,gca,red,lspace,bspace);
1983 %     graph1D_jd(x_disp,Nsomega_ce_ghz_m(:,1),0,0,'','','',NaN,xlim,ylim,style,marker,color3,width1);
1984 %     graph1D_jd(x_disp,Nsomega_ce_ghz(:,2:N),0,0,'','','',NaN,xlim,ylim,style,marker,color2,width1);
1985 %     graph1D_jd(x_disp,Nsomega_ce_ghz_m(:,2:N),0,0,'','','',NaN,xlim,ylim,style,marker,color3,width1);
1986 %     graph1D_jd(x_disp,Nsomega_ce_ghz_p,0,0,'','','',NaN,xlim,ylim,style,marker,color3,width1);
1987 %     graph1D_jd(xlim,[omega_rf_ghz,omega_rf_ghz],0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color1,width1);
1988 %     %
1989 % end
1990 % %
1991 % if opt(4) >=1
1992 %     %
1993 %     figure(9),clf
1994 %     %
1995 %     ylim = NaN;
1996 %     ylab = 'N_{||}';
1997 %     tit = '';
1998 % %
1999 %     [ax] = graph1D_jd(x_disp(:,1),syNpar(:,1),0,0,'','','',NaN,xlim,ylim,style,marker,color2,width1,opt.font);
2000 %     [ax] = graph1D_jd([x_disp(:,1),x_disp(:,1)],[syNpar(:,1)-sdNpar',syNpar(:,1)+sdNpar'],0,0,'','','',NaN,xlim,ylim,style,marker,color3,width1,opt.font);
2001 %     [ax] = graph1D_jd(x_disp,syNpar,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color1,width2,opt.font,gca,red,lspace,bspace);
2002 %     %
2003 % end
2004 
2005