raypostprocessing_yp

PURPOSE ^

C3PO - Post processing of the ray calculation

SYNOPSIS ^

function [iabs_out,rho_ray,theta_ray,x_ray,y_ray,z_ray,Npar_ray,Nperp_ray,epol_pmz_ray,phi_xyz_ray,pmode_ray,D_ray,s_ray,alphaphi_ray,B_ray,BRHO_ray,BP_ray,Bz_ray,Te_ray,ne_ray,zTi_ray,zni_ray,r_ray,salpha_ray,gradrhon_ray,tau_ray,eikval_ray,krho_ray,m_ray,kz_ray,n_ray,wpe2,wce2] = raypostprocessing_yp(tout,yout,yproc,rayparam,equil_fit,fluct_fit,display_text,display_mode,p_opt,filename,subtitle)

DESCRIPTION ^

C3PO - Post processing of the ray calculation

 Post processing of the ray calculation

 INPUT :
    - tout: ray time (w*t) [1,nray)
    - yout: ray coordinate (psin,theta,phi,kpsin,m,n,omega,s) [8,nray] (Nota: krho = krho*ap, rho = (R(:,1)-Rp)/ap
    - yproc: ray coordinate (psin,theta,phi,kpsin,m,n,omega,s) [8,nray] (Nota: krho = krho*ap, rho = (R(:,1)-Rp)/ap
    - rayparam: ray-tracing parameter
    - equil_fit: interpolated magnetic equilibrium structure
     - fluct_fit: interpolated plasma fluctuation structure

 OUTPUT : 

 By Yves Peysson (CEA/IRFM, yves.peysson@cea.fr) and Joan Decker (CEA/IRFM, joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [iabs_out,rho_ray,theta_ray,x_ray,y_ray,z_ray,Npar_ray,Nperp_ray,epol_pmz_ray,phi_xyz_ray,pmode_ray,D_ray,s_ray,alphaphi_ray,B_ray,BRHO_ray,BP_ray,Bz_ray,Te_ray,ne_ray,zTi_ray,zni_ray,r_ray,salpha_ray,gradrhon_ray,tau_ray,eikval_ray,krho_ray,m_ray,kz_ray,n_ray,wpe2,wce2] = raypostprocessing_yp(tout,yout,yproc,rayparam,equil_fit,fluct_fit,display_text,display_mode,p_opt,filename,subtitle)
0002 %C3PO - Post processing of the ray calculation
0003 %
0004 % Post processing of the ray calculation
0005 %
0006 % INPUT :
0007 %    - tout: ray time (w*t) [1,nray)
0008 %    - yout: ray coordinate (psin,theta,phi,kpsin,m,n,omega,s) [8,nray] (Nota: krho = krho*ap, rho = (R(:,1)-Rp)/ap
0009 %    - yproc: ray coordinate (psin,theta,phi,kpsin,m,n,omega,s) [8,nray] (Nota: krho = krho*ap, rho = (R(:,1)-Rp)/ap
0010 %    - rayparam: ray-tracing parameter
0011 %    - equil_fit: interpolated magnetic equilibrium structure
0012 %     - fluct_fit: interpolated plasma fluctuation structure
0013 %
0014 % OUTPUT :
0015 %
0016 % By Yves Peysson (CEA/IRFM, yves.peysson@cea.fr) and Joan Decker (CEA/IRFM, joan.decker@cea.fr)
0017 %
0018 if nargin < 5,
0019     error('Not enough input arguments in raypostprocessing_yp.m');
0020 end
0021 if nargin < 6,
0022     fluct_fit = [];
0023 end
0024 if nargin < 7,
0025     display_text = 0;
0026 end
0027 if nargin < 8,
0028     display_mode = 0;
0029 end
0030 if nargin < 9,
0031     p_opt = -1;
0032 end
0033 if nargin < 10,
0034     filename = '';
0035 end
0036 if nargin < 11,
0037     subtitle = '';
0038 end
0039 if ~isstruct(fluct_fit),
0040     fluct_fit = [];
0041 end
0042 %
0043 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0044 %
0045 tn_ray = tout;
0046 rho_ray = yout(:,1);
0047 theta_ray = yout(:,2);
0048 z_ray = yout(:,3);
0049 krho_ray = yout(:,4); 
0050 m_ray = yout(:,5); 
0051 kz_ray = yout(:,6);
0052 %
0053 if ~isinf(equil_fit.Rp),%for backward compatibility
0054     phi_ray = yout(:,2)/equil_fit.Rp;
0055     n_ray = yout(:,6)*equil_fit.Rp;
0056 else
0057     phi_ray = NaN(size(yout(:,3)));%cylindral magnetic configuration
0058     n_ray = NaN(size(yout(:,6)));%cylindral magnetic configuration
0059 end
0060 %
0061 omega_ray = yout(:,7);
0062 s_ray = yout(:,8)*equil_fit.ap;
0063 %
0064 omega = yout(1,7);
0065 %
0066 if ~isempty(yproc) & size(yproc,2) >= 26,
0067     D_ray = yproc(:,1);
0068     Npar_ray = yproc(:,2);
0069     Nperp_ray = yproc(:,3);
0070     epol_pmz_ray = [yproc(:,4) + i*yproc(:,5),yproc(:,6)+i*yproc(:,7),yproc(:,8)+i*yproc(:,9)];
0071     phi_xyz_ray = [yproc(:,10),yproc(:,11),yproc(:,12)];
0072     alphaphi_ray = yproc(:,13)*omega/clum;
0073     x_ray = yproc(:,14)*equil_fit.ap;
0074     y_ray = yproc(:,15)*equil_fit.ap;
0075     B_ray = yproc(:,16);
0076     BRHO_ray = yproc(:,17).*B_ray;
0077     BP_ray = yproc(:,18).*B_ray;
0078     Bz_ray = yproc(:,19).*B_ray;
0079     r_ray = yproc(:,20)*equil_fit.ap;
0080     salpha_ray = yproc(:,21);
0081     gradrhon_ray = yproc(:,22);
0082     eikval_ray = yproc(:,23);
0083     Upsilonn_ray = yproc(:,24);%(1 + x/Rp)/ap where Upsilon = 1 + x/Rp
0084     Te_ray = yproc(:,25);
0085     ne_ray = yproc(:,26);
0086     %
0087     ns = length(equil_fit.zZi);
0088     %
0089     for is = 1:ns,
0090         zTi_ray(:,is) = yproc(:,26+is);
0091         zni_ray(:,is) = yproc(:,26+ns+is);
0092     end  
0093     %
0094     tau_ray = yproc(:,27+2*ns);
0095     %
0096     for ii=1:length(Npar_ray),
0097         [Nperpp,Nperpm,Kxyz,ep_xyz,em_xyz,ep_pmz,em_pmz,ep_pyk,em_pyk,phip_xyz,phim_xyz,phip_pmz,phim_pmz,phip_pyk,phim_pyk,sigmap,sigmam,Dp,Dm,ps,ys] = colddisp_dke_jd(Npar_ray(ii),omega,ne_ray(ii),zni_ray(ii,:),equil_fit.zZi,equil_fit.zmi,B_ray(ii));
0098         %
0099         pmode_ray(ii,1) = abs(Nperpp - Nperp_ray(ii)) <= abs(Nperpm - Nperp_ray(ii)); 
0100         %
0101         Npar_acc = colddisp_dke_jd(NaN,omega,ne_ray(ii),zni_ray(ii,:),equil_fit.zZi,equil_fit.zmi,B_ray(ii));%accessibility value
0102         %
0103         Npar_ray(ii) = Npar_ray(ii) + i*real(Npar_acc);
0104         %
0105         wpe2(ii) = ps(1)^2;
0106         wce2(ii) = ys(1)^2;
0107     end
0108     %
0109     iabs_out = length(Npar_ray);
0110     %
0111 else
0112     %
0113     [x_ray,y_ray,Npar_ray,Nperp_ray,epol_pmz_ray,phi_xyz_ray,D_ray,Te_ray,ne_ray,zTi_ray,zni_ray,BRHO_ray,BP_ray,Bz_ray,B_ray,r_ray,salpha_ray,gradrho_ray,Upsilon_ray,pmode_ray,wpe2,wce2] = rayparameters_yp(rayparam,equil_fit,fluct_fit,omega_ray,rho_ray,theta_ray,z_ray,krho_ray,m_ray,kz_ray);
0114     %
0115     alphaphi_ray = [];
0116     tau_ray = [];
0117     eikval_ray = [];
0118     Upsilonn_ray = Upsilon_ray/equil_fit.ap;
0119     gradrhon_ray = gradrho_ray*equil_fit.ap;
0120     %
0121 end
0122 %
0123 if ~isempty(yproc) & rayparam.tau_lim > 0,
0124     %
0125     iabs_out = length(Npar_ray);% ray stopped after absorption
0126     %
0127 else
0128     %
0129     iabs_out = min([min(find(6.5./sqrt(Te_ray) < abs(Npar_ray))),length(Npar_ray)]);
0130     %
0131 end
0132 %
0133 pmode_ray = logical(pmode_ray);
0134 %
0135 if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation. Calculate the exact theoretical solution
0136     cnorm = clum/omega_ray(1)/equil_fit.ap;
0137     %
0138     if isinf(equil_fit.Rp),
0139         Rpn = Inf;
0140     else
0141         Rpn = equil_fit.Rp/equil_fit.ap;%Normalized major radius
0142     end
0143     %
0144     s = linspace(0,10,1001);%Ray length
0145     omega_vac = ones(length(s),1)*omega_ray(1);
0146     %
0147     Mrot0 = [cos(theta_ray(1))*cos(z_ray(1)/equil_fit.Rp),-sin(theta_ray(1))*cos(z_ray(1)/equil_fit.Rp),-sin(z_ray(1)/equil_fit.Rp);cos(theta_ray(1))*sin(z_ray(1)/equil_fit.Rp),-sin(theta_ray(1))*sin(z_ray(1)/equil_fit.Rp),cos(z_ray(1)/equil_fit.Rp);sin(theta_ray(1)),cos(theta_ray(1)),0];
0148     N_ray0 = cnorm*[krho_ray(1);m_ray(1)*equil_fit.ap/r_ray(1);kz_ray(1)*equil_fit.ap/Upsilon_ray(1)];%Initial wave vector in contravariant coordinates
0149     N_vac0_c = Mrot0*N_ray0;%Initial wave vector in cartesian coordinates
0150     %
0151     z_vac(1) = z_ray(1)/equil_fit.ap;
0152     rho_vac(1) = rho_ray(1);
0153     theta_vac(1) = theta_ray(1);
0154     rn_vac(1) = r_ray(1)/equil_fit.ap;
0155     xn_vac(1) = x_ray(1)/equil_fit.ap;
0156     krho_vac(1) = krho_ray(1);
0157     m_vac(1) = m_ray(1);
0158     kz_vac(1) = kz_ray(1);
0159     %
0160     if isinf(equil_fit.Rp),
0161         n_vac(1) = NaN;
0162     else
0163         n_vac(1) = kz_vac(1)*equil_fit.Rp;
0164     end    
0165     %
0166     if isinf(equil_fit.Rp),
0167         X_vac0_c = [rho_vac(1)*cos(theta_vac(1));z_vac(1);rho_vac(1)*sin(theta_vac(1))];%Initial ray position in Cartesian coordinates
0168     else
0169         X_vac0_c = [(Rpn + rho_vac(1)*cos(theta_vac(1)))*cos(z_vac(1)/Rpn);(Rpn + rho_vac(1)*cos(theta_vac(1)))*sin(z_vac(1)/Rpn);rho_vac(1)*sin(theta_vac(1))];%Initial ray position in Cartesian coordinates
0170     end    
0171     %
0172     for ii = 2:length(s),
0173         X_vac_c = X_vac0_c(:,ii-1) + N_vac0_c*(s(ii) - s(ii-1));%Ray equation in cartesian coordinates at next position
0174         %
0175         if isinf(equil_fit.Rp),
0176             z_vac(ii,1) = X_vac_c(ii,2);%New ray equation in toroidal coordinates
0177             rho_vac(ii,1) = sqrt(X_vac_c(1)^2 + X_vac_c(3)^2);%New ray equation in toroidal coordinates
0178             theta_vac(ii,1) = asin(X_vac_c(3)/rho_vac(ii));%New ray equation in toroidal coordinates
0179         else
0180             z_vac(ii,1) = Rpn*(atan(X_vac_c(2)/X_vac_c(1)) + pi*(X_vac_c(1) < 0));%New ray equation in toroidal coordinates
0181             rho_vac(ii,1) = sqrt((X_vac_c(1)/cos(z_vac(ii)/Rpn) - Rpn)^2 + X_vac_c(3)^2);%New ray equation in toroidal coordinates
0182             theta_vac(ii,1) = asin(X_vac_c(3)/rho_vac(ii))*(sqrt(X_vac_c(1)^2+X_vac_c(2)^2) > Rpn) + (pi*sign(theta_vac(ii-1)) - asin(X_vac_c(3)/rho_vac(ii)))*(sqrt(X_vac_c(1)^2+X_vac_c(2)^2) < Rpn);%New ray equation in toroidal coordinates
0183         end
0184         %
0185         if rho_vac(ii) >=1,%Specular reflexion is enforced at plasma boundary
0186             Mrot0 = [cos(theta_vac(ii-1))*cos(z_vac(ii-1)/Rpn),cos(theta_vac(ii-1))*sin(z_vac(ii-1)/Rpn),+sin(theta_vac(ii-1));-sin(theta_vac(ii-1))*cos(z_vac(ii-1)/Rpn),-sin(theta_vac(ii-1))*sin(z_vac(ii-1)/Rpn),cos(theta_vac(ii-1));sin(z_vac(ii-1)/Rpn),-cos(z_vac(ii-1)/Rpn),0];
0187             Mrot1 = [-cos(theta_vac(ii-1))*cos(z_vac(ii-1)/Rpn),-cos(theta_vac(ii-1))*sin(z_vac(ii-1)/Rpn),-sin(theta_vac(ii-1));-sin(theta_vac(ii-1))*cos(z_vac(ii-1)/Rpn),-sin(theta_vac(ii-1))*sin(z_vac(ii-1)/Rpn),cos(theta_vac(ii-1));sin(z_vac(ii-1)/Rpn),-cos(z_vac(ii-1)/Rpn),0];
0188             N_vac0_c = inv(Mrot0)*(Mrot1*N_vac0_c);%Ray direction is updated
0189             X_vac_c = X_vac0_c(:,ii-1) + N_vac0_c*(s(ii) - s(ii-1));%Ray equation in cartesian coordinates
0190             %
0191             if isinf(equil_fit.Rp),
0192                 z_vac(ii,1) = X_vac_c(ii,2);%New ray equation in toroidal coordinates
0193                 rho_vac(ii,1) = sqrt(X_vac_c(1)^2 + X_vac_c(3)^2);%New ray equation in toroidal coordinates
0194                 theta_vac(ii,1) = asin(X_vac_c(3)/rho_vac(ii));%New ray equation in toroidal coordinates
0195             else
0196                 r_vac(ii) = Rpn*(atan(X_vac_c(2)/X_vac_c(1)) + pi*(X_vac_c(1) < 0));%New ray equation in toroidal coordinates
0197                 rho_vac(ii) = sqrt((X_vac_c(1)/cos(z_vac(ii)/Rpn) - Rpn)^2 + X_vac_c(3)^2);%New ray equation in toroidal coordinates
0198                 theta_vac(ii) = asin(X_vac_c(3)/rho_vac(ii))*(sqrt(X_vac_c(1)^2+X_vac_c(2)^2) > Rpn) + (pi*sign(theta_vac(ii-1)) - asin(X_vac_c(3)/rho_vac(ii)))*(sqrt(X_vac_c(1)^2+X_vac_c(2)^2) < Rpn);%New ray equation in toroidal coordinates
0199             end
0200         end
0201         %
0202         rn_vac(ii,1) = rho_vac(ii);%circular concentric magnetic flux surfaces
0203         %
0204         invMrot = inv([cos(theta_vac(ii))*cos(z_vac(ii)/Rpn),-sin(theta_vac(ii))*cos(z_vac(ii)/Rpn),-sin(z_vac(ii)/Rpn);cos(theta_vac(ii))*sin(z_vac(ii)/Rpn),-sin(theta_vac(ii))*sin(z_vac(ii)/Rpn),cos(z_vac(ii)/Rpn);sin(theta_vac(ii)),cos(theta_vac(ii)),0]);
0205         N_vac = invMrot*N_vac0_c;
0206         %
0207         krho_vac(ii,1) = N_vac(1)/cnorm;
0208         m_vac(ii,1) = N_vac(2)*rn_vac(ii,1)/cnorm;       
0209         %
0210         if ~isinf(equil_fit.Rp),
0211             Rn_vac(ii,1) = sqrt(X_vac_c(1)^2 + X_vac_c(2)^2);
0212             kz_vac(ii,1) = N_vac(3)*Rn_vac(ii,1)/cnorm/equil_fit.Rp;
0213             n_vac(ii,1) = kz_vac(ii,1)*equil_fit.Rp;
0214             phi_vac(ii,1) = z_vac(ii,1)/Rpn;
0215         else
0216             kz_vac(ii,1) = N_vac(3)/cnorm/equil_fit.ap;
0217             n_vac(ii,1) = NaN;
0218             phi_vac(ii,1) = NaN;
0219         end
0220         %
0221         X_vac0_c(:,ii) = X_vac_c;
0222     end
0223     %
0224     [dummy,dummy,Npar_vac,Nperp_vac,epol_pmz_vac,phi_xyz_vac,D_vac,dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy] = rayparameters_yp(rayparam,equil_fit,fluct_fit,omega_vac,rho_vac,theta_vac,z_vac*equil_fit.ap,krho_vac,m_vac,kz_vac);
0225     %
0226     eikval_vac = [];
0227     %
0228 end
0229 %
0230 % Build magnetic equilibrium for display
0231 %
0232 if display_mode,
0233     %
0234     ntheta_fit = 65;
0235     rho_fit = linspace(0,1,15);
0236     theta_fit = linspace(0,2*pi,ntheta_fit);
0237     x_a0_fit = ppval(equil_fit.x_fit.pp_a0,rho_fit);
0238     x_an_fit = ppval(equil_fit.x_fit.pp_an,rho_fit);
0239     x_bn_fit = ppval(equil_fit.x_fit.pp_bn,rho_fit);
0240     y_a0_fit = ppval(equil_fit.y_fit.pp_a0,rho_fit);
0241     y_an_fit = ppval(equil_fit.y_fit.pp_an,rho_fit);
0242     y_bn_fit = ppval(equil_fit.y_fit.pp_bn,rho_fit);
0243     %
0244     if ~isinf(equil_fit.Rp),% For toroidal display
0245         %
0246         phi_fit = linspace(0,2*pi,ntheta_fit);
0247         %
0248         for ii = 1:ntheta_fit,
0249             R_fit(:,ii) = equil_fit.Rp + x_a0_fit.' + x_an_fit*cos(equil_fit.x_fit.n*theta_fit(ii)) + x_bn_fit*sin(equil_fit.x_fit.n*theta_fit(ii));
0250             Z_fit(:,ii) = equil_fit.Zp + y_a0_fit.' + y_an_fit*cos(equil_fit.y_fit.n*theta_fit(ii)) + y_bn_fit*sin(equil_fit.y_fit.n*theta_fit(ii));
0251         end
0252         %
0253         Rmin = min(min(R_fit));
0254         Rmax = max(max(R_fit));  
0255         Zmin = min(min(Z_fit));
0256         Zmax = max(max(Z_fit));
0257         %
0258         phi_ray = z_ray/equil_fit.Rp;
0259         %
0260         X_ray = (x_ray + equil_fit.Rp).*cos(phi_ray);
0261         Y_ray = - (x_ray  + equil_fit.Rp).*sin(phi_ray);%Note : phi is clockwise from top
0262         Xmin_fit = Rmin*cos(phi_fit);   
0263         Ymin_fit = Rmin*sin(phi_fit); 
0264         Xmax_fit = Rmax*cos(phi_fit);   
0265         Ymax_fit = Rmax*sin(phi_fit); 
0266         Xp_fit = equil_fit.Rp*cos(phi_fit);   
0267         Yp_fit = equil_fit.Rp*sin(phi_fit); 
0268         ss = sign(diff(rho_ray));
0269         for iss = 1:length(ss),
0270             if ss(iss)==0,ss(iss) = ss(iss-1);end
0271         end
0272         %
0273         % Display results
0274         %
0275         flag_vacuum = 0;
0276         figure('Name','Poloidal ray trajectory')
0277         [ax] = graph1D_jd(x_ray(pmode_ray) + equil_fit.Rp,y_ray(pmode_ray) + equil_fit.Zp,0,0,'R(m)','Z(m)',['Poloidal Ray trajectory',subtitle],NaN,[Rmin,Rmax],[Zmin,Zmax],'none','.','r',2);hold on
0278         [ax] = graph1D_jd(x_ray(~pmode_ray) + equil_fit.Rp,y_ray(~pmode_ray) + equil_fit.Zp,0,0,'R(m)','Z(m)',['Poloidal Ray trajectory',subtitle],NaN,[Rmin,Rmax],[Zmin,Zmax],'none','.','b',2);hold on
0279     %
0280         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0281            [ax] = graph1D_jd(sqrt(X_vac0_c(1,:).^2+X_vac0_c(2,:).^2)*equil_fit.ap,X_vac0_c(3,:)*equil_fit.ap,0,0,'R(m)','Z(m)',['Ray trajectory',subtitle],NaN,[Rmin,Rmax],[Zmin,Zmax],'--','none','r',2);hold on
0282            legend('Ray-tracing','Straight lines model');
0283            flag_vacuum = 1;
0284         else
0285            [ax] = graph1D_jd(x_ray(iabs_out:end) + equil_fit.Rp,y_ray(iabs_out:end) + equil_fit.Zp,0,0,'R(m)','Z(m)',['Ray trajectory',subtitle],NaN,[Rmin,Rmax],[Zmin,Zmax],'-','none','k',2);hold on 
0286         end
0287         [ax] = graph1D_jd(R_fit,Z_fit,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],'-','none','b',1);drawnow
0288         [ax] = graph1D_jd(R_fit',Z_fit',0,0,'R(m)','Z(m)',['Ray trajectory',subtitle],NaN,[Rmin,Rmax],[Zmin,Zmax],'-','none','b',1);hold on
0289     %
0290         if display_text & flag_vacuum == 0,
0291           for ii = 1:length(x_ray)-2,
0292              if (ss(ii) ~= ss(ii+1)),
0293                 text(x_ray(ii) + equil_fit.Rp,x_ray(ii) + equil_fit.Zp,num2str(ii));
0294              end
0295           end
0296         end
0297         %
0298         axis('equal');drawnow  
0299         print_jd(p_opt,[filename,'_RayPoloidal']);
0300         %
0301         figure('Name','Toroidal ray trajectory')
0302         [ax] = graph1D_jd(X_ray(pmode_ray),Y_ray(pmode_ray),0,0,'X(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-Rmax,Rmax],[-Rmax,Rmax],'none','.','r',2);hold on
0303         [ax] = graph1D_jd(X_ray(~pmode_ray),Y_ray(~pmode_ray),0,0,'X(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-Rmax,Rmax],[-Rmax,Rmax],'none','.','b',2);hold on
0304         %
0305         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0306             [ax] = graph1D_jd(X_vac0_c(1,:)*equil_fit.ap,-X_vac0_c(2,:)*equil_fit.ap,0,0,'X(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-Rmax,Rmax],[-Rmax,Rmax],'--','none','r',2);hold on
0307             legend('Ray-tracing','Straight lines model');
0308         else
0309             [ax] = graph1D_jd(X_ray(iabs_out:end),Y_ray(iabs_out:end),0,0,'X(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-Rmax,Rmax],[-Rmax,Rmax],'-','none','k',2);hold on 
0310         end
0311         [ax] = graph1D_jd(Xmin_fit,Ymin_fit,0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],'-','none','b',1);drawnow
0312         [ax] = graph1D_jd(Xmax_fit,Ymax_fit,0,0,'X(m)','Y(m)','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],'-','none','b',1);hold on
0313         [ax] = graph1D_jd(Xp_fit,Yp_fit,0,0,'X(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-Rmax,Rmax],[-Rmax,Rmax],'--','none','b',1);hold on
0314         %
0315         if display_text & flag_vacuum == 0,
0316             for ii = 1:length(X_ray)-2,
0317                 if (ss(ii) ~= ss(ii+1)),
0318                     text(X_ray(ii),Y_ray(ii),num2str(ii));
0319                 end
0320             end
0321         end
0322         %
0323         axis('equal');drawnow  
0324         print_jd(p_opt,[filename,'_RayToroidal']);
0325         %
0326     else,% for cylindrical display
0327         %
0328         theta_fit = linspace(0,2*pi,ntheta_fit);
0329         %
0330         for ii = 1:ntheta_fit,
0331             x_fit(:,ii) = x_a0_fit.' + x_an_fit*cos(equil_fit.x_fit.n*theta_fit(ii)) + x_bn_fit*sin(equil_fit.x_fit.n*theta_fit(ii));
0332             y_fit(:,ii) = y_a0_fit.' + y_an_fit*cos(equil_fit.y_fit.n*theta_fit(ii)) + y_bn_fit*sin(equil_fit.y_fit.n*theta_fit(ii));
0333         end
0334         %
0335         xmin = min(min(x_fit));
0336         xmax = max(max(x_fit)); 
0337         Zmin = min(min(y_fit)) + equil_fit.Zp;
0338         Zmax = max(max(y_fit)) + equil_fit.Zp;
0339         %
0340         %phi_ray = z_ray/equil_fit.Rp;
0341         %
0342         X_ray = x_ray;
0343         Y_ray = - y_ray;%Note : phi is clockwise from top
0344         Xmin_fit = -equil_fit.ap*cos(theta_fit);   
0345         Ymin_fit = -equil_fit.ap*sin(theta_fit); 
0346         Xmax_fit = equil_fit.ap*cos(theta_fit);   
0347         Ymax_fit = equil_fit.ap*sin(theta_fit); 
0348         Zmin_fit = linspace(-equil_fit.ap,equil_fit.ap,ntheta_fit);
0349         Xp_fit = 0*cos(theta_fit);   
0350         Yp_fit = 0*sin(theta_fit);
0351         Zp_fit = 0*sin(theta_fit); 
0352         ss = sign(diff(rho_ray));
0353         for iss = 1:length(ss),
0354             if ss(iss)==0,ss(iss) = ss(iss-1);end
0355         end
0356         %
0357         % Display results
0358         %
0359         flag_vacuum = 0;
0360         figure('Name','Poloidal ray trajectory')
0361         [ax] = graph1D_jd(x_ray(pmode_ray),y_ray(pmode_ray) + equil_fit.Zp,0,0,'x(m)','Z(m)',['Poloidal Ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'none','.','r',2);hold on
0362         [ax] = graph1D_jd(x_ray(~pmode_ray),y_ray(~pmode_ray) + equil_fit.Zp,0,0,'x(m)','Z(m)',['Poloidal Ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'none','.','b',2);hold on
0363         %
0364         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0365            [ax] = graph1D_jd(sqrt(X_vac0_c(1,:).^2+X_vac0_c(2,:).^2)*equil_fit.ap,X_vac0_c(3,:)*equil_fit.ap,0,0,'x(m)','Z(m)',['Ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'--','none','r',2);hold on
0366            legend('Ray-tracing','Straight lines model');
0367            flag_vacuum = 1;
0368         else
0369            [ax] = graph1D_jd(x_ray(iabs_out:end),y_ray(iabs_out:end) + equil_fit.Zp,0,0,'x(m)','Z(m)',['Ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'-','none','k',2);hold on 
0370         end
0371         [ax] = graph1D_jd(x_fit,y_fit + equil_fit.Zp,0,0,'','','',NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'-','none','b',1);drawnow
0372         [ax] = graph1D_jd(x_fit',y_fit' + equil_fit.Zp,0,0,'x(m)','Z(m)',['Ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[Zmin,Zmax],'-','none','b',1);hold on
0373         %
0374         if display_text & flag_vacuum == 0,
0375             for ii = 1:length(x_ray)-2,
0376                 if (ss(ii) ~= ss(ii+1)),
0377                     text(x_ray(ii),x_ray(ii) + equil_fit.Zp,num2str(ii));
0378                 end
0379             end
0380         end
0381         %
0382         axis('equal');drawnow  
0383         print_jd(p_opt,[filename,'_RayPoloidal']);
0384         %
0385         figure('Name','Axial ray trajectory')
0386         [ax] = graph1D_jd(x_ray(pmode_ray),-z_ray(pmode_ray),0,0,'x(m)','Y(m)',['Axial ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[min(-z_ray),max(-z_ray)],'none','.','r',2);hold on
0387         [ax] = graph1D_jd(x_ray(~pmode_ray),-z_ray(~pmode_ray),0,0,'x(m)','Y(m)',['Axial ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[min(-z_ray),max(-z_ray)],'none','.','b',2);hold on
0388         %
0389         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0390             [ax] = graph1D_jd(X_vac0_c(1,:)*equil_fit.ap,-X_vac0_c(2,:)*equil_fit.ap,0,0,'x(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[min(z_ray),max(z_ray)],'--','none','r',2);hold on
0391             legend('Ray-tracing','Straight lines model');
0392         else
0393             [ax] = graph1D_jd(X_ray(iabs_out:end),Y_ray(iabs_out:end),0,0,'x(m)','Y(m)',['Axial ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[min(-z_ray),max(-z_ray)],'-','none','k',2);hold on 
0394         end
0395         %
0396         for ii = 1:ntheta_fit,
0397             [ax] = graph1D_jd([Xmin_fit(ii),Xmin_fit(ii)],[min(-z_ray),max(-z_ray)],0,0,'x(m)','Y(m)',['Toroidal ray trajectory',subtitle],NaN,[-equil_fit.ap,equil_fit.ap],[min(-z_ray),max(-z_ray)],'-','none','b',1);drawnow        
0398         end
0399         %
0400         if display_text & flag_vacuum == 0,
0401             for ii = 1:length(X_ray)-2,
0402                 if (ss(ii) ~= ss(ii+1)),
0403                     text(x_ray(ii),z_ray(ii),num2str(ii));
0404                 end
0405             end
0406         end
0407         %
0408         drawnow  
0409         print_jd(p_opt,[filename,'_RayAxial']);
0410         %
0411     end
0412     %
0413     if ~isinf(equil_fit.Rp),% For toroidal display
0414         figure('Name','Parallel wave refractive index')
0415         graph1D_jd(abs(phi_ray)/2/pi,[real(Npar_ray),imag(Npar_ray)*sign(real(Npar_ray(1)))],0,0,'(2*\phi_0+|\phi|)/2\pi','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'-','none',NaN,2);drawnow;hold off
0416         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0417           [ax] = graph1D_jd(phi_vac/2/pi,Npar_vac,0,0,'\phi/2\pi','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'--','none','r',2);hold on
0418           legend('Ray-tracing','Straight lines model');
0419         else
0420           graph1D_jd(abs(phi_ray)/2/pi,sign(real(Npar_ray(1)))*6.5./sqrt(Te_ray),0,0,'(2*\phi_0+|\phi|)/2\pi','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'-','none','r',2);drawnow;hold off
0421           legend('Ray-tracing','Accessibility','Landau damping');
0422         end
0423         print_jd(p_opt,[filename,'_Npar']);
0424     else,% for cylindrical display
0425         figure('Name','Parallel wave refractive index')
0426         graph1D_jd(abs(z_ray),[real(Npar_ray),imag(Npar_ray)*sign(real(Npar_ray(1)))],0,0,'|z|(m)','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'-','none',NaN,2);drawnow;hold off
0427         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0428           [ax] = graph1D_jd(abs(z_ray),Npar_vac,0,0,'z(m)','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'--','none','r',2);hold on
0429           legend('Ray-tracing','Straight lines model');
0430         else
0431           graph1D_jd(abs(z_ray),sign(real(Npar_ray(1)))*6.5./sqrt(Te_ray),0,0,'|z|(m)','N_{||}',['Parallel wave refractive index',subtitle],NaN,'',[min(real(Npar_ray)),max(real(Npar_ray))],'-','none','r',2);drawnow;hold off
0432           legend('Ray-tracing','Accessibility','Landau damping');
0433         end
0434         print_jd(p_opt,[filename,'_Npar']);
0435     end
0436     %
0437     if ~isinf(equil_fit.Rp),% For toroidal display
0438         figure('Name','Perpendicular wave refractive index')
0439         graph1D_jd(abs(phi_ray)/2/pi,Nperp_ray,0,0,'|\phi|/(2\pi)','N_{perp}',['Perpendicular wave refractive index',subtitle],NaN,'','','-','none','k',2);drawnow;hold off
0440         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0441            [ax] = graph1D_jd(phi_vac/2/pi,Nperp_vac,0,0,'\phi/2\pi','N_{perp}',['Perpendicular wave refractive index',subtitle],NaN,'','','--','none','r',2);hold on
0442            legend('Ray-tracing','Straight lines model');
0443         end
0444         print_jd(p_opt,[filename,'_Nperp']);
0445     else,% for cylindrical display
0446         figure('Name','Perpendicular wave refractive index')
0447         graph1D_jd(abs(z_ray),Nperp_ray,0,0,'z(m)','N_{perp}',['Perpendicular wave refractive index',subtitle],NaN,'','','-','none','k',2);drawnow;hold off
0448         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0449            [ax] = graph1D_jd(abs(z_ray),Nperp_vac,0,0,'|z|(m)','N_{perp}',['Perpendicular wave refractive index',subtitle],NaN,'','','--','none','r',2);hold on
0450            legend('Ray-tracing','Straight lines model','location','NorthWest');
0451         end
0452         print_jd(p_opt,[filename,'_Nperp']);    
0453     end   
0454     %
0455     if ~isinf(equil_fit.Rp),% For toroidal display
0456         figure('Name','Poloidal mode number')
0457         graph1D_jd(abs(phi_ray)/2/pi,yout(:,5),0,0,'|\phi|/(2\pi)','m',['Poloidal mode number',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0458         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0459            [ax] = graph1D_jd(phi_vac/2/pi,m_vac,0,0,'\phi/2\pi','m',['Poloidal mode number',subtitle],NaN,'','','--','none','r',2);hold on
0460            legend('Ray-tracing','Straight lines model');
0461         end
0462         print_jd(p_opt,[filename,'_PoloidalMode']);
0463     else,% for cylindrical display
0464         figure('Name','Poloidal mode number')
0465         graph1D_jd(abs(z_ray),yout(:,5),0,0,'z(m)','m',['Poloidal mode number',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0466         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0467            [ax] = graph1D_jd(abs(z_ray),m_vac,0,0,'|z|(m)','m',['Poloidal mode number',subtitle],NaN,'','','--','none','r',2);hold on
0468            legend('Ray-tracing','Straight lines model');
0469         end
0470         print_jd(p_opt,[filename,'_PoloidalMode']);        
0471     end
0472     %
0473     if ~isinf(equil_fit.Rp),% For toroidal display
0474         figure('Name','Toroidal mode number')
0475         graph1D_jd(abs(phi_ray)/2/pi,yout(:,6)*equil_fit.Rp,0,0,'|\phi|/(2\pi)','n',['Toroidal mode number',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0476         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0477            [ax] = graph1D_jd(phi_vac/2/pi,n_vac,0,0,'\phi/2\pi','n',['Toroidal mode number',subtitle],NaN,'',[0,max(yout(1,6))*1.1],'--','none','r',2);hold on
0478            legend('Ray-tracing','Straight lines model');
0479         end
0480         print_jd(p_opt,[filename,'_ToroidalMode']);
0481     else,% for cylindrical display
0482         figure('Name','Axial wave vector')
0483         graph1D_jd(abs(z_ray),yout(:,6),0,0,'|z|(m)','kz (m-1)',['Axial wave vector',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0484         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0485            [ax] = graph1D_jd(abs(z_ray),kz_vac,0,0,'|z|(m)','kz (m-1)',['Axial wave vector',subtitle],NaN,'',[0,max(yout(1,6))*1.1],'--','none','r',2);hold on
0486            legend('Ray-tracing','Straight lines model');
0487         end
0488         print_jd(p_opt,[filename,'_AxialWaveVector']); 
0489     end
0490     %
0491     if ~isinf(equil_fit.Rp),% For toroidal display
0492         figure('Name','Poloidal angle')
0493         graph1D_jd(abs(phi_ray)/2/pi,yout(:,2)/2/pi,0,0,'|\phi|/(2\pi)','\theta/2\pi',['Poloidal angle',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0494         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0495            [ax] = graph1D_jd(phi_vac/2/pi,theta_vac/2/pi,0,0,'\phi/2\pi','(\theta/2\pi)',['Poloidal angle',subtitle],NaN,'','','--','none','r',2);hold on
0496            legend('Ray-tracing','Straight lines model');
0497         end
0498         print_jd(p_opt,[filename,'_Theta']);
0499     else,% for cylindrical display
0500         figure('Name','Poloidal angle')
0501         graph1D_jd(abs(z_ray),yout(:,2)/2/pi,0,0,'|z|(m)','\theta/2\pi',['Poloidal angle',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0502         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0503            [ax] = graph1D_jd(abs(z_ray),theta_vac/2/pi,0,0,'|z|(m)','\theta/(2\pi)',['Poloidal angle',subtitle],NaN,'','','--','none','r',2);hold on
0504            legend('Ray-tracing','Straight lines model');
0505         end
0506         print_jd(p_opt,[filename,'_Theta']);
0507     end
0508     %
0509     if ~isinf(equil_fit.Rp),% For toroidal display
0510         figure('Name','Radial position')
0511         graph1D_jd(abs(phi_ray)/2/pi,yout(:,1),0,0,'|\phi|/(2\pi)','rho',['Normalized radial position',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0512         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0513            [ax] = graph1D_jd(phi_vac/2/pi,rho_vac,0,0,'\phi/(2\pi)','rho',['Normalized radial position',subtitle],NaN,'','','--','none','r',2);hold on
0514            legend('Ray-tracing','Straight lines model');
0515         end
0516         print_jd(p_opt,[filename,'_Rho']);
0517     else,% for cylindrical display
0518         figure('Name','Radial position')
0519         graph1D_jd(abs(z_ray),yout(:,1),0,0,'|z|(m)','rho',['Normalized radial position',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0520         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0521            [ax] = graph1D_jd(abs(z_ray),rho_vac,0,0,'|z|(m)','rho',['Normalized radial position',subtitle],NaN,'','','--','none','r',2);hold on
0522            legend('Ray-tracing','Straight lines model');
0523         end
0524         print_jd(p_opt,[filename,'_Rho']);   
0525     end
0526     %
0527     if ~isinf(equil_fit.Rp),% For toroidal display
0528         figure('Name','Normalized radial wave vector')
0529         graph1D_jd(abs(phi_ray)/2/pi,yout(:,4),0,0,'|\phi|/(2\pi)','krho*ap',['Normalized radial wavenumber',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0530         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0531            [ax] = graph1D_jd(phi_vac/2/pi,krho_vac,0,0,'\phi/2\pi','krho*ap',['Normalized radial wavenumber',subtitle],NaN,'','','--','none','r',2);hold on
0532            legend('Ray-tracing','Straight lines model');
0533         end
0534         print_jd(p_opt,[filename,'_RadialWaveVector']);
0535     else,% for cylindrical display
0536         figure('Name','Normalized radial wave vector')
0537         graph1D_jd(abs(z_ray),yout(:,4),0,0,'|z|(m)','krho*ap',['Normalized radial wave vector',subtitle],NaN,'','','-','none','k',2);drawnow;hold on
0538         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0539            [ax] = graph1D_jd(abs(z_ray),krho_vac,0,0,'|z|(m)','krho*ap',['Normalized radial wave vector',subtitle],NaN,'','','--','none','r',2);hold on
0540            legend('Ray-tracing','Straight lines model');
0541         end
0542         print_jd(p_opt,[filename,'_RadialWaveVector']);    
0543     end
0544     %
0545     if ~isinf(equil_fit.Rp),% For toroidal display
0546         figure('Name','Dispersion relation'),
0547         graph1D_jd(abs(phi_ray(pmode_ray))/2/pi,abs(D_ray(pmode_ray)),0,1,'|\phi|/(2\pi)','D',['Dispersion relation',subtitle],NaN,'','','none','.','r',2);
0548         graph1D_jd(abs(phi_ray(~pmode_ray))/2/pi,abs(D_ray(~pmode_ray)),0,1,'|\phi|/(2\pi)','D',['Dispersion relation',subtitle],NaN,'','','none','.','b',2);drawnow
0549         if isempty(phi_ray(pmode_ray))
0550            legend({'mode m'})
0551         elseif isempty(phi_ray(~pmode_ray))
0552            legend({'mode p'})
0553         else
0554            legend({'mode p','mode m'})
0555         end
0556         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0557           [ax] = graph1D_jd(phi_vac/2/pi,D_vac,0,0,'\phi/2\pi','D',['Dispersion relation',subtitle],NaN,'','','--','none','r',2);hold on
0558           legend('Ray-tracing','Straight lines model');
0559         end
0560         print_jd(p_opt,[filename,'_Disp']);
0561     else,% for cylindrical display
0562         figure('Name','Dispersion relation'),
0563         graph1D_jd(abs(z_ray(pmode_ray)),abs(D_ray(pmode_ray)),0,1,'|z|(m)','D',['Dispersion relation',subtitle],NaN,'','','none','.','r',2);
0564         graph1D_jd(abs(z_ray(~pmode_ray)),abs(D_ray(~pmode_ray)),0,1,'|z|(m)','D',['Dispersion relation',subtitle],NaN,'','','none','.','b',2);drawnow
0565         if isempty(z_ray(pmode_ray))
0566            legend({'mode m'})
0567         elseif isempty(z_ray(~pmode_ray))
0568            legend({'mode p'})
0569         else
0570            legend({'mode p','mode m'})
0571         end
0572         if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0573           [ax] = graph1D_jd(abs(z_vac),D_vac,0,0,'|z|(m)','D',['Dispersion relation',subtitle],NaN,'','','--','none','r',2);hold on
0574           legend('Ray-tracing','Straight lines model');
0575         end
0576         print_jd(p_opt,[filename,'_Disp']);        
0577     end      
0578     %
0579     if ~isinf(equil_fit.Rp),% For toroidal display
0580         if ~isempty(eikval_ray),
0581             figure('Name','Validity of the eikonal approximation'),
0582             graph1D_jd(abs(phi_ray(pmode_ray))/2/pi,abs(eikval_ray(pmode_ray)),0,1,'|\phi|/(2\pi)','\nabla\cdotN/N^2',['Validity of the eikonal approximation (< 1)',subtitle],NaN,'','','none','.','r',2);
0583             graph1D_jd(abs(phi_ray(~pmode_ray))/2/pi,abs(eikval_ray(~pmode_ray)),0,1,'|\phi|/(2\pi)','\nabla\cdotN/N^2',['Validity of the eikonal approximation (< 1)',subtitle],NaN,'','','none','.','b',2);drawnow
0584             if isempty(phi_ray(pmode_ray))
0585                 legend({'mode m'})
0586             elseif isempty(phi_ray(~pmode_ray))
0587                 legend({'mode p'})
0588             else
0589                 legend({'mode p','mode m'})
0590             end
0591             print_jd(p_opt,[filename,'_Eikval']);
0592         end
0593         %
0594         if exist('eikval_vac'),
0595            if ~isempty(eikval_vac),
0596                 if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0597                     [ax] = graph1D_jd(phi_vac/2/pi,eikval_vac,0,0,'\phi/2\pi','\nabla\cdotN/N^2',['\nabla\cdotN/N^2',subtitle],NaN,'','','--','none','r',2);hold on
0598                     legend('Ray-tracing','Straight lines model');
0599                 end
0600                 print_jd(p_opt,[filename,'_Eikval']);
0601            end
0602         end
0603     else,% for cylindrical display
0604         if ~isempty(eikval_ray),
0605             figure('Name','Validity of the eikonal approximation'),
0606             graph1D_jd(abs(z_ray(pmode_ray)),abs(eikval_ray(pmode_ray)),0,1,'|z|(m)','\nabla\cdotN/N^2',['Validity of the eikonal approximation (< 1)',subtitle],NaN,'','','none','.','r',2);
0607             graph1D_jd(abs(z_ray(~pmode_ray)),abs(eikval_ray(~pmode_ray)),0,1,'|z|(m)','\nabla\cdotN/N^2',['Validity of the eikonal approximation (< 1)',subtitle],NaN,'','','none','.','b',2);drawnow
0608             if isempty(z_ray(pmode_ray))
0609                 legend({'mode m'})
0610             elseif isempty(z_ray(~pmode_ray))
0611                 legend({'mode p'})
0612             else
0613                 legend({'mode p','mode m'})
0614             end
0615             print_jd(p_opt,[filename,'_Eikval']);
0616         end
0617         %
0618         if exist('eikval_vac'),
0619            if ~isempty(eikval_vac),
0620                 if rem(abs(1-real(Npar_ray(1))^2-Nperp_ray(1)^2),1) <= 100*eps,% Vacuum propagation.
0621                     [ax] = graph1D_jd(abs(z_vac),eikval_vac,0,0,'|z|(m)','\nabla\cdotN/N^2',['\nabla\cdotN/N^2',subtitle],NaN,'','','--','none','r',2);hold on
0622                     legend('Ray-tracing','Straight lines model');
0623                 end
0624                 print_jd(p_opt,[filename,'_Eikval']);
0625            end
0626         end
0627         
0628     end
0629 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.