acces_lh

0001 
0002 
0003 e     = 1.602176462000000e-19;
0004 e2    = e*e;
0005 me    = 9.109381880000001e-31;
0006 mp    = 1.672648500000000e-27;
0007 epsi0 = 8.854187817620391e-12;
0008 x     = gene.x;
0009 dens     = [equil.pne;equil.pzni];
0010 Bfield   = prof.bphi;
0011 omega_rf =  wavestructs{1}.launch.omega_rf;
0012 ms       = [me;mp*equil.zmi(:)];
0013 ns       = dens;
0014 qs       = e*[-1;equil.zZi(:)];    
0015 bphimax  = min(squeeze(equi.BPHI(1,:,:))');
0016 bphi     = double(squeeze(equi.BPHI(1,:,:))');
0017 br       = double(squeeze(equi.BR(1,:,:))');
0018 bz       = double(squeeze(equi.BZ(1,:,:))');
0019 bphi     = sqrt(bphi.^2 + br.^2 + bz.^2);
0020 bphi     = bphi';
0021 R        = double(squeeze(equi.R));
0022 Z        = double(squeeze(equi.Z));
0023 %
0024 
0025 for k=1:length(x)
0026   wps(:,k) = sqrt(ns(:,k).*qs.^2/epsi0./ms);
0027   wcs(:,k) = qs*Bfield(k)./ms;
0028   wcsmax(:,k) = qs*bphimax(k)./ms;
0029   
0030   ps(:,k)  = wps(:,k)/omega_rf(1);
0031   ys(:,k)  = wcs(:,k)/omega_rf(1);
0032   val      = ys(:,k);
0033   val(val.^2 == 1) = val(val.^2 == 1) + eps;
0034   ys(:,k)  = val;
0035   ysmax(:,k)  = wcsmax(:,k)/omega_rf(1);
0036   val      = ysmax(:,k);
0037   val(val.^2 == 1) = val(val.^2 == 1) + eps;
0038   ysmax(:,k)  = val;
0039 %
0040 %   Dielectric tensor elements
0041 %
0042   Kperp = 1 - sum(ps(:,k).^2./(1 - ys(:,k).^2));
0043   Kpar  = 1 - sum(ps(:,k).^2);
0044   Kxy   = sum(ps(:,k).^2.*ys(:,k)./(1 - ys(:,k).^2));
0045   Kperpmax = 1 - sum(ps(:,k).^2./(1 - ysmax(:,k).^2));
0046   Kparmax  = 1 - sum(ps(:,k).^2);
0047   Kxymax   = sum(ps(:,k).^2.*ysmax(:,k)./(1 - ysmax(:,k).^2));
0048 %
0049 %   Accessibility condition
0050 %
0051   Npar(k) = sqrt(Kperp + (sqrt(4*Kperp*Kpar*Kxy^2*(Kxy^2 - (Kperp - Kpar)^2)) - Kxy^2*(Kperp + Kpar))/(Kperp - Kpar)^2);%warning: this solution is valid in the LH frequency range only
0052   Nparmax(k) = sqrt(Kperpmax + (sqrt(4*Kperpmax*Kparmax*Kxymax^2*(Kxymax^2 - (Kperpmax - Kparmax)^2)) - Kxymax^2*(Kperpmax + Kparmax))/(Kperpmax - Kparmax)^2);%warning: this solution is valid in the LH frequency range only
0053 %
0054 % etude en theta
0055 %
0056   for m=1:65
0057      wcstheta(:,k)    = qs*bphi(k,m)./ms;
0058      ystheta(:,k)     = wcstheta(:,k)/omega_rf(1);
0059      val              = ystheta(:,k);
0060      val(val.^2 == 1) = val(val.^2 == 1) + eps;
0061      ystheta(:,k)     = val;
0062 %
0063 %   Dielectric tensor elements
0064 %
0065      Kperpt   = 1 - sum(ps(:,k).^2./(1 - ystheta(:,k).^2));
0066      Kxyt     = sum(ps(:,k).^2.*ystheta(:,k)./(1 - ystheta(:,k).^2));
0067 %
0068 %   Accessibility condition
0069 %
0070      Npart(k,m) = sqrt(Kperpt + (sqrt(4*Kperpt*Kpar*Kxyt^2*(Kxyt^2 - (Kperp - Kpar)^2)) - Kxyt^2*(Kperpt + Kpar))/(Kperpt - Kpar)^2);%warning: this solution is valid in the LH frequency range only
0071    end
0072 end
0073 if ~isstruct(cons) 
0074   npar0 = angle(cons);
0075   Plh   = abs(cons);
0076 else
0077   npar0=angle(cons.lh_cons);
0078   Plh   = abs(cons.lh_cons);  
0079 end
0080 nlh = length(Plh);
0081 acces = 1;
0082 for k=1:nlh
0083   if Plh(k) > 0
0084     if sum(any((-Npart+npar0(k))<0))
0085       acces = acces-1;
0086     end
0087   end
0088 end 
0089 if acces < 1
0090 
0091 
0092   h = findobj(0,'type','figure','tag','acces_lh');
0093   if isempty(h)
0094        h=figure('tag','acces_lh');
0095   else
0096        figure(h);
0097   end   
0098   clf
0099   set(h,'defaultaxesfontsize',12,'defaultaxesfontweight','bold','defaultaxesfontname','times','defaultlinelinewidth',3,'color',[1 1 1])
0100 
0101   col=['r','b','k','m'];
0102   for k=1:nlh    
0103     if Plh(k) > 0
0104       subplot(1,nlh,k)
0105       contourf(R,Z,-Npart+npar0(k),0.01)
0106       colorbar
0107       title(['antenna ',int2str(k),' accessibility boundary, n// > nacc'])
0108       axis(['equal'])
0109     end    
0110   end    
0111   hold on
0112   plot(double(geo.R),double(geo.Z))
0113   hold off
0114 
0115   disp('warning : accessibility problem for the raytracing of LUKE')
0116 
0117 end