rayparameters_yp

PURPOSE ^

SYNOPSIS ^

function [x,y,Npar,Nperp,epol_pmz,phi_xyz,D,Te,ne,zTi,zni,Brho,BP,Bz,B,r,salpha,gradrho,Upsilon,pmode,wpe2,wce2] = rayparameters_yp(rayparam,equil_fit,fluct_fit,omega,rho,theta,z,krho,m,kz)

DESCRIPTION ^

 Calculates ray position, Npar and Nperp 

 INPUT

 OUTPUT: 

 By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x,y,Npar,Nperp,epol_pmz,phi_xyz,D,Te,ne,zTi,zni,Brho,BP,Bz,B,r,salpha,gradrho,Upsilon,pmode,wpe2,wce2] = rayparameters_yp(rayparam,equil_fit,fluct_fit,omega,rho,theta,z,krho,m,kz)
0002 %
0003 % Calculates ray position, Npar and Nperp
0004 %
0005 % INPUT
0006 %
0007 % OUTPUT:
0008 %
0009 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0010 %
0011 if isfield(equil_fit,'Te_fit'),Te = ppval(equil_fit.Te_fit.pp_f,rho);else Te = zeros(size(rho));end
0012 if isfield(equil_fit,'ne_fit'),ne = ppval(equil_fit.ne_fit.pp_f,rho);else ne = zeros(size(rho));end
0013 if isfield(equil_fit,'zTi_fit'),zTi = ppval(equil_fit.zTi_fit.pp_f,rho);else zTi = zeros(length(rho),length(equil_fit.zZi));end
0014 if isfield(equil_fit,'zni_fit'),zni = ppval(equil_fit.zni_fit.pp_f,rho);else zni = zeros(length(rho),length(equil_fit.zZi));end
0015 %
0016 Rp = equil_fit.Rp;
0017 Zp = equil_fit.Zp;
0018 %
0019 x_a0 = ppval_yp(equil_fit.x_fit.pp_a0,rho);
0020 x_an = ppval_yp(equil_fit.x_fit.pp_an,rho);
0021 x_bn = ppval_yp(equil_fit.x_fit.pp_bn,rho);
0022 %
0023 y_a0 = ppval_yp(equil_fit.y_fit.pp_a0,rho);
0024 y_an = ppval_yp(equil_fit.y_fit.pp_an,rho);
0025 y_bn = ppval_yp(equil_fit.y_fit.pp_bn,rho);
0026 %
0027 Bx_a0 = ppval(equil_fit.Bx_fit.pp_a0,rho);
0028 Bx_an = ppval(equil_fit.Bx_fit.pp_an,rho);
0029 Bx_bn = ppval(equil_fit.Bx_fit.pp_bn,rho);
0030 %
0031 By_a0 = ppval(equil_fit.By_fit.pp_a0,rho);
0032 By_an = ppval(equil_fit.By_fit.pp_an,rho);
0033 By_bn = ppval(equil_fit.By_fit.pp_bn,rho);
0034 %
0035 Bz_a0 = ppval(equil_fit.Bz_fit.pp_a0,rho);
0036 Bz_an = ppval(equil_fit.Bz_fit.pp_an,rho);
0037 Bz_bn = ppval(equil_fit.Bz_fit.pp_bn,rho);
0038 %
0039 BP_a0 = ppval(equil_fit.BP_fit.pp_a0,rho);
0040 BP_an = ppval(equil_fit.BP_fit.pp_an,rho);
0041 BP_bn = ppval(equil_fit.BP_fit.pp_bn,rho);
0042 %
0043 B_a0 = ppval(equil_fit.B_fit.pp_a0,rho);
0044 B_an = ppval(equil_fit.B_fit.pp_an,rho);
0045 B_bn = ppval(equil_fit.B_fit.pp_bn,rho);
0046 %
0047 r_a0 = ppval(equil_fit.r_fit.pp_a0,rho);
0048 r_an = ppval(equil_fit.r_fit.pp_an,rho);
0049 r_bn = ppval(equil_fit.r_fit.pp_bn,rho);
0050 %
0051 calpha_a0 = ppval(equil_fit.calpha_fit.pp_a0,rho);
0052 calpha_an = ppval(equil_fit.calpha_fit.pp_an,rho);
0053 calpha_bn = ppval(equil_fit.calpha_fit.pp_bn,rho);
0054 %
0055 salpha_a0 = ppval(equil_fit.salpha_fit.pp_a0,rho);
0056 salpha_an = ppval(equil_fit.salpha_fit.pp_an,rho);
0057 salpha_bn = ppval(equil_fit.salpha_fit.pp_bn,rho);
0058 %
0059 gradrho_a0 = ppval(equil_fit.gradrho_fit.pp_a0,rho);
0060 gradrho_an = ppval(equil_fit.gradrho_fit.pp_an,rho);
0061 gradrho_bn = ppval(equil_fit.gradrho_fit.pp_bn,rho);
0062 %
0063 x = zeros(length(rho),1);
0064 y = zeros(length(rho),1);
0065 %
0066 r = zeros(length(rho),1);
0067 %
0068 calpha = zeros(length(rho),1);
0069 salpha = zeros(length(rho),1);
0070 gradrho = zeros(length(rho),1);
0071 %
0072 Bx = zeros(length(rho),1);
0073 By = zeros(length(rho),1);
0074 Bz = zeros(length(rho),1);
0075 %
0076 Brho = zeros(length(rho),1);
0077 BP = zeros(length(rho),1);
0078 B = zeros(length(rho),1);
0079 %
0080 for ii = 1:length(rho),
0081    x(ii) = x_a0(ii) + x_an(ii,:)*cos(equil_fit.x_fit.n*theta(ii)) + x_bn(ii,:)*sin(equil_fit.x_fit.n*theta(ii));
0082    y(ii) = y_a0(ii) + y_an(ii,:)*cos(equil_fit.y_fit.n*theta(ii)) + y_bn(ii,:)*sin(equil_fit.y_fit.n*theta(ii));
0083    %
0084    r(ii) = r_a0(ii) + r_an(ii,:)*cos(equil_fit.r_fit.n*theta(ii)) + r_bn(ii,:)*sin(equil_fit.r_fit.n*theta(ii));
0085    %
0086    calpha(ii) = calpha_a0(ii) + calpha_an(ii,:)*cos(equil_fit.calpha_fit.n*theta(ii)) + calpha_bn(ii,:)*sin(equil_fit.calpha_fit.n*theta(ii));
0087    salpha(ii) = salpha_a0(ii) + salpha_an(ii,:)*cos(equil_fit.salpha_fit.n*theta(ii)) + salpha_bn(ii,:)*sin(equil_fit.salpha_fit.n*theta(ii));
0088    gradrho(ii) = gradrho_a0(ii) + gradrho_an(ii,:)*cos(equil_fit.gradrho_fit.n*theta(ii)) + gradrho_bn(ii,:)*sin(equil_fit.gradrho_fit.n*theta(ii));
0089    %
0090    Bx(ii) = Bx_a0(ii) + Bx_an(ii,:)*cos(equil_fit.Bx_fit.n*theta(ii)) + Bx_bn(ii,:)*sin(equil_fit.Bx_fit.n*theta(ii));
0091    By(ii) = By_a0(ii) + By_an(ii,:)*cos(equil_fit.By_fit.n*theta(ii)) + By_bn(ii,:)*sin(equil_fit.By_fit.n*theta(ii));
0092    Bz(ii) = Bz_a0(ii) + Bz_an(ii,:)*cos(equil_fit.Bz_fit.n*theta(ii)) + Bz_bn(ii,:)*sin(equil_fit.Bz_fit.n*theta(ii));
0093    %
0094    BP(ii) = BP_a0(ii) + BP_an(ii,:)*cos(equil_fit.BP_fit.n*theta(ii)) + BP_bn(ii,:)*sin(equil_fit.BP_fit.n*theta(ii));
0095    B(ii) = B_a0(ii) + B_an(ii,:)*cos(equil_fit.B_fit.n*theta(ii)) + B_bn(ii,:)*sin(equil_fit.B_fit.n*theta(ii));
0096 end
0097 %
0098 if isinf(equil_fit.Rp),
0099     Upsilon = 1.0;
0100 else
0101     Upsilon = 1 + x/equil_fit.Rp;%toroidal factor (1 -> cylindrical limit)
0102 end
0103 %
0104 Bs = BP;%Bs is always positive by definition
0105 %
0106 if exist('fluct_fit'),
0107     if isstruct(fluct_fit),
0108         for ii = 1:length(rho),
0109             [fluctval] = fluctval_yp(equil_fit,fluct_fit,rho(ii),theta(ii),phi(ii)); 
0110             %
0111             ne(ii) = ne(ii) + ne(ii)*fluctval.ner;
0112             Te(ii) = Te(ii) + Te(ii)*fluctval.Ter;
0113             zni(ii,:) = zni(ii,:) + zni(ii,:).*fluctval.znir;
0114             zTi(ii,:) = zTi(ii,:) + zTi(ii,:).*fluctval.zTir;
0115             %
0116             cta = cos(theta(ii)).*calpha(ii) + sin(theta(ii)).*salpha(ii);
0117             sta = -cos(theta(ii)).*salpha(ii) + sin(theta(ii)).*calpha(ii);
0118             %
0119             Brho(ii) = Brho(ii) + B(ii)*(cta*fluctval.Bxr + sta*fluctval.Byr);%normalized to unperturbed magnetic equilibrium
0120             Bs(ii) = Bs(ii) + B(ii)*(cta*fluctval.Byr - sta*fluctval.Bxr);%normalized to unperturbed magnetic equilibrium
0121             Bz(ii) = Bz(ii) + B(ii)*fluctval.Bzr;%normalized to unperturbed magnetic equilibrium
0122             %
0123             B(ii) = sqrt(Brho(ii).^2 + Bs(ii).^2 + Bz(ii).^2);%perturbed magnetic field
0124             %
0125         end
0126        %
0127     end
0128 end
0129 %
0130 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0131 %
0132 Npar = clum*(krho.*gradrho.*Brho./B + m.*(calpha.*Bs./B - salpha.*Brho./B)./r + kz.*Bz./Upsilon./B)./omega;
0133 Nperp = sqrt(clum*clum*(krho.*krho.*gradrho.*gradrho - 2*krho.*m.*gradrho.*salpha./r + (m./r).^2 + (kz./Upsilon).^2)./omega./omega - Npar.^2);
0134 %
0135 D = zeros(length(rho),1);
0136 pmode = zeros(length(rho),1);
0137 epol_pmz = zeros(length(rho),3);
0138 phi_xyz = zeros(length(rho),3);
0139 wpe2 = zeros(length(rho),1);
0140 wce2 = zeros(length(rho),1);
0141 %
0142 for ii = 1:length(rho),
0143     if rayparam.tensortype == 0,
0144         [dummy,dummy,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,ps,ys]...
0145             = colddisp_dke_jd(Npar(ii),omega(ii),ne(ii),zni(ii,:)',equil_fit.zZi,equil_fit.zmi,B(ii),Nperp(ii));
0146 %
0147         [Nperpp,Nperpm] = colddisp_dke_jd(Npar(ii),omega(ii),ne(ii),zni(ii,:)',equil_fit.zZi,equil_fit.zmi,B(ii));
0148 %
0149     else
0150         error('not yet implemented')
0151     end
0152 %
0153     pmode(ii) = abs(Nperpp - Nperp(ii)) <= abs(Nperpm - Nperp(ii)); 
0154 %
0155     N = [Nperp(ii);0;Npar(ii)];
0156     N2 = Npar(ii)^2+Nperp(ii)^2;
0157     D(ii) = det(N*N.'/N2- eye(3) + Kxyz/N2);
0158 %
0159     epol_pmz(ii,:) = ep_pmz.';
0160     phi_xyz(ii,:) = phip_xyz.';
0161     %
0162     wpe2(ii) = ps(1)^2;
0163     wce2(ii) = ys(1)^2;
0164 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.