RF input parameters for the Drift Kinetic equation solver (dke_4_1yp.m).
Input:
Output:
Warning: none
by J.Decker MIT/RLE <jodecker@mit.edu> and Y.Peysson CEA/DRFC <yves.peysson@cea.fr> 04/29/20030001 function [waveparam,xyprop_dke_rf,xyP0_2piRp_in_rf,xyP0_2piRp_in_rf_coll,xyP0_2piRp_in_rf_noabs,yb_rf,yP0_2piRp,xys_rf] = rfwave_dke_jd(dkeparam,display,equil,equilDKE,gridDKE,Zbouncecoef,mksa,waves); 0002 % 0003 % RF input parameters for the Drift Kinetic equation solver (dke_4_1yp.m). 0004 % 0005 % Input: 0006 % 0007 % Output: 0008 % 0009 % Warning: none 0010 % 0011 %by J.Decker MIT/RLE <jodecker@mit.edu> and Y.Peysson CEA/DRFC <yves.peysson@cea.fr> 04/29/2003 0012 % 0013 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp; 0014 % 0015 %dke_mode = dkeparam.dke_mode; 0016 % 0017 if isfield(dkeparam,'npresmin'), 0018 npresmin = dkeparam.npresmin;% minimum number of resonant points for numerical stability. 0019 else 0020 npresmin = 5; 0021 end 0022 % 0023 display_mode = display.display_mode; 0024 % 0025 xpsin = gridDKE.xpsin; 0026 xpsin_dke = gridDKE.xpsin_dke; 0027 %xdpsin_dke = gridDKE.xdpsin_dke; 0028 %r_dke = gridDKE.rdke; 0029 xpsin_S_dke = gridDKE.xpsin_S_dke; 0030 % 0031 % Xx = equilDKE.Xx; 0032 % Xy = equilDKE.Xy; 0033 % XBx = equilDKE.XBx; 0034 % XBy = equilDKE.XBy; 0035 % XBphi = equilDKE.XBphi; 0036 % xB0 = equilDKE.xB0; 0037 % Rp = equilDKE.Rp; 0038 % ap = equilDKE.ap; 0039 % % 0040 % xqtilde = Zbouncecoef.xqtilde; 0041 % xqhat = Zbouncecoef.xqhat; 0042 % % 0043 % ne_ref = mksa.ne_ref; 0044 % lnc_e_ref = mksa.lnc_e_ref; 0045 % 0046 equil_id = equil.id; 0047 method = 'linear';%'spline';%interpolation technique in psi 0048 % 0049 %nr_dke = length(xpsin_dke); 0050 nr = length(xpsin); 0051 % 0052 % equilibrium data with Fourier-enhanced poloidal precision and original radial grid 0053 % 0054 [equilDKE_Fourier] = equilibrium_jd(equil); 0055 % 0056 % theta = equilDKE_Fourier.mtheta; 0057 % ppsin = equilDKE_Fourier.xpsin; 0058 % ptx = equilDKE_Fourier.Xx; 0059 % pty = equilDKE_Fourier.Xy; 0060 % ptBx = equilDKE_Fourier.XBx; 0061 % ptBy = equilDKE_Fourier.XBy; 0062 % ptBphi = equilDKE_Fourier.XBphi; 0063 % pB0 = equilDKE_Fourier.xB0; 0064 % px0 = equilDKE_Fourier.xx0; 0065 % pBT0 = equilDKE_Fourier.xBT0; 0066 % pBp0 = equilDKE_Fourier.xBp0; 0067 % prho = equilDKE_Fourier.xrho; 0068 % pTe = equilDKE_Fourier.xTe; 0069 % pne = equilDKE_Fourier.xne; 0070 % pzTi = equilDKE_Fourier.xzTi; 0071 % pzni = equilDKE_Fourier.xzni; 0072 % zZi = equilDKE_Fourier.zZi; 0073 % zmi = equilDKE_Fourier.zmi; 0074 % pZeff = equilDKE_Fourier.xZeff; 0075 % Rp = equilDKE_Fourier.Rp; 0076 % Zp = equilDKE_Fourier.Zp; 0077 % % 0078 % ptr = sqrt(ptx.^2 + pty.^2); %minor radius on the (psi,stheta) grid 0079 % ptBP = sqrt(ptBx.^2 + ptBy.^2); %poloidal magnetic field on the (psi,stheta) grid 0080 % ptB = sqrt(ptBP.^2 + ptBphi.^2); %magnetic field on the (psi,stheta) grid 0081 % 0082 ibw = 0; 0083 rays = {}; 0084 % vsrays = {}; 0085 % 0086 n0_rf_list = []; 0087 % 0088 for iw = 1:length(waves),% 0089 % 0090 % ====================================================================== 0091 % treats each wave individually 0092 % ====================================================================== 0093 % 0094 wave = waves{iw}; 0095 % 0096 if wave.equil_id ~= equil_id 0097 error('DKE and WAVE equilibria are different') 0098 end 0099 % 0100 if isfield(wave,'n_rf_list'), 0101 n_rf_list = wave.n_rf_list; 0102 elseif isfield(dkeparam,'n_rf_list'), 0103 n_rf_list = dkeparam.n_rf_list; 0104 elseif isfield(wave,'waveparam') && isfield(wave.waveparam,'n_rf_list'), 0105 n_rf_list = wave.waveparam.n_rf_list; 0106 else 0107 disp('Warning : no n_rf_list data found, enforced to 0:3') 0108 n_rf_list = 0:3; 0109 end 0110 % 0111 if isfield(wave,'model') && strcmp(wave.model,'FW'),%full wave simulation 0112 waveparam = wave; 0113 waveparam.n_rf_list = n_rf_list; 0114 xyprop_dke_rf = zeros(nr,0); 0115 xyP0_2piRp_in_rf = zeros(nr,0); 0116 xys_rf = zeros(nr,0); 0117 yb_rf = zeros(1,0); 0118 yP0_2piRp = zeros(1,0); 0119 % 0120 return 0121 end 0122 % 0123 omega_rf = wave.omega_rf;% The wave frequency (in rad/s) [1,1] 0124 % 0125 if isfield(wave,'dsmin'), 0126 dsmin = wave.dsmin;%minimum interval for ray fragments 0127 else 0128 dsmin = Inf; 0129 end 0130 % 0131 %rays = wave.rays;%list of rays {1,nb} 0132 opt_rf = wave.opt_rf;%option for FLR calculations [1,1] 0133 opt_disp = wave.opt_disp; % option for solving the wave equation, which gives (ray.sepol_pmz, ray.sphi, ray.sNperp) 0134 % as a function of the ray location, the equilibrium and rays_Npar. 0135 %(0): given by ray-tracing code (requires the three parameters above to be given by RT) 0136 %(1): solved in DKE, cold plasma model -colddisp_dke_jd- 0137 %(2): solved in DKE, warm plasma model -warmdisp_dke_jd- (first order thermal corrections) 0138 %(3): solved in DKE, non-relativistic hot plasma model -R2D2- (Warning: so far, electrons only) 0139 %(4): solved in DKE, relativistic hot plasma model -R2D2- (Warning: so far, electrons only) 0140 % 0141 for iray = 1:length(wave.rays),% 0142 % 0143 % ====================================================================== 0144 % treats each ray individually 0145 % ====================================================================== 0146 % 0147 ray = wave.rays{iray}; 0148 % ====================================================================== 0149 % cut invalid ray path !!!!! CRONOS error 0150 % ====================================================================== 0151 % 0152 if ~isempty(ray.sepol_pmz), 0153 ismax = min([find(isnan(ray.sepol_pmz(3,:)))-1,Inf]); 0154 else 0155 ismax = Inf; 0156 end 0157 % 0158 % 0159 % ====================================================================== 0160 % calculates missing ray info if needed 0161 % ====================================================================== 0162 % 0163 % The following information is necessary for including the 0164 % contribution of the ray in DKE: 0165 % - P0: Power launched in the ray (W) [1,1] 0166 % - ss: Propagation distance parameter for each ray [1,ns] 0167 % - spsin: Ray location (normalized magnetic flux) along the propagation [1,ns] 0168 % - stheta: Ray location (poloidal angle) along the propagation [1,ns] 0169 % - sNpar: Ray parallel wave number [1,ns] 0170 % - sNperp: Ray perpendicular wave number [1,ns] 0171 % - sdNpar: Ray spectral width in Npar [1,ns] 0172 % if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar) 0173 % - sepol_pmz: polarization vector [3,ns] complex 0174 % - sphi_xyz: normalized power flow [3,ns] 0175 % 0176 % if these quantities are not calculated by the ray-tracing code, 0177 % they will be determined from the other RT output according to the 0178 % following scheme: 0179 % 0180 % - sR, sphi, sZ -> ss 0181 % - sR, sZ -> spsin, stheta 0182 % - skR, sn, skZ -> sNpar, sNperp 0183 % - sNpar -> sdNpar 0184 % - opt_disp, sNpar, omega_rf, equil -> sepol_pmz, sNperp, sphi_xyz 0185 % 0186 if isfield(wave,'rayparam') && isfield(wave.rayparam,'colldamp'),%backward compatibility 0187 colldamp = wave.rayparam.colldamp; 0188 else 0189 colldamp = 0;%no collisional damping 0190 end 0191 % 0192 ray = rayprocess_jd(ray,equilDKE_Fourier,opt_disp,omega_rf,n_rf_list,method,colldamp); 0193 % 0194 % P0_2piRp = ray.P0_2piRp;% Power launched in the ray (W) [1,1] 0195 % % 0196 % ss = ray.ss(1:min(ismax,length(ray.ss)));% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns] 0197 % spsin = ray.spsin(1:min(ismax,length(ray.spsin)));% Ray location (normalized magnetic flux) along the propagation [1,ns] 0198 % stheta = ray.stheta(1:min(ismax,length(ray.stheta)));% Ray location (poloidal angle) along the propagation [1,ns] 0199 % % 0200 % sNpar = ray.sNpar(1:min(ismax,length(ray.sNpar)));% Ray parallel wave number [1,ns] 0201 % sNperp = ray.sNperp(1:min(ismax,length(ray.sNperp)));% Ray perpendicular wave number [1,ns] 0202 % sdNpar = ray.sdNpar(1:min(ismax,length(ray.sdNpar)));% Ray spectral width in Npar [1,ns] 0203 % % if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar) 0204 % % 0205 % sepol_pmz = ray.sepol_pmz(:,1:min(ismax,size(ray.sepol_pmz,2)));% polarization vector [3,ns] complex 0206 % sphi_xyz = ray.sphi_xyz(:,1:min(ismax,size(ray.sphi_xyz,2)));% normalized power flow [3,ns] 0207 % salphaphi_lin = ray.salphaphi_lin(:,1:min(ismax,length(ray.salphaphi_lin)));% Linear absorption coefficient (m-1) along the propagation [1,ns] 0208 % % 0209 % % Store useful information for DKE in a vsray matrix 0210 % % 0211 % vsray = [ss;spsin;stheta;sNpar;sNperp;sdNpar;sepol_pmz;sphi_xyz;salphaphi_lin]; 0212 % 0213 % Regroup vsray matrices in a list 0214 % 0215 % ibw = ibw + 1; 0216 % vsrays{ibw} = vsray; 0217 % yomega_rf(ibw) = omega_rf; 0218 % yopt_rf(ibw) = opt_rf; 0219 % yP0_2piRp(ibw) = P0_2piRp; 0220 ibw = ibw + 1; 0221 ray.omega_rf = omega_rf; 0222 ray.opt_rf = opt_rf; 0223 ray.dsmin = dsmin;%minimum interval for ray fragments 0224 % 0225 % Minimum spectral width for numerical stability. LH only for now 0226 % 0227 somega_ce = qe*ray.sB/me; 0228 sy = somega_ce./omega_rf; 0229 sNpar = real(ray.sNpar); 0230 sbetath = sqrt(ray.sTe/mc2); 0231 % 0232 snres = zeros(1,length(sy)); 0233 % 0234 snres = snres + (abs(sNpar) < 1 ).*ceil(sqrt(1 - sNpar.^2)./sy); 0235 snres = snres + (abs(sNpar) >= 1).*round(1./sy); 0236 % 0237 syn = snres.*sy; 0238 % 0239 spparmin = abs((abs(sNpar).*syn - sqrt(syn.^2 + sNpar.^2 - 1))./(1 - sNpar.^2))./sbetath; 0240 spparmin1 = NaN(1,length(sy)); 0241 spparmin1(syn ~= 0) = abs(1 - syn(syn ~= 0).^2)./(2*sbetath(syn ~= 0).*syn(syn ~= 0)); 0242 spparmin(sNpar == 1) = spparmin1(sNpar == 1); 0243 % 0244 spparmax = (abs(sNpar).*syn + sqrt(syn.^2 + sNpar.^2 - 1))./(1 - sNpar.^2)./sbetath; 0245 spparmax(sNpar >= 1) = Inf; 0246 % 0247 spparmax(spparmax > spparmin + 1) = spparmin(spparmax > spparmin + 1) + 1; 0248 % 0249 sgammamin = sqrt(1 + sbetath.^2.*spparmin.^2); 0250 sgammamax = sqrt(1 + sbetath.^2.*spparmax.^2); 0251 % 0252 sdNpardpparmin = sbetath.*abs(sNpar).*sqrt(2.*sgammamin.*syn - 1 - sgammamin.^2.*(1 - sNpar.^2))./sgammamin./abs(sgammamin - syn); 0253 sdNpardpparmax = sbetath.*abs(sNpar).*sqrt(2.*sgammamax.*syn - 1 - sgammamax.^2.*(1 - sNpar.^2))./sgammamax./abs(sgammamax - syn); 0254 % 0255 if any(imag([sdNpardpparmin,sdNpardpparmax]) ~= 0), 0256 info_dke_yp(2,['WARNING: imaginary dnparmin found']); 0257 end 0258 % 0259 sdNparmindpparmin = min([sdNpardpparmin;sdNpardpparmax]); 0260 % 0261 sdNparmin = sdNparmindpparmin*npresmin*dkeparam.pnmax_S/dkeparam.np_S; 0262 % 0263 sdNparmin = min([sdNparmin;abs(sNpar)]); 0264 % 0265 smaskmin = imag(ray.sdNpar) == 0 & ray.sdNpar < sdNparmin; 0266 % 0267 if any(smaskmin), 0268 info_dke_yp(2,['WARNING: wave #',num2str(iw),' ray #',num2str(iray),' dNpar increased for ',num2str(sum(smaskmin)),'/',num2str(length(smaskmin)),' data points for numerical stability.']); 0269 info_dke_yp(4,['Average increase factor : ',num2str(mean(sdNparmin(smaskmin)./ray.sdNpar(smaskmin))),'; Max dNpar value, old : ',num2str(max(ray.sdNpar)),' , new : ',num2str(max(sdNparmin(smaskmin)))]); 0270 ray.sdNpar(smaskmin) = sdNparmin(smaskmin); 0271 end 0272 % 0273 rays{ibw} = ray;% minimum number of resonant points for numerical stability. 0274 end 0275 % 0276 n0_rf_list = min([n0_rf_list,n_rf_list]):max([n0_rf_list,n_rf_list]); 0277 % 0278 end 0279 % 0280 % ====================================================================== 0281 % Discretization of the rays on the flux-surfaces 0282 % ====================================================================== 0283 % 0284 [waveparam,xyprop_dke_rf,xys_rf,xyP0_2piRp_in_rf,xyP0_2piRp_in_rf_coll,xyP0_2piRp_in_rf_noabs,yb,yP0_2piRp] = propagation_jd(rays,xpsin_S_dke,xpsin_dke,method,display_mode); 0285 %disp('-------------------------------------------------------------------------------------------------------------------'); 0286 % 0287 if isfield(dkeparam,'n_rf_list'),% in this case the LUKE n_rf_list is imposed by dkeparam 0288 if display.display_mode > 0, 0289 disp(['List of cyclotron harmonics in LUKE imposed by dkeparam : ',num2str(dkeparam.n_rf_list)]) 0290 end 0291 waveparam.n_rf_list = dkeparam.n_rf_list; 0292 else 0293 if display.display_mode > 0, 0294 disp(['List of cyclotron harmonics in LUKE compiled from all waves : ',num2str(n0_rf_list)]) 0295 end 0296 waveparam.n_rf_list = n0_rf_list; 0297 end 0298 % 0299 nb_rf = size(yb,2);%number of actual rays 0300 %nbb_rf = size(xyprop_dke_rf,2);%number of ray subdivisions 0301 % 0302 if nb_rf == 0%no rays 0303 waveparam.model = 'RT'; 0304 waveparam.bomega_rf = ''; 0305 xyprop_dke_rf = zeros(nr,0); 0306 xyP0_2piRp_in_rf = zeros(nr,0); 0307 xys_rf = zeros(nr,0); 0308 yb_rf = zeros(1,0); 0309 % 0310 return 0311 end 0312 % 0313 yb_rf = yb(1,:);%position of ray starts 0314 % for ibb = 1:nbb_rf 0315 % ib_rf(ibb) = max(find(yb_rf <= ibb));%new ray number 0316 % end 0317 % ib_rf_old = ib_old(ib_rf);%old ray number (may be different if some rays did not enter the plasma) 0318 % % 0319 % % parameters that do not vary along the ray 0320 % % 0321 % xyomega_rf = ones(nr,1)*yomega_rf(ib_rf_old);% wave frequency 0322 % xyopt_rf = ones(nr,1)*yopt_rf(ib_rf_old);% FLR option 0323 % % 0324 % % parameters that vary along the ray and are so far calculated for dke positions only 0325 % % 0326 % % Note: for positions that are not visited by the ray, we set xyP0_2piRp_in_rf = 0 0327 % % 0328 % xymask_prop_dke = xyprop_dke_rf ~= 0;%DKE points visited by the ray 0329 % if dke_mode == 1 % interpolation for the FP positions (only necessary if dke_mode == 1) 0330 % % 0331 % xyvray = zeros(nr,nbb_rf,12); 0332 % xyfinc_rf = zeros(nr,nbb_rf); 0333 % xyP0_2piRp_in_rf = zeros(nr,nbb_rf); 0334 % % 0335 % for ibb = 1:nbb_rf 0336 % xmask_prop_fp = [(min(r_dke(xymask_prop_dke(:,ibb))) - 1):(max(r_dke(xymask_prop_dke(:,ibb))) + 1)];%corresponding FP points 0337 % % 0338 % if nr_dke == 1, 0339 % xyvray(xmask_prop_fp,ibb,:) = repmat(xyvray_dke(xymask_prop_dke(:,ibb),ibb,:),[3,1,1]); 0340 % xyfinc_rf(xmask_prop_fp,ibb) = repmat(xyfinc_dke_rf(xymask_prop_dke(:,ibb),ibb),[3,1,1]); 0341 % xyP0_2piRp_in_rf(xmask_prop_fp,ibb) = yP0_2piRp(ib_rf_old(ibb));% initial incident power 0342 % else 0343 % xyvray(xmask_prop_fp,ibb,:) = interp1(xpsin_dke(xymask_prop_dke(:,ibb)),xyvray_dke(xymask_prop_dke(:,ibb)),xpsin(xmask_prop_fp),method); 0344 % xyfinc_rf(xmask_prop_fp,ibb) = interp1(xpsin_dke(xymask_prop_dke(:,ibb)),xyfinc_dke_rf(xymask_prop_dke(:,ibb)),xpsin(xmask_prop_fp),method); 0345 % xyP0_2piRp_in_rf(xmask_prop_fp,ibb) = yP0_2piRp(ib_rf_old(ibb));% initial incident power 0346 % end 0347 % % 0348 % end 0349 % else 0350 % xyvray = xyvray_dke; 0351 % xyfinc_rf = xyfinc_dke_rf; 0352 % xyP0_2piRp_in_rf = xymask_prop_dke.*(ones(nr,1)*yP0_2piRp(ib_rf_old));% initial incident power 0353 % end 0354 % % 0355 % xys_rf = xyvray(:,:,1); 0356 % xythetab_rf = xyvray(:,:,3); 0357 % xyNpar_in = xyvray(:,:,4); 0358 % xyNperp_rf = xyvray(:,:,5); 0359 % xydNpar_in = xyvray(:,:,6); 0360 % xyepolp_rf = xyvray(:,:,7); 0361 % xyepolm_rf = xyvray(:,:,8); 0362 % xyepolz_rf = xyvray(:,:,9); 0363 % xyphix_rf = xyvray(:,:,10); 0364 % xyphiy_rf = xyvray(:,:,11); 0365 % xyphiz_rf = xyvray(:,:,12); 0366 % xyalphaphi_lin_rf = xyvray(:,:,13); 0367 % % 0368 % % Spectral properties 0369 % % 0370 % xyNparmin_rf = NaN*ones(nr,nbb_rf); 0371 % xyNparmax_rf = NaN*ones(nr,nbb_rf); 0372 % xyNpar_rf = NaN*ones(nr,nbb_rf); 0373 % xydNpar_rf = NaN*ones(nr,nbb_rf); 0374 % % 0375 % xymaskGauss_rf = find(imag(xydNpar_in) == 0); 0376 % xymaskSquare_rf = find(imag(xydNpar_in) ~= 0); 0377 % % 0378 % xyNpar_rf(xymaskGauss_rf) = xyNpar_in(xymaskGauss_rf); 0379 % xydNpar_rf(xymaskGauss_rf) = xydNpar_in(xymaskGauss_rf); 0380 % % 0381 % xyNparmin_rf(xymaskSquare_rf) = xyNpar_in(xymaskSquare_rf) - imag(xydNpar_in(xymaskSquare_rf))/2; 0382 % xyNparmax_rf(xymaskSquare_rf) = xyNpar_in(xymaskSquare_rf) + imag(xydNpar_in(xymaskSquare_rf))/2; 0383 % % 0384 % wpe_ref = sqrt(ne_ref*qe^2/e0/me); 0385 % % 0386 % xyabsphi_xyz_rf = sqrt(abs(xyphix_rf).^2 + abs(xyphiy_rf).^2 + abs(xyphiz_rf).^2); 0387 % % 0388 % xyabsphi_xyz_rf(xyabsphi_xyz_rf == 0) = eps;%avoid singularities 0389 % xydNpar_rf(xydNpar_rf == 0) = eps;%avoid singularities 0390 % % 0391 % [xyBb_rf,xyBPb_rf,xyrb_rf,xycosb_rf,xyxb_rf] = mg_interp_jd(theta,Xx,Xy,XBx,XBy,XBphi,xythetab_rf,method); 0392 % % 0393 % xyB0 = repmat(xB0.',[1,nbb_rf]); 0394 % xyqhat = repmat(xqhat.',[1,nbb_rf]); 0395 % xyqtilde = repmat(xqtilde.',[1,nbb_rf]); 0396 % % 0397 % xyFhat = xyrb_rf.*xyB0./(ap*xyBPb_rf.*xyqhat.*xycosb_rf.*(1 + xyxb_rf/Rp)); 0398 % xyFtilde = xyrb_rf.*xyBb_rf./(ap*xyBPb_rf.*xyqtilde.*xycosb_rf); 0399 % % 0400 % xyD0_rf_ref = 2*pi*xyfinc_rf.*xyP0_2piRp_in_rf.*xyFhat.*xyFtilde./(xyrb_rf.*me*lnc_e_ref.*xyomega_rf*wpe_ref^2.*xyabsphi_xyz_rf); 0401 % % 0402 % waveparam = zeros(nr,nbb_rf,13); 0403 % % 0404 % waveparam(:,:,1) = xyomega_rf;%RF wave frequency 0405 % waveparam(:,:,2) = xyopt_rf;%Options for Larmor radius effects and/or EBW calc. 0406 % waveparam(:,:,3) = xythetab_rf;%Poloidal location of RF beam (relevant for bounce_mode > 0 only) 0407 % waveparam(:,:,4) = xyD0_rf_ref;%RF quasi-linear diffusion coefficient (nhuth_ref*pth_ref^2) 0408 % waveparam(:,:,5) = xyNparmin_rf;%Square N// spectrum: Lower limit 0409 % waveparam(:,:,6) = xyNparmax_rf;%Square N// spectrum: Upper limit 0410 % waveparam(:,:,7) = xyNpar_rf;%Gaussian N// spectrum: peak n// 0411 % waveparam(:,:,8) = xydNpar_rf;%Gaussian N// spectrum: width in n// 0412 % waveparam(:,:,9) = xyNperp_rf;%RF wave : n_perp 0413 % waveparam(:,:,10) = xyepolp_rf;%RF wave : polarization plus 0414 % waveparam(:,:,11) = xyepolm_rf;%RF wave : polarization minus 0415 % waveparam(:,:,12) = xyepolz_rf;%RF wave : polarization parallel 0416 % waveparam(:,:,13) = xyalphaphi_lin_rf;%RF wave : Linear absorption coefficient 0417 %