graph_comp_RT_jd

 This function compares C3PO results using numerical and analytical forms for the equilibrium.

0001 function graph_comp_RT_jd(rays,legs,equil_fit,filename,opt)
0002 %
0003 % This function compares C3PO results using numerical and analytical forms for the equilibrium.
0004 %
0005 if nargin == 0,
0006     disp('No ray to be displayed.')
0007     return
0008 else
0009     nb = length(rays);
0010 end
0011 %
0012 if nargin < 5,
0013     p_opt = -1;
0014     ntheta_fit = 65;
0015     nrho_fit = 15;
0016     propvar = 0;% (0) : ray.ss, (1) : ray.sphi/(2*pi)
0017 else
0018     p_opt = opt.p_opt;
0019     ntheta_fit = opt.ntheta_fit;
0020     nrho_fit = opt.nrho_fit;
0021     propvar = opt.propvar;% (0) : ray.ss, (1) : ray.sphi/(2*pi)
0022 end
0023 if nargin < 4,
0024     filename = 'Fig';
0025 end
0026 if nargin < 3,
0027     equil_fit = '';
0028 elseif ~isempty(equil_fit) && (~isstruct(equil_fit) || (length(equil_fit) ~= 1 && length(equil_fit) ~= nb)),
0029     error('equils_fit must be a structure array of length 1 or nb.')
0030 end
0031 if nargin < 2 || isempty(legs),
0032     legs = cell(2,nb);
0033     for ib = 1:nb,
0034         legs{1,ib} = [ 'Ray # ',num2str(ib)];% for figure legend
0035         legs{2,ib} = [ 'Ray_',num2str(ib)];% for filename
0036     end
0037 elseif size(legs,1) == 1,%automatic filenames
0038     for ib = 1:nb,
0039         legs{2,ib} = [ 'Ray_',num2str(ib)];% for filename
0040     end
0041 end
0042 %
0043 % For poloidal ray trajectory comparison
0044 %
0045 rho_fit = linspace(0,1,nrho_fit);
0046 theta_fit = linspace(0,2*pi,ntheta_fit);
0047 phi_fit = linspace(0,2*pi,ntheta_fit);
0048 %
0049 neq = length(equil_fit);
0050 %
0051 Rmin = NaN(1,neq);
0052 Rmax = NaN(1,neq);
0053 Zmin = NaN(1,neq);
0054 Zmax = NaN(1,neq);
0055 %
0056 for ieq = 1:neq,
0057     equil_fit(ieq).x_a0_fit = ppval(equil_fit(ieq).x_fit.pp_a0,rho_fit);
0058     equil_fit(ieq).x_an_fit = ppval(equil_fit(ieq).x_fit.pp_an,rho_fit);
0059     equil_fit(ieq).x_bn_fit = ppval(equil_fit(ieq).x_fit.pp_bn,rho_fit);
0060     equil_fit(ieq).y_a0_fit = ppval(equil_fit(ieq).y_fit.pp_a0,rho_fit);
0061     equil_fit(ieq).y_an_fit = ppval(equil_fit(ieq).y_fit.pp_an,rho_fit);
0062     equil_fit(ieq).y_bn_fit = ppval(equil_fit(ieq).y_fit.pp_bn,rho_fit);
0063     %
0064     equil_fit(ieq).R_fit = NaN(nrho_fit,ntheta_fit);
0065     equil_fit(ieq).Z_fit = NaN(nrho_fit,ntheta_fit);
0066     %
0067     for ii = 1:ntheta_fit,
0068         equil_fit(ieq).R_fit(:,ii) = equil_fit(ieq).Rp + equil_fit(ieq).x_a0_fit.' + ...
0069             equil_fit(ieq).x_an_fit*cos(equil_fit(ieq).x_fit.n*theta_fit(ii)) + equil_fit(ieq).x_bn_fit*sin(equil_fit(ieq).x_fit.n*theta_fit(ii));
0070         equil_fit(ieq).Z_fit(:,ii) = equil_fit(ieq).Zp + equil_fit(ieq).y_a0_fit.' + ...
0071             equil_fit(ieq).y_an_fit*cos(equil_fit(ieq).y_fit.n*theta_fit(ii)) + equil_fit(ieq).y_bn_fit*sin(equil_fit(ieq).y_fit.n*theta_fit(ii));
0072     end
0073     %
0074     Rmin(ieq) = min(min(equil_fit(ieq).R_fit));
0075     Rmax(ieq) = max(max(equil_fit(ieq).R_fit));
0076     Zmin(ieq) = min(min(equil_fit(ieq).Z_fit));
0077     Zmax(ieq) = max(max(equil_fit(ieq).Z_fit));
0078     %
0079 end
0080 %
0081 Rmin = min(Rmin);
0082 Rmax = max(Rmax);
0083 Zmin = min(Zmin);
0084 Zmax = max(Zmax);
0085 %
0086 style = '-';
0087 style1 = '--';
0088 style2 = 'none';
0089 marker = 'none';
0090 marker2 = '.';
0091 colors = {'r','b','g','m','c','y'};
0092 width = 0.5;
0093 width2 = 2;
0094 siz = 20+14i;
0095 red = 0.9;
0096 lspace = 0.7;
0097 bspace = 0.5;
0098 %
0099 ncol = length(colors);
0100 %
0101 if neq > 0,
0102     %
0103     Xmin_fit = Rmin*cos(phi_fit);   
0104     Ymin_fit = Rmin*sin(phi_fit); 
0105     Xmax_fit = Rmax*cos(phi_fit);   
0106     Ymax_fit = Rmax*sin(phi_fit); 
0107     %
0108     % --- all rays ---
0109     %
0110     figure('Name','Poloidal ray trajectories')
0111     axis('equal');
0112     for ib = 1:nb,
0113         ieq = min([ib,neq]);%equil_fit number
0114         graph1D_jd(rays{ib}.sx + equil_fit(ieq).Rp,rays{ib}.sy + equil_fit(ieq).Zp,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],style,marker,colors{rem(ib,ncol)},width2);
0115     end
0116     graph1D_jd(equil_fit(1).R_fit,equil_fit(1).Z_fit,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],style,marker,'k',width);
0117     graph1D_jd(equil_fit(1).R_fit.',equil_fit(1).Z_fit.',0,0,'R(m)','Z(m)','Poloidal Ray trajectories',legs(1,:),[Rmin,Rmax],[Zmin,Zmax],style,marker,'k',width,siz,gca,red,lspace,bspace);
0118     print_jd(p_opt,[filename,'_poloidal_trajectories']);
0119     %
0120     figure('Name','Toroidal ray trajectories')
0121     axis('equal');
0122     for ib = 1:nb,
0123         ieq = min([ib,neq]);%equil_fit number
0124         graph1D_jd((rays{ib}.sx + equil_fit(ieq).Rp).*cos(-rays{ib}.sphi),(rays{ib}.sx + equil_fit(ieq).Rp).*sin(-rays{ib}.sphi),0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style,marker,colors{rem(ib,ncol)},width2);
0125     end
0126     graph1D_jd(Xmin_fit,Ymin_fit,0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style,marker,'k',width);
0127     graph1D_jd(Xmax_fit,Ymax_fit,0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style,marker,'k',width);
0128     graph1D_jd(equil_fit(1).Rp*cos(phi_fit),equil_fit(1).Rp*sin(phi_fit),0,0,'X(m)','Y(m)','Toroidal Ray trajectories',legs(1,:),[-Rmax,Rmax],[-Rmax,Rmax],style1,marker,'k',width,siz,gca,red,lspace,bspace);
0129     print_jd(p_opt,[filename,'_toroidal_trajectories']);
0130     %
0131     % --- each ray, with polarization ---
0132     %
0133     leg = {'mode p','mode m'};
0134     for ib = 1:nb,
0135         figure('Name',['Poloidal ray trajectory - ',legs{2,ib}])
0136         axis('equal');
0137         ieq = min([ib,neq]);%equil_fit number
0138         graph1D_jd(rays{ib}.sx(rays{ib}.spmode == 1) + equil_fit(ieq).Rp,rays{ib}.sy(rays{ib}.spmode == 1) + equil_fit(ieq).Zp,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],style2,marker2,'r',width2);
0139         graph1D_jd(rays{ib}.sx(rays{ib}.spmode ~= 1) + equil_fit(ieq).Rp,rays{ib}.sy(rays{ib}.spmode ~= 1) + equil_fit(ieq).Zp,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],style2,marker2,'b',width2);
0140         graph1D_jd(equil_fit(ieq).R_fit,equil_fit(ieq).Z_fit,0,0,'','','',NaN,[Rmin,Rmax],[Zmin,Zmax],style,marker,'k',width);
0141         graph1D_jd(equil_fit(ieq).R_fit.',equil_fit(ieq).Z_fit.',0,0,'R(m)','Z(m)',['Poloidal Ray trajectory - ',legs{2,ib}],leg,[Rmin,Rmax],[Zmin,Zmax],style,marker,'k',width,siz,gca,red,lspace,bspace);
0142         print_jd(p_opt,[filename,'_poloidal_trajectory_',legs{2,ib}]);
0143         %
0144         figure('Name',['Toroidal ray trajectory - ',legs{2,ib}])
0145         axis('equal');
0146         ieq = min([ib,neq]);%equil_fit number
0147         graph1D_jd((rays{ib}.sx(rays{ib}.spmode == 1) + equil_fit(ieq).Rp).*cos(-rays{ib}.sphi(rays{ib}.spmode == 1)),(rays{ib}.sx(rays{ib}.spmode == 1) + equil_fit(ieq).Rp).*sin(-rays{ib}.sphi(rays{ib}.spmode == 1)),0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style2,marker2,'r',width2);
0148         graph1D_jd((rays{ib}.sx(rays{ib}.spmode ~= 1) + equil_fit(ieq).Rp).*cos(-rays{ib}.sphi(rays{ib}.spmode ~= 1)),(rays{ib}.sx(rays{ib}.spmode ~= 1) + equil_fit(ieq).Rp).*sin(-rays{ib}.sphi(rays{ib}.spmode ~= 1)),0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style2,marker2,'b',width2);
0149         graph1D_jd(Xmin_fit,Ymin_fit,0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style,marker,'k',width);
0150         graph1D_jd(Xmax_fit,Ymax_fit,0,0,'','','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],style,marker,'k',width);
0151         graph1D_jd(equil_fit(ieq).Rp*cos(phi_fit),equil_fit(ieq).Rp*sin(phi_fit),0,0,'X(m)','Y(m)',['Toroidal Ray trajectory - ',legs{2,ib}],leg,[-Rmax,Rmax],[-Rmax,Rmax],style1,marker,'k',width,siz,gca,red,lspace,bspace);
0152         print_jd(p_opt,[filename,'_toroidal_trajectory_',legs{2,ib}]);
0153     end
0154 end
0155 %
0156 xlim = NaN;
0157 if propvar == 0,
0158     propfield = 'ss';
0159     proplab = 'ray length (m)';
0160     propfac = 1;
0161 else
0162     propfield = 'sphi';
0163     proplab = '|\phi|/(2\pi)';
0164     propfac = 2*pi;
0165 end
0166 %
0167 % --- radial position ---
0168 %
0169 figure('Name','Normalized radial position')
0170 ylim = [0,1];
0171 for ib = 1:nb,
0172     ieq = min([ib,neq]);%equil_fit number
0173     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.srho,0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0174 end
0175 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.srho,0,0,proplab,'\rho','Normalized radial position',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0176 print_jd(p_opt,[filename,'_RhoComp']);
0177 %
0178 % --- poloidal position ---
0179 %
0180 figure('Name','Normalized poloidal position')
0181 ylim = NaN;
0182 for ib = 1:nb,
0183     ieq = min([ib,neq]);%equil_fit number
0184     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.stheta/(2*pi),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0185 end
0186 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.stheta/(2*pi),0,0,proplab,'\theta/(2\pi)','Normalized poloidal position',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0187 print_jd(p_opt,[filename,'_ThetaComp']);
0188 %
0189 % --- toroidal position ---
0190 %
0191 figure('Name','Normalized toroidal position')
0192 ylim = NaN;
0193 for ib = 1:nb,
0194     ieq = min([ib,neq]);%equil_fit number
0195     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.sphi/(2*pi),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0196 end
0197 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.sphi/(2*pi),0,0,proplab,'\phi/(2\pi)','Normalized toroidal position',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0198 print_jd(p_opt,[filename,'_PhiComp']);
0199 %
0200 % --- parallel wave refractive index ---
0201 %
0202 figure('Name','Parallel wave refractive index')
0203 ylim = NaN;
0204 for ib = 1:nb,
0205     ieq = min([ib,neq]);%equil_fit number
0206     graph1D_jd(abs(rays{ib}.(propfield))/propfac,real(rays{ib}.sNpar),0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0207 end
0208 graph1D_jd(abs(rays{nb}.(propfield))/propfac,real(rays{nb}.sNpar),0,0,proplab,'N_{||}','Parallel wave refractive index',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0209 print_jd(p_opt,[filename,'_NparComp']);
0210 %
0211 % --- perpendicular wave refractive index ---
0212 %
0213 figure('Name','Perpendicular wave refractive index')
0214 ylim = NaN;
0215 for ib = 1:nb,
0216     ieq = min([ib,neq]);%equil_fit number
0217     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.sNperp,0,0,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0218 end
0219 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.sNperp,0,0,proplab,'N_{\perp}','Perpendicular wave refractive index',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0220 print_jd(p_opt,[filename,'_NperpComp']);
0221 %
0222 % --- dispersion relation ---
0223 %
0224 figure('Name','Dispersion relation')
0225 ylim = NaN;
0226 for ib = 1:nb,
0227     ieq = min([ib,neq]);%equil_fit number
0228     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.sD,0,1,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0229 end
0230 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.sD,0,1,proplab,'det(D)','Dispersion relation',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0231 print_jd(p_opt,[filename,'_DispComp']);
0232 %
0233 % --- cold plasma mode ---
0234 %
0235 figure('Name','Cold plasma mode')
0236 ylim = NaN;
0237 for ib = 1:nb,
0238     ieq = min([ib,neq]);%equil_fit number
0239     graph1D_jd(abs(rays{ib}.(propfield))/propfac,rays{ib}.spmode,0,1,'','','',NaN,xlim,ylim,style,marker,colors{rem(ib,ncol)},width2);
0240 end
0241 graph1D_jd(abs(rays{nb}.(propfield))/propfac,rays{nb}.spmode,0,1,proplab,'p mode','Cold plasma mode',NaN,xlim,ylim,style,marker,colors{rem(nb,ncol)},width2,siz,gca,red,lspace,bspace);
0242 print_jd(p_opt,[filename,'_PmodeComp']);
0243 %