LH_acces_yp

   Plot of a ray in the contour plot of the n Stix-Golant criterion
    
    Input:

      - equil: equilibrium structure
      - ray: ray structure
     - mask: mask for linear absorption
     - iray: ray index

    Output:

         - na: accessible refractive index on the (psi,theta) grid [npsi,mtheta]

by Y.PEYSSON (yves.peysson@cea.fr) and J. Decker (joan.decker@epfl.ch)

0001 function [na] = LH_acces_yp(equil,ray,mask,iray) 
0002 %
0003 %   Plot of a ray in the contour plot of the n Stix-Golant criterion
0004 %
0005 %    Input:
0006 %
0007 %      - equil: equilibrium structure
0008 %      - ray: ray structure
0009 %     - mask: mask for linear absorption
0010 %     - iray: ray index
0011 %
0012 %    Output:
0013 %
0014 %         - na: accessible refractive index on the (psi,theta) grid [npsi,mtheta]
0015 %
0016 %by Y.PEYSSON (yves.peysson@cea.fr) and J. Decker (joan.decker@epfl.ch)
0017 %
0018 if nargin < 2, error('Not enough input arguments in LH_acces_yp');end
0019 if nargin < 3,
0020     mask = ones(size(ray.sx));
0021     iray = NaN;  
0022 end
0023 %
0024 [na,E,R,Z,psin,B] = stixgolant_yp(equil,ray.omega_rf);
0025 %
0026 if ~isfield(equil,'consistency'),
0027     equil = equilconsistency_yp(equil,'',0);
0028 end
0029 %
0030 [qe,me,mp,mn,e0] = pc_dke_yp;
0031 %
0032 Bfield = sqrt(equil.ptBx.^2 + equil.ptBy.^2 + equil.ptBPHI.^2);
0033 %
0034 ns = [equil.pne;equil.pzni];
0035 qs = qe*[-1;equil.zZi(:)];    
0036 %
0037 ms = [me;mp*equil.zmi(:)];
0038 %
0039 wps = repmat(sqrt(ns.*(qs*ones(1,length(ns))).^2/e0./ms),[1,1,length(equil.theta)]);
0040 wcs = repmat(qs,[1,length(equil.pne),length(equil.theta)]).*shiftdim(repmat(Bfield,[1,1,length(qs)]),2)./repmat(ms,[1,length(equil.pne),length(equil.theta)]);
0041 %
0042 ps = wps/ray.omega_rf;
0043 ys = wcs/ray.omega_rf;
0044 %
0045 ys(ys.^2 == 1) = ys(ys.^2 == 1) + eps;
0046 %
0047 %   Dielectric tensor elements
0048 %
0049 Kperp = squeeze(1.0 - sum(ps.^2./(1 - ys.^2)));
0050 Kpar = squeeze(1.0 - sum(ps.^2));
0051 Kxy = squeeze(sum(ps.^2.*ys./(1.0 - ys.^2)));
0052 %
0053 ntheta = size(Kpar,2);
0054 xR = [flipud(equil.ptx(:,(ntheta-1)/2));equil.ptx(:,1)]' + equil.Rp;%[HFS:LFS] @ Bmin
0055 xTe = [flipud(equil.pTe(:,1));equil.pTe(:,1)]' + equil.Rp;%[HFS:LFS] @ Bmin
0056 %
0057 epsipar = [flipud(Kpar(:,(ntheta-1)/2));Kpar(:,1)]';%[HFS:LFS] @ Bmin
0058 epsiperp = [flipud(Kperp(:,(ntheta-1)/2));Kperp(:,1)]';%[HFS:LFS] @ Bmin
0059 epsixy = [flipud(Kxy(:,(ntheta-1)/2));Kxy(:,1)]';%[HFS:LFS] @ Bmin
0060 xna = [flipud(na(:,(ntheta-1)/2));na(:,1)]';%[HFS:LFS] @ Bmin
0061 %
0062 npar0 = abs(real(ray.sNpar(1)));%launched N||0
0063 %
0064 nnpar = 10000;
0065 nparmax = 10;%max N|| is 10
0066 M = linspace(1,nparmax,nnpar);
0067 x0 = equil.ptx(:,1)/equil.ptx(end,1); 
0068 xx0 = [flipud(x0(:,1));x0(:,1)]';%[HFS:LFS] @ Bmin
0069 [MM,xx] = meshgrid(M,xx0);
0070 MM = MM';
0071 xx = xx';
0072 qq = interp1(equil.consistency.rho,equil.consistency.q,xx,'spline','extrap');
0073 ap = equil.ptx(end,1);
0074 %
0075 P4 = ones(nnpar,1)*epsiperp;
0076 P2 = (ones(nnpar,1)*(epsipar + epsiperp)).*(MM.*MM - (ones(nnpar,1)*epsiperp)) + (ones(nnpar,1)*epsixy).^2;
0077 P0 = (ones(nnpar,1)*epsipar).*((MM.*MM - P4).^2 - (ones(nnpar,1)*epsixy).^2);
0078 K2 = (-P2 + sqrt(P2.*P2 - 4*P0.*P4))./P4/2 - (MM-npar0).*(MM-npar0).*qq.*qq*(equil.Rp/ap)*(equil.Rp/ap)./xx./xx;%electromagnetic limit
0079 %
0080 MMM = zeros(size(K2));
0081 MMM(K2>0) = MM(K2>0);
0082 xxna = ones(nnpar,1)*xna;
0083 maskna = MMM>xxna;% accessibility condition
0084 MMMM = MMM.*maskna;
0085 maxMMM = max(MMM);
0086 MMM(MMM==0) = MMM(MMM==0) + Inf;
0087 minMMM = min(MMM);
0088 maxMMM(maxMMM == 0) = maxMMM(maxMMM == 0)*NaN;
0089 minMMM(minMMM == 0) = minMMM(minMMM == 0)*NaN;
0090 maxMMMM = max(MMMM);
0091 MMMM(MMMM==0) = MMMM(MMMM==0) + Inf;
0092 minMMMM = min(MMMM);
0093 maxMMMM(maxMMMM == 0) = maxMMMM(maxMMMM == 0)*NaN;
0094 %
0095 boundaryMMMM = [minMMMM',maxMMMM'];
0096 boundaryR = [xR',xR'];
0097 %
0098 figure,
0099 graph1D_jd(boundaryR(:),boundaryMMMM(:),0,0,'R (m)','N_{||}','Z = Z_{p}',NaN,NaN,NaN,'none','o','k',2);
0100 graph1D_jd(xR,6.5./sqrt(xTe),0,0,'R (m)','N_{||}','Z = Z_{p}',NaN,NaN,NaN,'-.','none','r',2);
0101 graph1D_jd(xR,xna,0,0,'R (m)','N_{||}','Z = Z_{p}',NaN,NaN,NaN,'--','none','b',2);
0102 legend('Propagation boundaries','N_{||Landau} = 6.5/\surd{T_e}','N_{||acc.} (Stix-Golant)')
0103 title(['N_{||0} = ',num2str(npar0),' @ B = B_{min}'])
0104 %
0105 if 0,
0106     M = linspace(1,nparmax,nnpar);
0107     x0 = equil.ptx(:)'/equil.ptx(end,1); 
0108     [MM,xx] = meshgrid(M,x0);
0109     MM = MM';
0110     xx = xx';
0111     qq = interp1(equil.consistency.rho,equil.consistency.q,xx,'spline','extrap');
0112     ap = equil.ptx(end,1);
0113     %
0114     P4f = ones(nnpar,1)*Kperp(:)';
0115     P2f = (ones(nnpar,1)*(Kpar(:)' + Kperp(:)')).*(MM.*MM - (ones(nnpar,1)*Kperp(:)')) + (ones(nnpar,1)*Kxy(:)').^2;
0116     P0f = (ones(nnpar,1)*Kpar(:)').*((MM.*MM - P4f).^2 - (ones(nnpar,1)*Kxy(:)').^2);
0117     K2f = (-P2f + sqrt(P2f.*P2f - 4*P0f.*P4f))./P4f/2 - (MM-npar0).*(MM-npar0).*qq.*qq*(equil.Rp/ap)*(equil.Rp/ap)./xx./xx;%electromagnetic limit
0118     %
0119     MMM = zeros(size(K2f));
0120     MMM(K2f>0) = MM(K2f>0);
0121     xxna = ones(nnpar,1)*na(:)';
0122     mask = MMM>xxna;% accessibility condition
0123     MMMM = MMM.*mask;
0124     maxMMM = max(MMM);
0125     MMM(MMM==0) = MMM(MMM==0) + Inf;
0126     minMMM = min(MMM);
0127     maxMMM(maxMMM == 0) = maxMMM(maxMMM == 0)*NaN;
0128     minMMM(minMMM == 0) = minMMM(minMMM == 0)*NaN;
0129     maxMMMM = max(MMMM);
0130     MMMM(MMMM==0) = MMMM(MMMM==0) + Inf;
0131     minMMMM = min(MMMM);
0132     maxMMMM(maxMMMM == 0) = maxMMMM(maxMMMM == 0)*NaN;
0133     %
0134     figure,
0135     plot(equil.ptx(:)'+equil.Rp,minMMMM','b',equil.ptx(:)'+equil.Rp,maxMMMM','b');hold on
0136     plot(equil.ptx(:)'+equil.Rp,6.5./sqrt(equil.pte(:)'),'r-.');hold on
0137     plot(equil.ptx(:)'+equil.Rp,na(:)','g--');hold on
0138     plot(ray.sx+equil.Rp,abs(real(ray.sNpar)),'k.')
0139 
0140 end    
0141 %
0142 %Display mode
0143 %
0144 figure('Name','Npar accessible');
0145 subplot(121),graph1D_jd([flipud(R(:,(length(equil.theta)+1)/2));R(:,1)],[flipud(na(:,(length(equil.theta)+1)/2));na(:,1)],0,0,'R (m)','N_{||acc.}','Z = Z_{p}',NaN,NaN,NaN,'-','none','r',2);
0146 grid on
0147 subplot(122),graph2D_jd(R,Z,real(na'),'R (m)','Z (m)',strrep(equil.id,'_','-'),NaN,NaN,1,NaN,NaN,0);colorbar;hold on;
0148 subplot(122),graph2D_jd(R,Z,psin','R (m)','Z (m)',[strrep(equil.id,'_','-'),', ray #',int2str(iray)],NaN,NaN,0,NaN,NaN,0,'--','k');
0149 axis('equal')
0150 set(gca, 'CLim', [1, max(max(na))]);
0151 c = colorbar;
0152 c.Label.String = 'N_{||acc.}';
0153 hold on;
0154 graph1D_jd(ray.sx(mask)+equil.Rp,ray.sy(mask)+equil.Zp,0,0,'','','',NaN,NaN,NaN,'-','none','r',2);
0155 hold off
0156 %
0157 figure('Name','Energy of Resonant Electrons');
0158 subplot(121),graph1D_jd([flipud(R(:,(length(equil.theta)+1)/2));R(:,1)],[flipud(E(:,(length(equil.theta)+1)/2));E(:,1)],0,0,'R (m)','Resonant fast electrons, E_{c} (keV)','Z = Z_{p}',NaN,NaN,NaN,'-','none','r',2);
0159 grid on
0160 subplot(122),graph2D_jd(R,Z,real(E'),'R (m)','Z (m)',strrep(equil.id,'_','-'),NaN,NaN,1,NaN,NaN,0);colorbar;hold on;
0161 subplot(122),graph2D_jd(R,Z,psin','R (m)','Z (m)',[strrep(equil.id,'_','-'),', ray #',int2str(iray)],NaN,NaN,0,NaN,NaN,0,'--','k');
0162 c = colorbar;
0163 c.Label.String = 'Resonant electrons, E_{c} (keV)';
0164 axis('equal')
0165 hold on;
0166 graph1D_jd(ray.sx(mask)+equil.Rp,ray.sy(mask)+equil.Zp,0,0,'','','',NaN,NaN,NaN,'-','none','r',2);
0167 hold off
0168 %