Description of rayprocess

0001 function ray = rayprocess_jd(ray,equilDKE_Fourier,opt_disp,omega_rf,n_rf_list,method,colldamp);
0002 %
0003 % This function calculates missing ray info if needed.
0004 % Also calculates the correct Nperp, polarization and power flow if required
0005 %
0006 % created 03/01/05
0007 % by J. Decker and Y. Peysson
0008 %
0009 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;
0010 %
0011 % equilibrium data with Fourier-enhanced poloidal precision and original radial grid
0012 %
0013 %
0014 if nargin < 7,
0015     colldamp = 0;%no collisional damping
0016 end
0017 %
0018 theta = equilDKE_Fourier.mtheta;
0019 ppsin = equilDKE_Fourier.xpsin;
0020 ptx = equilDKE_Fourier.Xx;
0021 pty = equilDKE_Fourier.Xy;
0022 ptBx = equilDKE_Fourier.XBx;
0023 ptBy = equilDKE_Fourier.XBy;
0024 ptBphi = equilDKE_Fourier.XBphi;
0025 pB0 = equilDKE_Fourier.xB0;
0026 px0 = equilDKE_Fourier.xx0;
0027 pBT0 = equilDKE_Fourier.xBT0;
0028 pBp0 = equilDKE_Fourier.xBp0;
0029 prho = equilDKE_Fourier.xrho;
0030 pTe = equilDKE_Fourier.xTe;
0031 pne = equilDKE_Fourier.xne;
0032 pzTi = equilDKE_Fourier.xzTi;
0033 pzni = equilDKE_Fourier.xzni;
0034 zZi = equilDKE_Fourier.zZi;
0035 zmi = equilDKE_Fourier.zmi;
0036 pZeff = equilDKE_Fourier.xZeff;
0037 Rp = equilDKE_Fourier.Rp;
0038 Zp = equilDKE_Fourier.Zp;
0039 %
0040 ptr = sqrt(ptx.^2 + pty.^2); %minor radius on the (psi,stheta) grid
0041 ptBP = sqrt(ptBx.^2 + ptBy.^2); %poloidal magnetic field on the (psi,stheta) grid
0042 ptB = sqrt(ptBP.^2 + ptBphi.^2); %magnetic field on the (psi,stheta) grid
0043 %
0044 if ~isinf(Rp), %toroidal FP calculations
0045     if ~isfield(ray,'sphi') | ~isfield(ray,'sn'),% check for toroidal ray tracing
0046         error('cylindrical ray tracing and toroidal FP codes are incompatible');
0047     else
0048         sphi = ray.sphi;% Ray location (toroidal angle) along the propagation [1,ns]
0049         sn = ray.sn;% Ray wave number (toroidal) along the propagation [1,ns]
0050     end
0051     tor_mode = 1;
0052 else
0053     if ~isfield(ray,'sz') | ~isfield(ray,'skz'),% check for cylindrical ray tracing
0054         error('toroidal ray tracing and cylindrical FP codes are incompatible');
0055     else
0056         sz = ray.sz;% Ray longitudinal position
0057         skz = ray.skz;% Ray wave vector along the axis [1,ns]
0058     end
0059     tor_mode = 0;
0060 end
0061 %
0062 sx = ray.sx;% Ray location (major radius) along the propagation [1,ns]
0063 sy = ray.sy;% Ray location (vertical) along the propagation [1,ns]
0064 ss = ray.ss;% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0065 spsin = ray.spsin;% Ray location (normalized magnetic flux) along the propagation [1,ns]
0066 stheta = ray.stheta;% Ray location (poloidal angle) along the propagation [1,ns]
0067 %
0068 skx = ray.skx;% Ray wave vector (major radius) along the propagation [1,ns]
0069 sky = ray.sky;% Ray wave vector (vertical) along the propagation [1,ns]
0070 %
0071 if isfield(ray,'skrho') && isfield(ray,'sm')
0072     skrho = ray.skrho;% Ray wave vector (radial) along the propagation [1,ns]
0073     sm = ray.sm;% Ray wave number (poloidal) along the propagation [1,ns]
0074 else
0075     skrho = [];
0076     sm = [];
0077 end
0078 %
0079 %
0080 sNpar = ray.sNpar;% Ray parallel wave number [1,ns]
0081 sNperp = ray.sNperp;% Ray perpendicular wave number [1,ns]
0082 sdNpar = ray.sdNpar;% Ray spectral width in Npar [1,ns]
0083 %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar)
0084 %
0085 if isfield(ray,'swidth'),
0086     swidth = ray.swidth;
0087 else
0088     swidth = [];
0089 end
0090 %
0091 sepol_pmz = ray.sepol_pmz;% polarization vector [3,ns] complex
0092 sphi_xyz = ray.sphi_xyz;% normalized power flow [3,ns]
0093 %
0094 % Calculates ss (distance along the ray) if missing.
0095 %
0096 if isempty(ss),
0097     if tor_mode == 1,
0098         if isempty(sx) | isempty(sphi) | isempty(sy)
0099             error('Ray-tracing data incomplete')
0100         end
0101         ns = length(sx);
0102         ss(1) = 0;
0103         if ns > 1,
0104             sR = Rp + sx;
0105             sRh = (sR(1:ns-1) + sR(2:ns))/2;
0106             ds = sqrt(diff(sx).^2 + (sRh.*diff(sphi)).^2 + diff(sy).^2);
0107             ss(2:ns) = cumsum(ds);
0108         end
0109     else
0110         if isempty(sx) | isempty(sy) | isempty(sz)
0111             error('Ray-tracing data incomplete')
0112         end
0113         ns = length(sx);
0114         ss(1) = 0;
0115         if ns > 1,
0116             ds = sqrt(diff(sx).^2 + diff(sy).^2 + diff(sz).^2);
0117             ss(2:ns) = cumsum(ds);
0118         end
0119     end
0120 else 
0121     ns = length(ss);
0122 end
0123 %
0124 %mask = [1,find(diff(ss) > 0) + 1];
0125 %
0126 % Calculates spsin and stheta (normalized flux and poloidal angle, location along the ray) and sR sZ sphi if missing.
0127 %
0128 if isempty(spsin) | isempty(stheta),
0129     if isempty(sx) | isempty(sy)
0130         error('Ray-tracing data incomplete')
0131     end
0132     %
0133     % ray position in (r,theta) coordinates
0134     %
0135     sx(sx == 0) = eps; %to avoid singularities
0136     %sy(sy == 0) = eps; %to avoid singularities
0137     sr = sqrt(sx.^2 + sy.^2); %minor radius
0138     ssy = sign(sy);
0139     ssy(ssy == 0) = 1;
0140     stheta = atan(sy./sx) + pi - pi*(sx > 0).*ssy; %poloidal angle [0,2*pi]
0141     %
0142     if ns > 1,
0143         %
0144         % interpolation for r at ray poloidal locations
0145         %
0146         pstr = interp1(theta,ptr',stheta,method)';%r(psi,st)
0147         %
0148         % interpolation for psi at ray position (along cst theta lines)
0149         %
0150         %spsin = interp1(pstr,ppsin,sr,method);%psi(s)
0151         method2 = 'linear';
0152         for is = 1:ns 
0153             spsin(is) = interp1(pstr(:,is),ppsin,sr(is),method2,1);%psi(s)
0154         end
0155         %
0156         % remove ray parts outlise the FS.
0157         %
0158         %imask_1 = min(find(~isnan(spsin)));
0159         %imask_2 = min([find(isnan(spsin(imask_1:ns))),ns - imask_1 + 2]);
0160         %mask = [imask_1:imask_1 + imask_2 - 2];
0161         %
0162         % positive values only for spsin
0163         %
0164     else
0165         spsin = ppsin;
0166     end
0167     %
0168     if min(spsin) < -0
0169         error('The psin interpolation has produced negative values')
0170     end
0171     spsin(spsin <= 0) = eps;
0172 else
0173     if isempty(sx) | isempty(sy),
0174         sx = interp2(theta,ppsin,ptx,stheta,sypsin,method);
0175         sy = interp2(theta,ppsin,pty,stheta,sypsin,method);
0176         if tor_mode == 1,
0177             if ~isempty(sNpar)
0178                 ssigmah = sign(real(sNpar(1:ns-1,:) + sNpar(2:ns,:)));
0179             elseif ~isempty(sn)
0180                 ssigmah = sign(sn(1:ns-1,:) + sn(2:ns,:));
0181             else
0182                 error('propagation data is imcomplete')
0183             end
0184             sphi(1,:) = 0;
0185             sR = Rp + sx;
0186             sRh = (sR(1:ns-1,:) + sR(2:ns,:))/2;
0187             dsphi = sqrt(diff(ss).^2 - diff(sx).^2 - diff(sy).^2).*ssigmah./sRh;
0188             sphi(2:ns) = cumsum(dsphi);
0189         else
0190             if ~isempty(sNpar)
0191                 ssigmah = sign(real(sNpar(1:ns-1,:) + sNpar(2:ns,:)));
0192             elseif ~isempty(skz)
0193                 ssigmah = sign(skz(1:ns-1,:) + skz(2:ns,:));
0194             else
0195                 error('propagation data is imcomplete')
0196             end
0197             sz = sqrt(diff(ss).^2 - diff(sx).^2 - diff(sy).^2).*ssigmah;
0198         end            
0199     end
0200     sr = sqrt(sx.^2 + sy.^2); %minor radius
0201     %
0202 %    is0 = min(find(spsin < 1)); %position where the ray enters the plasma
0203 %    ismax = min([find(spsin >= 1 & ss > ss(is0)) - 1,ns]); %position where the ray leaves the plasma (if applicable)
0204 %    mask = is0:ismax;
0205 end    
0206 %
0207 if isfield(ray,'salphaphi_lin'),
0208     salphaphi_lin = ray.salphaphi_lin;
0209 else
0210     salphaphi_lin = NaN*ones(1,ns);
0211     %opt_disp = opt_disp + i;
0212 end
0213 %
0214 if isfield(ray,'salphaphi_coll') && colldamp == 1,
0215     salphaphi_coll = ray.salphaphi_coll;
0216 else
0217     salphaphi_coll = NaN*ones(1,ns);
0218 end
0219 %
0220 if isfield(ray,'snhui'),
0221     snhui = ray.snhui;
0222 else
0223     snhui = NaN*ones(1,ns);
0224 end
0225 % %
0226 % % remove unwanted ray parts
0227 % %
0228 % ns = length(mask);
0229 % %
0230 % spsin = spsin(mask);%part of ray outside the last flux-surface
0231 % stheta = stheta(mask);
0232 % %
0233 % sr = sr(mask);
0234 % %
0235 % sx = sx(mask);
0236 % sy = sy(mask);
0237 % %
0238 % if tor_mode == 1,
0239 %     sphi = sphi(mask);
0240 %     if ~isempty(sn),
0241 %         sn = sn(mask);
0242 %     end
0243 % else
0244 %     sz = sz(mask);
0245 %     if ~isempty(skz),
0246 %         skz = skz(mask);
0247 %     end
0248 % end
0249 % %
0250 % ss = ss(mask);
0251 % ss = ss - ss(1);%reinitialize ss
0252 % %
0253 % if ~isempty(skx) & ~isempty(skZ)
0254 %     skx = skx(mask);
0255 %     sky = sky(mask);
0256 % end
0257 % if ~isempty(skrho) & ~isempty(sm)
0258 %     skrho = skrho(mask);
0259 %     sm = sm(mask);
0260 % end
0261 % if ~isempty(sNpar)
0262 %     sNpar = sNpar(mask);
0263 % end
0264 % if ~isempty(sNperp)
0265 %     sNperp = sNperp(mask);
0266 % end
0267 % if ~isempty(sdNpar)
0268 %     sdNpar = sdNpar(mask);
0269 % end
0270 % if ~isempty(swidth)
0271 %     swidth = swidth(mask);
0272 % end
0273 % if ~isempty(sepol_pmz)
0274 %     sepol_pmz = sepol_pmz(:,mask);
0275 % end
0276 % if ~isempty(sphi_xyz)
0277 %     sphi_xyz = sphi_xyz(:,mask);
0278 % end
0279 % if ~isempty(salphaphi_lin)
0280 %     salphaphi_lin = salphaphi_lin(:,mask);
0281 % end
0282 %
0283 % Calculates Equilibrium properties along the ray
0284 %
0285 if ns > 1,
0286     %
0287     if isfield(ray,'sBx') & isfield(ray,'sBy') & isfield(ray,'sBphi'),
0288         sBx = ray.sBx;
0289         sBy = ray.sBy;
0290         sBphi = ray.sBphi;            
0291     else
0292         sBx = interp2(theta,ppsin,ptBx,stheta,spsin,method);
0293         sBy = interp2(theta,ppsin,ptBy,stheta,spsin,method);
0294         sBphi = interp2(theta,ppsin,ptBphi,stheta,spsin,method);
0295     end
0296     %
0297     if isfield(ray,'sB')
0298         sB = ray.sB;
0299     else
0300         sB = interp2(theta,ppsin,ptB,stheta,spsin,method);
0301     end
0302     %
0303     if isfield(ray,'sne') & isfield(ray,'szni') & isfield(ray,'sTe') & isfield(ray,'szTi'),
0304         sne = ray.sne;
0305         szni = ray.szni;
0306         sTe = ray.sTe;
0307         szTi = ray.szTi;
0308     else
0309         sne = interp1(ppsin,pne,spsin,method);
0310         szni = interp1(ppsin,pzni',spsin,method)';
0311         sTe = interp1(ppsin,pTe,spsin,method);
0312         szTi = interp1(ppsin,pzTi',spsin,method)';
0313     end
0314     %
0315     if isfield(ray,'srho'),
0316         srho = ray.srho;
0317     else
0318         srho = interp1(ppsin,prho,spsin,method);
0319     end
0320 else
0321     sBx = ptBx;
0322     sBphi = ptBphi;
0323     sBy = ptBy;
0324     sB = ptB;
0325     sne = pne;
0326     szni = pzni;
0327     sTe = pTe;
0328     szTi = pzTi;
0329     srho = prho;
0330 end    
0331 %
0332 % Calculates Npar and Nperp from kx,n and kZ if needed.
0333 %
0334 if isempty(sNpar),
0335     if tor_mode == 1,
0336         if isempty(skx) | isempty(sn) | isempty(sky) 
0337             error('Ray-tracing data incomplete')
0338         end
0339     %    if isempty(sR)
0340     %        sR = interp2(theta,ppsin,ptR,st,spsin,method);
0341     %    end
0342         %
0343         sR = Rp + sx;
0344         skpar = (skx.*sBx + sn.*sBphi./sR + sky.*sBy)./sqrt(sBx.^2 + sBphi.^2 + sBy.^2);
0345         sNpar = skpar*clum/omega_rf;
0346         %
0347         if isempty(sNperp), 
0348             skperp = sqrt(skx.^2 + (sn./sR).^2 + sky.^2 - skpar.^2);
0349             sNperp = skperp*clum/omega_rf;
0350         end
0351     else
0352         if isempty(skx) | isempty(sky) | isempty(skz) 
0353             error('Ray-tracing data incomplete')
0354         end
0355         %
0356         skpar = (skx.*sBx + sky.*sBy + skz.*sBz)./sqrt(sBx.^2 + sBy.^2 + sBz.^2);
0357         sNpar = skpar*clum/omega_rf;
0358         %
0359         if isempty(sNperp), 
0360             skperp = sqrt(skx.^2 + sky.^2 + skz.^2 - skpar.^2);
0361             sNperp = skperp*clum/omega_rf;
0362         end        
0363     end
0364 end
0365 %
0366 if real(opt_disp) == 0, %sepol_pmz, sNperp and phi_xyz given by ray-tracing
0367     if isempty(sepol_pmz) | isempty(sNperp) | isempty(sphi_xyz) 
0368         error('Ray-tracing data incomplete')
0369     end
0370 else
0371 %
0372 % Solves the dispersion relation for sepol_pmz, sNperp and phi_xyz
0373 %
0374     if real(opt_disp) == 1, % cold plasma dispersion tensor for sepol_pmz, sNperp and phi_xyz (colddisp_dke_jd)
0375         for is = 1:ns
0376             [sNperpp(is),sNperpm(is),sKxyz(1:3,1:3,is),...
0377                 sep_xyz(1:3,is),sem_xyz(1:3,is),sep_pmz(1:3,is),sem_pmz(1:3,is),sep_pyk(1:3,is),sem_pyk(1:3,is),...
0378                 sphip_xyz(1:3,is),sphim_xyz(1:3,is),sphip_pmz(1:3,is),sphim_pmz(1:3,is),sphip_pyk(1:3,is),sphim_pyk(1:3,is)]...
0379             = colddisp_dke_jd(real(sNpar(is)),omega_rf,sne(is),szni(:,is),zZi,zmi,sB(is));
0380         end
0381         if length(sNperp) == ns,
0382             smaskp = (abs(sNperpp - sNperp) <= abs(sNperpm - sNperp));
0383             sNperp = sNperpp.*smaskp + sNperpm.*~smaskp;
0384             sepol_pmz = sep_pmz.*(ones(3,1)*smaskp) + sem_pmz.*(ones(3,1)*~smaskp);
0385             sphi_xyz = sphip_xyz.*(ones(3,1)*smaskp) + sphim_xyz.*(ones(3,1)*~smaskp);
0386         elseif length(sNperp) == 1
0387             %mask0 = (sNperp == 0)%O-mode;
0388             %%%%%%%%%% to be implemented %%%%%%%%%%
0389             error('model not yet implemented');
0390         else
0391             error('No information was given to choose between the two cold modes')
0392         end
0393     elseif real(opt_disp) == 2, % non-relativistic hot plasma dispersion tensor for sepol_pmz, sNperp and phi_xyz (hotdisp_dke_jd)
0394         for is=1:ns,
0395             disp(['Hot plasma dispersion relation: ',num2str(is),'/',num2str(ns)])
0396             [sNperp(is),e_xyz,sepol_pmz(:,is),e_pyk,D,X,X_s,X_sn,salphaphi_lin(is),sphiT_xyz(:,is),sphiP_xyz(:,is)]...
0397                 = hotdisp_dke_jd(real(sNpar(is)),omega_rf,sTe(is),sne(is),szTi(:,is),szni(:,is),zZi,zmi,sB(is),sNperp(is),1);
0398         end
0399         %
0400         smask = ~isnan(sNperp);
0401         ns = sum(smask);
0402         %
0403         sNperp = sNperp(smask);
0404         sepol_pmz = sepol_pmz(:,smask);
0405         salphaphi_lin = salphaphi_lin(smask);
0406         sphiT_xyz = sphiT_xyz(:,smask);
0407         sphiP_xyz = sphiP_xyz(:,smask);
0408         %
0409         sphi_xyz = real(sphiP_xyz + sphiT_xyz);
0410         salphaphi_lin = real(salphaphi_lin)*omega_rf/clum;
0411         %
0412         ss = ss(smask);
0413         spsin = spsin(smask);
0414         stheta = stheta(smask);
0415         %
0416         sr = sr(smask);
0417         sx = sx(smask);
0418         sy = sy(smask);
0419         %
0420         if tor_mode == 1,
0421             sphi = sphi(smask);
0422             if ~isempty(sn)
0423                 sn = sn(smask);
0424             end               
0425         else
0426             sz = sz(smask);
0427             if ~isempty(skz)
0428                 skz = skz(smask);
0429             end               
0430         end
0431         %
0432         sNpar = sNpar(smask);
0433         sdNpar = sdNpar(smask);
0434         %
0435         if ~isempty(skx) && ~isempty(sky)
0436             skx = skx(smask);
0437             sky = sky(smask);
0438         end        
0439         %
0440         if ~isempty(skrho) && ~isempty(sm)
0441             skrho = skrho(smask);
0442             sm = sm(smask);
0443         end        
0444         %
0445         if ~isempty(swidth)
0446             swidth = swidth(smask);
0447         end        
0448         %
0449         sB = sB(smask);
0450         sBx = sBx(smask);
0451         sBy = sBy(smask);
0452         sBphi = sBphi(smask);
0453         %
0454         sne = sne(smask);
0455         szni = szni(:,smask);
0456         sTe = sTe(smask);
0457         szTi = szTi(:,smask);
0458         %
0459         srho = srho(smask);
0460         %
0461     elseif real(opt_disp) == 3, % non-relativistic hot plasma dispersion tensor for sepol_pmz, sNperp and phi_xyz (R2D2)
0462         [smask_ht,sNperp_ht,sepol_pmz_ht,sPflow_ht,sPabs_at] = ...
0463             disp_hot_jd(ss,sB,sne,sTe,real(sNpar),omega_rf,sNperp,0); 
0464 %
0465         ns = length(smask_ht);
0466         %
0467         sNperp = sNperp_ht(smask_ht);
0468         sepol_pmz = sepol_pmz_ht(:,smask_ht);
0469 %         sphi_xyz(1,:) = abs(2*sPflow_ht(1,smask_ht)/e0/clum);
0470 %         sphi_xyz(2,:) = zeros(1,ns);
0471 %         sphi_xyz(3,:) = zeros(1,ns);
0472         sphi_xyz = abs(2*sPflow_ht(:,smask_ht)/e0/clum);
0473         salphaphi_lin = 2*sPabs_at(smask_ht)/e0/clum;
0474         %
0475         ss = ss(smask_ht);
0476         spsin = spsin(smask_ht);
0477         stheta = stheta(smask_ht);
0478         %
0479         sr = sr(smask_ht);
0480         sx = sx(smask_ht);
0481         sy = sy(smask_ht);
0482         %
0483         if tor_mode == 1,
0484             sphi = sphi(smask_ht);
0485             if ~isempty(sn)
0486                 sn = sn(smask_ht);
0487             end               
0488         else
0489             sz = sz(smask_ht);
0490             if ~isempty(skz)
0491                 skz = skz(smask_ht);
0492             end               
0493         end
0494         %
0495         sNpar = sNpar(smask_ht);
0496         sdNpar = sdNpar(smask_ht);
0497         %
0498         if ~isempty(skx) && ~isempty(sky)
0499             skx = skx(smask_ht);
0500             sky = sky(smask_ht);
0501         end        
0502         %
0503         if ~isempty(skrho) && ~isempty(sm)
0504             skrho = skrho(smask_ht);
0505             sm = sm(smask_ht);
0506         end        
0507         %
0508         if ~isempty(swidth)
0509             swidth = swidth(smask_ht);
0510         end        
0511         %
0512         sB = sB(smask_ht);
0513         sBx = sBx(smask_ht);
0514         sBy = sBy(smask_ht);
0515         sBphi = sBphi(smask_ht);
0516         %
0517         sne = sne(smask_ht);
0518         szni = szni(:,smask_ht);
0519         sTe = sTe(smask_ht);
0520         szTi = szTi(:,smask_ht);
0521         %
0522         srho = srho(smask_ht);
0523         %
0524     elseif real(opt_disp) == 4, % relativistic hot plasma dispersion tensor for sepol_pmz, sNperp and phi_xyz (R2D2)
0525         %%%%%%%%%% to be implemented %%%%%%%%%%
0526         error('model not yet implemented');
0527     elseif real(opt_disp) == 5, % electrostatic approximation for EBW calculations
0528         %
0529         [sNperp,sepol_pmz,sphi_xyz,salphaphi_lin_nr,salphaphi_lin_wr,salphaphi_lin_fr] = ES_ebw_disp_jd(real(sNpar),omega_rf,sTe,sne,sB,sNperp);%,0,imag(opt_disp)
0530         %
0531         salphaphi_lin = real(salphaphi_lin_fr)*omega_rf/clum;
0532         opt_disp = real(opt_disp);
0533         %
0534         smask = ~isnan(sNperp);
0535         ns = sum(smask);
0536         %
0537         sNperp = sNperp(smask);
0538         sepol_pmz = sepol_pmz(:,smask);
0539         salphaphi_lin = salphaphi_lin(smask);
0540         sphi_xyz = sphi_xyz(:,smask);
0541         %
0542         ss = ss(smask);
0543         spsin = spsin(smask);
0544         stheta = stheta(smask);
0545         %
0546         sr = sr(smask);
0547         sx = sx(smask);
0548         sy = sy(smask);
0549         %
0550         if tor_mode == 1,
0551             sphi = sphi(smask);
0552             if ~isempty(sn)
0553                 sn = sn(smask);
0554             end               
0555         else
0556             sz = sz(smask);
0557             if ~isempty(skz)
0558                 skz = skz(smask);
0559             end               
0560         end
0561         %
0562         sNpar = sNpar(smask);
0563         sdNpar = sdNpar(smask);
0564         %
0565         if ~isempty(skx) && ~isempty(sky)
0566             skx = skx(smask);
0567             sky = sky(smask);
0568         end        
0569         %
0570         if ~isempty(skrho) && ~isempty(sm)
0571             skrho = skrho(smask);
0572             sm = sm(smask);
0573         end        
0574         %
0575         if ~isempty(swidth)
0576             swidth = swidth(smask);
0577         end        
0578         %
0579         sB = sB(smask);
0580         sBx = sBx(smask);
0581         sBy = sBy(smask);
0582         sBphi = sBphi(smask);
0583         %
0584         sne = sne(smask);
0585         szni = szni(:,smask);
0586         sTe = sTe(smask);
0587         szTi = szTi(:,smask);
0588         %
0589         srho = srho(smask);
0590         %
0591     else
0592         error('invalid disp option');
0593     end
0594 end
0595 %
0596 % Calculates sdNpar if needed.
0597 %
0598 if isempty(sdNpar) || any(isnan(sdNpar)),
0599     if ~isempty(swidth),
0600         sphi_norm = sqrt(abs(sphi_xyz(1,:).^2 + sphi_xyz(2,:).^2 + sphi_xyz(3,:).^2));
0601         sdNpar = clum.*abs(sphi_xyz(1,:))./sphi_norm/omega_rf./swidth;
0602     else
0603         sdNpar = 0.01*ones(1,ns);%assume a gaussian spectrum with this arbitrary width
0604     end
0605 end
0606 %
0607 if abs(imag(opt_disp)) >= 1,%calculation of the linear absorption coefficient
0608     %
0609     rel_opt = imag(opt_disp) > 0;
0610     %
0611     sa = sne*qe^2/e0/me/omega_rf^2;
0612     sb = sB*qe/me/omega_rf;
0613     sbeta = sqrt(sTe/mc2);
0614     %
0615     salphaphi_lin = zeros(1,ns);
0616     %
0617     for n = n_rf_list,
0618         if abs(imag(opt_disp)) == 1,
0619             salphaphi_lin = salphaphi_lin + omega_rf/clum*alphaphi_fr_jd(sa,sb,sbeta,real(sNpar),real(sNperp),sepol_pmz,n,rel_opt);
0620         else
0621             salphaphi_lin = salphaphi_lin + omega_rf/clum*alphaphi_fr_jd(sa,sb,sbeta,real(sNpar) + sign(real(sNpar)).*ss*abs(imag(opt_disp) - 1)/2,real(sNperp),sepol_pmz,n,rel_opt);
0622         end
0623     end
0624     %
0625    
0626 end
0627 %
0628 if colldamp == 1,%calculation of the collisional damping absorption coefficient
0629     %
0630     if ~isfield(ray,'sTe')
0631         if exist('sTe'),
0632             ray.sTe = sTe;
0633         else 
0634             colldamp = 0;
0635             disp('WARNING: Collisional damping is turned-off because ray.sTe missing.')
0636         end
0637     end
0638     %
0639     if ~isfield(ray,'sne')
0640         if exist('sne'),
0641             ray.sne = sne;
0642         else
0643             colldamp = 0;
0644             disp('WARNING: Collisional damping is turned-off because ray.sne missing.')
0645         end
0646     end
0647     %
0648     if ~isfield(ray,'szni')
0649         if exist('szni'),
0650             ray.szni = szni;
0651         else
0652             colldamp = 0;
0653             disp('WARNING: Collisional damping is turned-off because ray.szni missing.')
0654         end
0655     end
0656     %
0657     if ~isfield(ray,'swpe2') 
0658         colldamp = 0;
0659         disp('WARNING: Collisional damping is turned-off because ray.swpe2 missing.')
0660     end
0661     %
0662     if ~isfield(ray,'swce2') 
0663         colldamp = 0;
0664         disp('WARNING: Collisional damping is turned-off because ray.swce2 missing.')
0665     end
0666     %
0667     salphaphi_coll = zeros(1,ns);
0668     %
0669     [xalphaphi_coll,snhui] = alphaphi_cd_yp(ray,equilDKE_Fourier,omega_rf);
0670     salphaphi_coll = salphaphi_coll + omega_rf/clum*xalphaphi_coll;%calculation of the non-resonant collisional absorption coefficient
0671 end
0672 %
0673 ray.sx = sx;% Ray location (major radius) along the propagation [1,ns]
0674 ray.sy = sy;% Ray location (vertical) along the propagation [1,ns]
0675 ray.ss = ss;% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0676 ray.spsin = spsin;% Ray location (normalized magnetic flux) along the propagation [1,ns]
0677 ray.stheta = stheta;% Ray location (poloidal angle) along the propagation [1,ns]
0678 %
0679 if tor_mode == 1,
0680     ray.sphi = sphi;% Ray location (toroidal angle) along the propagation [1,ns]
0681     if ~isempty(sn)
0682         ray.sn = sn;% Ray wave number (toroidal) along the propagation [1,ns]
0683     end               
0684 else
0685     ray.sz = sz;% Ray location (toroidal angle) along the propagation [1,ns]
0686     if ~isempty(skz)
0687         ray.skz = skz;
0688     end               
0689 end
0690 %
0691 ray.skx = skx;% Ray wave vector (major radius) along the propagation [1,ns]
0692 ray.sky = sky;% Ray wave vector (vertical) along the propagation [1,ns]
0693 %
0694 ray.skrho = skrho;% Ray wave vector (radial) along the propagation [1,ns]
0695 ray.sm = sm;% Ray wave number (poloidal) along the propagation [1,ns]
0696 %
0697 ray.sNpar = sNpar;% Ray parallel wave number [1,ns]
0698 ray.sNperp = real(sNperp);% Ray perpendicular wave number [1,ns]
0699 ray.sdNpar = sdNpar;% Ray spectral width in Npar [1,ns]
0700 %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar)
0701 %
0702 ray.sepol_pmz = sepol_pmz;% polarization vector [3,ns] complex
0703 ray.sphi_xyz = sphi_xyz;% normalized power flow [3,ns]
0704 %
0705 % Additional ray parameters
0706 %
0707 ray.sr = sr;
0708 %
0709 ray.sB = sB;
0710 ray.sBx = sBx;
0711 ray.sBy = sBy;
0712 ray.sBphi = sBphi;
0713 %
0714 ray.sne = sne;
0715 ray.szni = szni;
0716 ray.sTe = sTe;
0717 ray.szTi = szTi;
0718 %
0719 ray.srho = srho;
0720 %
0721 ray.salphaphi_lin = salphaphi_lin;
0722 %
0723 if colldamp == 1,
0724     ray.salphaphi_coll = salphaphi_coll;
0725     ray.snhui = snhui;
0726 end
0727 %
0728 ray.swidth = swidth;
0729 %
0730 mask = [1,1 + find(diff(ray.ss) > 0)];
0731 %
0732 fnames = fieldnames(ray);
0733 %
0734 ns = length(ray.ss);
0735 %
0736 for iname = 1:length(fnames),
0737     %
0738     data = ray.(fnames{iname});
0739     %
0740     sdata = size(data);
0741     %
0742     for is = 1:length(sdata),
0743         %
0744         if sdata(is) == ns,
0745             data = shiftdim(data,is - 1);
0746             data = data(mask,:);
0747             data = shiftdim(data,length(sdata) - is + 1);
0748         end
0749         %
0750         ray.(fnames{iname}) = data;
0751         %
0752     end
0753     %
0754 end
rayprocess_jd

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
