stixgolant_yp

    Calculation of the Stix-Golant LH ray accessibility
    
    Input:

      - equil: equilibrium structure
      - omaga_lh : lh wave frequency = 2*pi*f (Hz)

    Output:

         - na: accessible refractive index on the (psi,theta) grid [npsi,mtheta]
         - E: resonant electron energies at na (psi,theta) grid [npsi,mtheta]
         - R: Major radius (psi,theta) grid [npsi,mtheta]
         - Z: Vertical position (psi,theta) grid [npsi,mtheta]
         - psin: %-normalized poloidal magnetic flux (psi,theta) grid [npsi,mtheta]
         - B: magnetic field (psi,theta) grid [npsi,mtheta]


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

0001 function [na,E,R,Z,psin,B] = stixgolant_yp(equil,omega_lh)
0002 %
0003 %    Calculation of the Stix-Golant LH ray accessibility
0004 %
0005 %    Input:
0006 %
0007 %      - equil: equilibrium structure
0008 %      - omaga_lh : lh wave frequency = 2*pi*f (Hz)
0009 %
0010 %    Output:
0011 %
0012 %         - na: accessible refractive index on the (psi,theta) grid [npsi,mtheta]
0013 %         - E: resonant electron energies at na (psi,theta) grid [npsi,mtheta]
0014 %         - R: Major radius (psi,theta) grid [npsi,mtheta]
0015 %         - Z: Vertical position (psi,theta) grid [npsi,mtheta]
0016 %         - psin: %-normalized poloidal magnetic flux (psi,theta) grid [npsi,mtheta]
0017 %         - B: magnetic field (psi,theta) grid [npsi,mtheta]
0018 %
0019 %
0020 %by Y.PEYSSON (yves.peysson@cea.fr) and J. Decker (joan.decker@epfl.ch)
0021 %
0022 %
0023 [qe,me,mp,~,e0] = pc_dke_yp;%Universal physics constant
0024 %
0025 R = equil.Rp + equil.ptx;
0026 Z = equil.Zp + equil.pty;
0027 B = sqrt(equil.ptBPHI.^2 + equil.ptBx.^2 + equil.ptBy.^2);%magnetic field
0028 psin = equil.psi_apRp(:)*ones(1,length(equil.theta))/equil.psi_apRp(end);%normalized poloidal magnetic flux
0029 %
0030 %Calculation of epsiperp
0031 %
0032 omega_cen = (qe/me).*B/omega_lh;
0033 omega_pen = qe*sqrt(equil.pne(:)*ones(1,length(equil.theta))/e0/me)/omega_lh;
0034 omega_cin = shiftdim(repmat((qe*equil.zZi./mp./equil.zmi)',1,length(equil.psi_apRp),length(equil.theta)),1).*repmat(B,1,1,length(equil.zZi))/omega_lh;
0035 omega_pin = shiftdim(repmat((qe*equil.zZi./sqrt(e0.*mp.*equil.zmi))',1,length(equil.psi_apRp),length(equil.theta)),1).*shiftdim(repmat(sqrt(equil.pzni),1,1,length(equil.theta)),1)/omega_lh;
0036 %
0037 Xperp = -omega_pen.^2./(1.0 - omega_cen.^2) - sum(omega_pin.^2./(1.0 - omega_cin.^2),3);
0038 Xpar = -omega_pen.^2 - sum(omega_pin.^2,3);
0039 Xx = omega_pen.^2.*omega_cen./(1.0 - omega_cen.^2) + sum(omega_pin.^2.*omega_cin./(1.0 - omega_cin.^2),3);
0040 %
0041 Kperp = 1 + Xperp;
0042 Kpar = 1 + Xpar;
0043 Kx = Xx;
0044 %
0045 % Cold fast-slow wave conversion when P2^2 - 4P0P4 = 0 sets the Stix-Golant LH ray accessibility condition
0046 %
0047 a = (Kpar + Kperp).^2 - 4*Kpar.*Kperp;
0048 b = 8*Kpar.*Kperp.^2 - 2*Kperp.*(Kpar + Kperp).^2 + 2*(Kpar + Kperp).*Kx.^2; 
0049 c = -4*Kpar.*Kperp.^3 + Kperp.^2.*(Kpar + Kperp).^2 + 4*Kpar.*Kperp.*Kx.^2 - 2*Kperp.*(Kpar + Kperp).*Kx.^2 + Kx.^4;
0050 delta = b.^2 - 4*a.*c;
0051 %
0052 na = sqrt((-b + sqrt(delta))./2./a);%N||acc
0053 %
0054 E = 511*(1.0./sqrt(1-1.0./na.^2)-1);%Energy of the resonant electrons