init_RT_jd

PURPOSE ^

SYNOPSIS ^

function xY = init_RT_jd(equil_fit,xrho,xtheta,xz,xomega,xNpar,xm,xkz,mmode,kmode,fluct_fit);

DESCRIPTION ^

 This function solves the dispersion relation from initial positions and 
 wave characteristics to initiate ray-tracing calculations.

   INPUTS :

       - equil_fit: fitted axisymmetric magnetic equilibrium structure
       - xrho: normalized minor radius positions [1,1]
       - xtheta: poloidal positions (rad) [1,1]
       - xz: axial or toroidal positions [m] [1,1]
       - xomega: wave frequencies (rad/s) [1,1]
       - xNpar: parallel indices of refraction [1,1]
       - xm: poloidal mode numbers [1,1]
       - xkz: axial or toroidal wave vector [1,1]
       - mmode: cold plasma mode [1,1] : (-1) m (1) p 
           p is the slow mode when kperp > 0 (ex : LH slow wave)
           if mmode has an imaginary part, it is outward propagation, otherwise inward
       - sigma: direction of the wave propagation (outward, or inward)
       - kmode: kinetic (1) or cold (0) plasma model
       - fluct_fit: fluctuations of the magnetic equilibrium

   OUTPUTS :

       - xY: canonical coordinates vector [8,nx] decomposed into
               - xpsin [1,1]
               - xtheta (rad) [1,1]
               - xz (rad) [1,1]
               - xkrhon  [1,1]
               - xm [1,1]
               - xkz [1,1]
               - xomega/clum (m^-1) [1,1]

   By J. Decker DRFC-CEA <joan.decker@cea.fr> and Y. Peysson DRFC-CEA <yves.peysson@cea.fr>

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function xY = init_RT_jd(equil_fit,xrho,xtheta,xz,xomega,xNpar,xm,xkz,mmode,kmode,fluct_fit);
0002 %
0003 % This function solves the dispersion relation from initial positions and
0004 % wave characteristics to initiate ray-tracing calculations.
0005 %
0006 %   INPUTS :
0007 %
0008 %       - equil_fit: fitted axisymmetric magnetic equilibrium structure
0009 %       - xrho: normalized minor radius positions [1,1]
0010 %       - xtheta: poloidal positions (rad) [1,1]
0011 %       - xz: axial or toroidal positions [m] [1,1]
0012 %       - xomega: wave frequencies (rad/s) [1,1]
0013 %       - xNpar: parallel indices of refraction [1,1]
0014 %       - xm: poloidal mode numbers [1,1]
0015 %       - xkz: axial or toroidal wave vector [1,1]
0016 %       - mmode: cold plasma mode [1,1] : (-1) m (1) p
0017 %           p is the slow mode when kperp > 0 (ex : LH slow wave)
0018 %           if mmode has an imaginary part, it is outward propagation, otherwise inward
0019 %       - sigma: direction of the wave propagation (outward, or inward)
0020 %       - kmode: kinetic (1) or cold (0) plasma model
0021 %       - fluct_fit: fluctuations of the magnetic equilibrium
0022 %
0023 %   OUTPUTS :
0024 %
0025 %       - xY: canonical coordinates vector [8,nx] decomposed into
0026 %               - xpsin [1,1]
0027 %               - xtheta (rad) [1,1]
0028 %               - xz (rad) [1,1]
0029 %               - xkrhon  [1,1]
0030 %               - xm [1,1]
0031 %               - xkz [1,1]
0032 %               - xomega/clum (m^-1) [1,1]
0033 %
0034 %   By J. Decker DRFC-CEA <joan.decker@cea.fr> and Y. Peysson DRFC-CEA <yves.peysson@cea.fr>
0035 %
0036 if nargin < 9,
0037     error('Not enough input arguments')
0038 end
0039 %
0040 if nargin < 10,
0041     kmode = 0;
0042 end
0043 %
0044 if nargin < 11,
0045     fluct_fit = [];
0046 end
0047 %
0048 %if abs(sigma) ~= 1,error('Parameter ''sigma'' for the wave direction atlaunch must be +1 or +1');end
0049 %
0050 if length(xtheta) ~= 1 | length(xz) ~= 1 | length(xomega) ~= 1 | length(xNpar) ~= 1 | length(xm) ~= 1 | length(xkz) ~= 1,,
0051     error('Input argument dimensions are not compatible');
0052 end
0053 %
0054 zmi = equil_fit.zmi;
0055 zZi = equil_fit.zZi;
0056 %
0057 if isfield(equil_fit,'Te_fit'),xTe = ppval(equil_fit.Te_fit.pp_f,xrho);else xTe = zeros(size(xrho));end
0058 if isfield(equil_fit,'ne_fit'),xne = ppval(equil_fit.ne_fit.pp_f,xrho);else xne = zeros(size(xrho));end
0059 if isfield(equil_fit,'zTi_fit'),zxTi = ppval(equil_fit.zTi_fit.pp_f,xrho);else zxTi = zeros(length(zZi),length(xrho));end
0060 if isfield(equil_fit,'zni_fit'),zxni = ppval(equil_fit.zni_fit.pp_f,xrho);else zxni = zeros(length(zZi),length(xrho));end
0061 %
0062 x_a0_fit = ppval(equil_fit.x_fit.pp_a0,xrho);
0063 x_an_fit = ppval(equil_fit.x_fit.pp_an,xrho);
0064 x_bn_fit = ppval(equil_fit.x_fit.pp_bn,xrho);
0065 %
0066 y_a0_fit = ppval(equil_fit.y_fit.pp_a0,xrho);
0067 y_an_fit = ppval(equil_fit.y_fit.pp_an,xrho);
0068 y_bn_fit = ppval(equil_fit.y_fit.pp_bn,xrho);
0069 %
0070 r_a0_fit = ppval(equil_fit.r_fit.pp_a0,xrho);
0071 r_an_fit = ppval(equil_fit.r_fit.pp_an,xrho);
0072 r_bn_fit = ppval(equil_fit.r_fit.pp_bn,xrho);
0073 %
0074 c_a0_fit = ppval(equil_fit.calpha_fit.pp_a0,xrho);
0075 c_an_fit = ppval(equil_fit.calpha_fit.pp_an,xrho);
0076 c_bn_fit = ppval(equil_fit.calpha_fit.pp_bn,xrho);
0077 %
0078 s_a0_fit = ppval(equil_fit.salpha_fit.pp_a0,xrho);
0079 s_an_fit = ppval(equil_fit.salpha_fit.pp_an,xrho);
0080 s_bn_fit = ppval(equil_fit.salpha_fit.pp_bn,xrho);
0081 %
0082 Bx_a0_fit = ppval(equil_fit.Bx_fit.pp_a0,xrho);
0083 Bx_an_fit = ppval(equil_fit.Bx_fit.pp_an,xrho);
0084 Bx_bn_fit = ppval(equil_fit.Bx_fit.pp_bn,xrho);
0085 %
0086 By_a0_fit = ppval(equil_fit.By_fit.pp_a0,xrho);
0087 By_an_fit = ppval(equil_fit.By_fit.pp_an,xrho);
0088 By_bn_fit = ppval(equil_fit.By_fit.pp_bn,xrho);
0089 %
0090 Bz_a0_fit = ppval(equil_fit.Bz_fit.pp_a0,xrho);
0091 Bz_an_fit = ppval(equil_fit.Bz_fit.pp_an,xrho);
0092 Bz_bn_fit = ppval(equil_fit.Bz_fit.pp_bn,xrho);
0093 %
0094 BP_a0_fit = ppval(equil_fit.BP_fit.pp_a0,xrho);
0095 BP_an_fit = ppval(equil_fit.BP_fit.pp_an,xrho);
0096 BP_bn_fit = ppval(equil_fit.BP_fit.pp_bn,xrho);
0097 %
0098 B_a0_fit = ppval(equil_fit.B_fit.pp_a0,xrho);
0099 B_an_fit = ppval(equil_fit.B_fit.pp_an,xrho);
0100 B_bn_fit = ppval(equil_fit.B_fit.pp_bn,xrho);
0101 %
0102 gradrho_a0_fit = ppval(equil_fit.gradrho_fit.pp_a0,xrho);
0103 gradrho_an_fit = ppval(equil_fit.gradrho_fit.pp_an,xrho);
0104 gradrho_bn_fit = ppval(equil_fit.gradrho_fit.pp_bn,xrho);
0105 %
0106 % Build interpolated magnetic fields
0107 %
0108 xx = calcval_yp(equil_fit,xtheta,x_a0_fit,x_an_fit,x_bn_fit);
0109 xy = calcval_yp(equil_fit,xtheta,y_a0_fit,y_an_fit,y_bn_fit);
0110 xr = calcval_yp(equil_fit,xtheta,r_a0_fit,r_an_fit,r_bn_fit);
0111 xs = calcval_yp(equil_fit,xtheta,s_a0_fit,s_an_fit,s_bn_fit);
0112 xc = calcval_yp(equil_fit,xtheta,c_a0_fit,c_an_fit,c_bn_fit);
0113 xgradrho = calcval_yp(equil_fit,xtheta,gradrho_a0_fit,gradrho_an_fit,gradrho_bn_fit);
0114 xBx = calcval_yp(equil_fit,xtheta,Bx_a0_fit,Bx_an_fit,Bx_bn_fit);
0115 xBy = calcval_yp(equil_fit,xtheta,By_a0_fit,By_an_fit,By_bn_fit);
0116 xBz = calcval_yp(equil_fit,xtheta,Bz_a0_fit,Bz_an_fit,Bz_bn_fit);
0117 xBP = calcval_yp(equil_fit,xtheta,BP_a0_fit,BP_an_fit,BP_bn_fit);
0118 xB = calcval_yp(equil_fit,xtheta,B_a0_fit,B_an_fit,B_bn_fit);
0119 %
0120 if isinf(equil_fit.Rp),
0121     Upsilon = 1.0;
0122 else
0123     Upsilon = 1 + xx/equil_fit.Rp;%toroidal factor (1 -> cylindrical limit)
0124 end
0125 %
0126 xBsn = xBP./xB;%xBsn is always positive by definition
0127 xBzn = xBz./xB;
0128 %
0129 % Plasma fluctuations or modulations
0130 %
0131 xBRHOn = 0;
0132 %
0133 if isstruct(fluct_fit),
0134    [fluctval] = fluctval_yp(equil_fit,fluct_fit,xrho,xtheta,xz/equil_fit.Rp);
0135    %
0136    xne = xne + xne*fluctval.ner;%tanh(x) is for limitating the perturbation between -1 and 1
0137    xTe = xTe + xTe*fluctval.Ter;
0138    zxni = zxni + zxni.*fluctval.znir;
0139    zxTi = zxTi + zxTi.*fluctval.zTir;
0140    %
0141    xcta = cos(xtheta).*xc + sin(xtheta).*xs;
0142    xsta = -cos(xtheta).*xs + sin(xtheta).*xc;
0143    %
0144    xBRHO = xB*(xcta*fluctval.Bxr + xsta*fluctval.Byr);%normalized to unperturbed magnetic equilibrium
0145    xBs = xB*(xBsn + xcta*fluctval.Byr - xsta*fluctval.Bxr);%normalized to unperturbed magnetic equilibrium
0146    xBz = xB*(xBzn + fluctval.Bzr);%normalized to unperturbed magnetic equilibrium
0147    %
0148    xB = sqrt(xBRHO.^2 + xBs.^2 + xBz.^2);%perturbed magnetic field
0149    %
0150    xBRHOn = xBRHO/xB;%normalized to perturbed magnetic equilibrium
0151    xBsn = xBs/xB;%normalized to perturbed magnetic equilibrium
0152    xBzn = xBz/xB;%normalized to perturbed magnetic equilibrium
0153    %
0154 end
0155 %
0156 % Resolution of the dispersion relation
0157 %
0158 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0159 %
0160 sm = [me;mp*zmi(:)];
0161 sZ = qe*[-1;zZi(:)];
0162 sxn = [xne;zxni(:)];
0163 sxT = [xTe;zxTi(:)];
0164 %
0165 xNperp = zeros(1,1);
0166 %
0167 % Case with complex Npar, in which case Ns = real(Npar) & Nz = imag(Npar)
0168 %
0169 if ~isreal(xNpar),
0170     xNs = real(xNpar);
0171     xNz = imag(xNpar);
0172     %
0173     xm = xr.*xNs.*xomega/clum./xc;
0174     xkz = Upsilon.*xNz.*xomega/clum;
0175     xNpar = xBsn.*xNs + xBzn.*xNz;
0176 end
0177 %
0178 % Cold plasma calculations
0179 %
0180 [Nperpp,Nperpm,dummy,...
0181     dummy,dummy,dummy,dummy,dummy,dummy,...
0182     phip_xyz,phim_xyz] = colddisp_dke_jd(xNpar,xomega,xne,zxni,zZi,zmi,xB);
0183 %
0184 if real(mmode) == 1,
0185     Nperpc = Nperpm;
0186     phic_xyz = phim_xyz;
0187 elseif real(mmode) == -1,
0188     Nperpc = Nperpp;
0189     phic_xyz = phip_xyz;
0190 else
0191     error('Wrong input in mmode')
0192 end
0193 %
0194 % Sign for choosing the inward propagating root (mmode real)
0195 %
0196 sigma = -sign(phic_xyz(1));
0197 %
0198 if ~isreal(mmode)
0199     sigma = - sigma;
0200 end
0201 %
0202 % Normalization
0203 %
0204 if kmode == 2,
0205     sa = sZ.^2.*sxn./sm/e0/xomega^2;
0206     sy = sZ*xB./sm/xomega;
0207     sbeta = sqrt(qe*sxT./sm/clum^2);
0208     %
0209     xNperp = kineticdisp_jd(sa,sy,sbeta,xNpar,Nperpc,0);%   Kinetic plasma calculations
0210 else
0211     xNperp = Nperpc;
0212 end
0213 %
0214 xlambda = xc.*xBsn - xs.*xBRHOn;
0215 %
0216 if abs(xBRHOn) < 1e-6,
0217     %
0218     if isnan(xkz),
0219         xkz = Upsilon.*(xNpar.*xomega/clum - xm.*xc.*xBsn./xr)./xBzn;%xm is specified at the plasma edge (tokamak case)
0220     elseif isnan(xm),
0221         xm = xr.*(xNpar.*xomega/clum - xkz.*xBzn./Upsilon)./xc./xBsn;%xkz is specified at the plasma edge (RFP case)
0222     end
0223     %
0224     xkrho = (xm.*xs./xr + sigma*sqrt(-(xm.*xc./xr).^2 - (xkz./Upsilon).^2 + (xNperp.^2 + xNpar.^2).*xomega.^2/clum^2))./xgradrho;
0225     %
0226 else
0227     %
0228     if isnan(xkz),%xm is specified at the plasma edge (tokamak case)
0229         %
0230         xdelta = xNpar.*xomega./clum - xm.*xlambda./xr;
0231         a = (1.0 - xBsn.*xBsn)./xBRHOn./xBRHOn;
0232         b = xBzn.*(xNpar.*xomega/clum - xm.*xBsn.*xc./xr)./xBRHOn./xBRHOn;
0233         c = (xNpar.*xomega/clum - xm.*xlambda./xr).^2*(1.0 - xBRHOn.*xBRHOn)./xBRHOn./xBRHOn + xm.*xm.*(1.0 - xlambda.*xlambda)./xr./xr - 2.0.*xm.*(xNpar.*xomega/clum - xm.*xlambda./xr).*(xs + xBRHOn.*xlambda)./xr./xBRHOn - (xNperp.*xomega/clum).^2;
0234         xkzp = Upsilon.*(b + sqrt(b.*b - a*c))/a;
0235         xkzm = Upsilon.*(b - sqrt(b.*b - a*c))/a;
0236         xkrhop = (-xkzp.*xBzn./Upsilon + xdelta)./xBRHOn./xgradrho;
0237         xkrhom = (-xkzm.*xBzn./Upsilon + xdelta)./xBRHOn./xgradrho;
0238         %
0239         if imag(xkrhop) ~= 0 & imag(xkrhom) ~= 0,
0240             error('The wave cannot propagate in the plasma. Let change the initial conditions.')
0241         else
0242             xkrhoroots = [xkrhop,xkrhom];
0243             xkzroots = [xkzp,xkzm];
0244             xkrho = xkrhoroots(find(xkrhoroots>0));%choice of the root for inward propagation
0245             xkz = xkzroots(find(xkrhoroots>0));
0246         end
0247         %
0248     elseif isnan(xm),%xkz is specified at the plasma edge (RFP case)
0249         %TODO: calculation for RFP
0250         %xdelta = xNpar.*xomega./clum - xkz.*xBzn./Uspilon;
0251         %a = ;
0252         %b = ;
0253         %c = ;
0254         %xmp = xr.*(b + sqrt(b.*b - a*c))/a;
0255         %xmm = xr.*(b - sqrt(b.*b - a*c))/a;
0256         %xkrhop = (-xmp.*xlambda./xr + xdelta)./xBRHOn./xgradrho;
0257         %xkrhom = (-xmm.*xlambda./xr + xdelta)./xBRHOn./xgradrho;
0258         %
0259         if imag(xkrhop) ~= 0 & imag(xkrhom) ~= 0,
0260             error('The wave cannot propagate in the plasma. Let change the initial conditions.')
0261         else   
0262             if imag(xkrhop) == 0,
0263                 xkrho = xkrhop;
0264                 xm = xmp;
0265             else
0266                 xkrho = xkrhom;
0267                 xm = xmm;
0268             end
0269         end
0270     end
0271 end
0272 %
0273 xY = [xrho;xtheta;xz;xkrho;xm;xkz;xomega];
0274 %
0275 % verification
0276 %
0277 if 1,
0278     xNpar2 = (xkrho*xgradrho.*xBRHOn + xm.*xlambda./xr + xkz.*xBzn./Upsilon)*clum./xomega;
0279     %
0280     %xNperp2 = sqrt((xkrho*xgradrho - xs.*xm./xr).^2 + (xBzn.*xc.*xm./xr - xBsn.*xkz./Upsilon).^2)*clum./xomega;
0281     %
0282     xNperp2 = sqrt(xkrho.^2*xgradrho.^2.*(1.0 - xBRHOn.^2) + xm.^2./xr.^2.*(1.0 - xlambda.^2) + xkz.^2./Upsilon.^2.*(1.0 - xBzn.^2) ... 
0283             - 2.0*xkrho.*xgradrho.*xkz./Upsilon.*xBRHOn.*xBzn ...
0284             - 2.0*xm./xr.*(xkrho.*xgradrho.*(xs + xBRHOn.*xlambda) + xlambda.*xkz.*xBzn./Upsilon))*clum./xomega;    
0285     %
0286     [Nperpp2,Nperpm2,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]...
0287         = colddisp_dke_jd(xNpar2,xomega,xne,zxni,zZi,zmi,xB,xNperp2);
0288    %disp(['xNpar,xNperp,xNpar2,xNperp2,Dp,Dm:',num2str(xNpar),',',num2str(xNperp),',',num2str(xNpar2),',',num2str(xNperp2),',',num2str(Dp),',',num2str(Dm)]);
0289 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.