propagation_jd

 This function couples the DKE code with a ray-tracing code.
 RF rays are distributed onto the DKE flux surfaces.
 Rays that "turn back" are divided into several subrays. 
 
   INPUTS: 

       - rays: list of RF rays. {1,nb_rf}. rays{ib} = ray is a structure 
                   that includes all characteristics along the ray path  
       - xpsin_S: flux radial grid [1,nr+1]
       - xpsin: distribution radial grid [1,nr]
       - method: interpolation method in psi (default = spline)
       - display_mode: (1) display calculation result

   OUTPUTS:

       - waveparam: strucure for all ray fragments
       - xyprop: RF rays propagation factor (1 if increasing psin, -1 if
                   decreasing psin, 0 if flux-surface (FS) is not visited by the ray)
       - xys: RF rays path length at center of FS
       - xys: RF rays linear power flow at center of FS
       - yb: list of ray starts in y grid and corresponding ray number. [2,nbb]

 by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)

0001 function [waveparam,xyprop,xys,xyPinc_2piRp_lin,xyPinc_2piRp_coll,xyPinc_2piRp_noabs,yb,yP0_2piRp] = propagation_jd(rays,xpsin_S,xpsin,method,display_mode)
0002 %
0003 % This function couples the DKE code with a ray-tracing code.
0004 % RF rays are distributed onto the DKE flux surfaces.
0005 % Rays that "turn back" are divided into several subrays.
0006 %
0007 %   INPUTS:
0008 %
0009 %       - rays: list of RF rays. {1,nb_rf}. rays{ib} = ray is a structure
0010 %                   that includes all characteristics along the ray path
0011 %       - xpsin_S: flux radial grid [1,nr+1]
0012 %       - xpsin: distribution radial grid [1,nr]
0013 %       - method: interpolation method in psi (default = spline)
0014 %       - display_mode: (1) display calculation result
0015 %
0016 %   OUTPUTS:
0017 %
0018 %       - waveparam: strucure for all ray fragments
0019 %       - xyprop: RF rays propagation factor (1 if increasing psin, -1 if
0020 %                   decreasing psin, 0 if flux-surface (FS) is not visited by the ray)
0021 %       - xys: RF rays path length at center of FS
0022 %       - xys: RF rays linear power flow at center of FS
0023 %       - yb: list of ray starts in y grid and corresponding ray number. [2,nbb]
0024 %
0025 % by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)
0026 %
0027 if nargin < 5
0028     display_mode = 0;
0029 end
0030 if nargin < 4
0031     method = 'linear';%interpolation technique in psi
0032 end
0033 %
0034 nb = length(rays);
0035 nr = length(xpsin_S) - 1;
0036 %
0037 xyprop = zeros(nr,0); %propagation factor (NaN is ray start, +1 is increasing psi, -1 is decreasing psi, 0 is non-applicable)
0038 xyPinc_2piRp_lin = zeros(nr,0);% power flow in the ray based on linear absorption calculation at center of flux-surface
0039 xyPinc_2piRp_coll = zeros(nr,0);% power flow in the ray based on non resonant collisional absorption calculation at center of flux-surface
0040 xyPinc_2piRp_noabs = zeros(nr,0);
0041 xys = zeros(nr,0);% ray path length at center of flux-surface
0042 yb = []; %list of ray starts
0043 yP0_2piRp = []; %list of initial powers
0044 %
0045 iy = 0; %counting global number of ray
0046 %
0047 yiy_rf = [];%Index of the ray number ray fragment belongs to
0048 yir_rf = [];%Index of the flux surface the ray fragment belongs to
0049 %ys_rf = [];%Ray path length coordinate at center of the ray fragment
0050 yds_rf = [];%incremental length of the ray fragment
0051 yomega_rf = [];%RF wave frequency
0052 yopt_rf = [];%Options for Larmor radius effects and/or EBW calc.
0053 ythetab_rf = [];%Poloidal location of RF beam (relevant for bounce_mode > 0 only)
0054 yB_rf =  [];%Magnetic field
0055 yNpar_rf = [];%Gaussian N// spectrum: peak n//
0056 ydNpar_rf = [];%Gaussian N// spectrum: width in n//
0057 yNperp_rf = [];%RF wave : n_perp
0058 yepolp_rf = [];%RF wave : polarization plus
0059 yepolm_rf = [];%RF wave : polarization minus
0060 yepolz_rf = [];%RF wave : polarization parallel
0061 yphinorm_rf = [];%RF wave : amplitude of the normalized energy flow density
0062 %
0063 % Each ray has to be treated separately, because their lengths differ.
0064 %
0065 time0 = clock;
0066 %
0067 for ib = 1:nb
0068     %
0069     ray = rays{ib};
0070     ss = ray.ss;
0071     spsin = ray.spsin;
0072     %
0073     dsmin = ray.dsmin;
0074     %
0075     ns = length(ss);
0076     %
0077     if ns > 1,
0078         sphinorm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0079     else
0080         sphinorm = abs(ray.sphi_xyz(1,:));% only the perpendicular contribution if to ray direction is given
0081     end
0082     %
0083     if isempty(find(spsin < 1)), 
0084         disp('Warning: one ray did not enter the plasma and was rejected')
0085         continue
0086     end
0087     %
0088     if min(spsin) < 0 | max(isnan(spsin)), 
0089         disp('Warning: one ray is ill defined (spsin < 0 or NaN) and was rejected')
0090         continue
0091     end
0092     %
0093     if max(spsin) > 1, 
0094         disp('Warning: ray positions outside the last FS are projected on that last FS')
0095         spsin(find(spsin > 1)) = 1;
0096     end
0097     %
0098     %is = max([1,min(find(spsin < 1)) - 1]);
0099     %
0100     yb = [yb,[iy+1;ib]];%keeps track of new ray indices.
0101     yP0_2piRp = [yP0_2piRp,ray.P0_2piRp];
0102     %
0103 %
0104 % Find ray subdivisions
0105 %
0106     sdpsin = diff(spsin);
0107     %
0108     svar = sign(sdpsin);
0109     srelvar = sign(sdpsin(1:ns-2).*sdpsin(2:ns-1));%relative psin variations between two consecutive steps
0110     %
0111     islim = find(srelvar == -1) + 1;% positions where ray turns around
0112 %
0113 % Treatment of stretches with constant psin
0114 %
0115     islim0 = find(srelvar == 0) + 1;% positions where ray stay at constant psin
0116     %
0117     start = 1;
0118     finish = 0;
0119     %
0120     for iislim0 = 1:length(islim0),
0121         %
0122         if start == 1,% beginning of stretch
0123             bdir_before = svar(islim0(iislim0) - 1);% 0 if ray start with constant psin stretch
0124             finish = 0;
0125             start = 0;
0126         end
0127         %
0128         if (iislim0 == length(islim0))...
0129          | (islim0(iislim0 + 1) > islim0(iislim0) + 1),%find end of stretch
0130             finish = 1;
0131         end
0132         %
0133         if finish == 1 & islim0(iislim0),% end of stretch
0134             bdir_after = svar(islim0(iislim0));% 0 if ray finishes with constant psin stretch
0135             if bdir_after*bdir_before == -1,
0136                 islim = [islim,islim0(iislim0)];% add this position to turning positions
0137             end
0138             start = 1;
0139         end
0140         %
0141     end
0142     %
0143     islim = sort([1,islim,ns]);
0144 %
0145 % Distribute subrays onto flux-surfaces
0146 %
0147     ny = length(islim) - 1;%number of subrays
0148     %
0149     for iiy = 1:ny,% for each subray
0150         %
0151         iy = iy + 1;
0152         %
0153         % stretch of ray that belongs to the subray
0154         %
0155         is0 = islim(iiy);
0156         is1 = islim(iiy+1);
0157         %
0158         ismask = is0:is1;
0159         ns_mask = length(ismask);
0160         %
0161         % find what flux-surface each fragment belongs to
0162         %
0163         xspsin_S = repmat(xpsin_S.',[1,ns_mask]);
0164         xspsin_ray = repmat(spsin(ismask),[nr,1]);
0165         %
0166         [sir,check] = find(xspsin_S(1:nr,:) < xspsin_ray & xspsin_S(2:nr+1,:) >= xspsin_ray);
0167         %
0168         sir = sir(:).';
0169         check = check(:).';
0170         %
0171         if length(check) ~= ns_mask | check ~= 1:ns_mask,
0172             error('The subray distribution failed');
0173         end
0174         %
0175         % find fragments that cross flux-surface limits
0176         %
0177         iscross = find(diff(sir) ~= 0);
0178         %
0179         % determines the limit of ray fragments
0180         %
0181         ir_cross = max([sir(iscross);sir(iscross + 1)]);%index of flux-surfaces being crossed
0182         %
0183         ssiscross = ss(iscross + is0 - 1) + (ss(iscross + is0) - ss(iscross + is0 - 1)).*...
0184             (xspsin_S(ir_cross) - spsin(iscross + is0 - 1))./...
0185             (spsin(iscross + is0) - spsin(iscross + is0 - 1));        
0186         %
0187 %         ssiscross = interp1([spsin(iscross + is0 - 1);spsin(iscross + is0)],...
0188 %                             [ss(iscross + is0 - 1);ss(iscross + is0)],...
0189 %                             xspsin_S(sir(iscross)+1),method);
0190         %
0191         % reduce number of fragments based on dsmin
0192         %
0193         if dsmin > 0,
0194             is_red = is0;
0195             ismask_red = [is0];
0196             while is_red < is1,
0197                 ss_red = ss(is_red);
0198                 is_red = min([is1,find(ss >= ss_red + dsmin)]);
0199                 ismask_red = [ismask_red,is_red];
0200             end
0201         else
0202             ismask_red = ismask;
0203         end
0204         %
0205         ss_S = [ss(ismask_red),ssiscross];%vector of fragment limits
0206         spsin_S = [spsin(ismask_red),xpsin_S(ir_cross)];%corresponding psin values
0207         %
0208         [ss_S,iss_S] = sort(ss_S);
0209         spsin_S = spsin_S(iss_S);
0210         %
0211         mask_S = [1,1+find(diff(ss_S) ~= 0)];
0212         ys_S = ss_S(mask_S);
0213         ypsin_S = spsin_S(mask_S);
0214         %
0215         if ns_mask > 1,
0216             ys = (ys_S(1:end-1) + ys_S(2:end))/2;
0217             ypsin = (ypsin_S(1:end-1) + ypsin_S(2:end))/2;
0218         else
0219             ys = ys_S;
0220             ypsin = ypsin_S;
0221         end            
0222         %
0223         npsin = length(ys);
0224         %
0225         xypsin_S = repmat(xpsin_S.',[1,npsin]);
0226         xypsin_ray = repmat(ypsin,[nr,1]);
0227         %
0228         [yir,check] = find(xypsin_S(1:nr,:) < xypsin_ray & xypsin_S(2:nr+1,:) >= xypsin_ray);
0229         %
0230         yir = yir(:).';
0231         check = check(:).';
0232         %
0233         if length(check) ~= npsin | check ~= 1:npsin,
0234             error('The subray fragmentation failed');
0235         end
0236         %
0237         yir_rf = [yir_rf,yir];
0238 %        ys_rf = [ys_rf,ys];
0239         %
0240         yds_rf = [yds_rf,diff(ys_S)];%incremental length of the ray path
0241         %
0242         yiy_rf = [yiy_rf,iy*ones(1,npsin)];%RF wave frequency
0243         yomega_rf = [yomega_rf,ray.omega_rf*ones(1,npsin)];%RF wave frequency
0244         yopt_rf = [yopt_rf,ray.opt_rf*ones(1,npsin)];%Options for Larmor radius effects and/or EBW calc.
0245         %
0246         % interpolation for the center of each fragment.
0247         %
0248         if ns_mask > 1,
0249             ythetab_rf = [ythetab_rf,interp1(ss(ismask),ray.stheta(ismask),ys,method)];%Poloidal location of RF beam (relevant for bounce_mode > 0 only)
0250             yB_rf = [yB_rf,interp1(ss(ismask),ray.sB(ismask),ys,method)];%Magnetic field
0251             yNpar_rf = [yNpar_rf,interp1(ss(ismask),ray.sNpar(ismask),ys,method)];%Gaussian N// spectrum: peak n//
0252             ydNpar_rf = [ydNpar_rf,interp1(ss(ismask),ray.sdNpar(ismask),ys,method)];%Gaussian N// spectrum: width in n//
0253             yNperp_rf = [yNperp_rf,interp1(ss(ismask),ray.sNperp(ismask),ys,method)];%RF wave : n_perp
0254             yepolp_rf = [yepolp_rf,interp1(ss(ismask),ray.sepol_pmz(1,ismask),ys,method)];%RF wave : polarization plus
0255             yepolm_rf = [yepolm_rf,interp1(ss(ismask),ray.sepol_pmz(2,ismask),ys,method)];%RF wave : polarization minus
0256             yepolz_rf = [yepolz_rf,interp1(ss(ismask),ray.sepol_pmz(3,ismask),ys,method)];%RF wave : polarization parallel
0257             yphinorm_rf = [yphinorm_rf,interp1(ss(ismask),sphinorm(ismask),ys,method)];%RF wave : amplitude of the normalized energy flow density
0258         else
0259             ythetab_rf = [ythetab_rf,ray.stheta(ismask)];%Poloidal location of RF beam (relevant for bounce_mode > 0 only)
0260             yB_rf = [yB_rf,ray.sB(ismask)];%Magnetic field
0261             yNpar_rf = [yNpar_rf,ray.sNpar(ismask)];%Gaussian N// spectrum: peak n//
0262             ydNpar_rf = [ydNpar_rf,ray.sdNpar(ismask)];%Gaussian N// spectrum: width in n//
0263             yNperp_rf = [yNperp_rf,ray.sNperp(ismask)];%RF wave : n_perp
0264             yepolp_rf = [yepolp_rf,ray.sepol_pmz(1,ismask)];%RF wave : polarization plus
0265             yepolm_rf = [yepolm_rf,ray.sepol_pmz(2,ismask)];%RF wave : polarization minus
0266             yepolz_rf = [yepolz_rf,ray.sepol_pmz(3,ismask)];%RF wave : polarization parallel
0267             yphinorm_rf = [yphinorm_rf,sphinorm(ismask)];%RF wave : amplitude of the normalized energy flow density
0268         end
0269         %
0270         % build xyprop and xyPinc_lin
0271         %
0272         ir0 = max(find(xpsin_S < spsin(is0)));%initial FS on ray trajectory
0273         ir1 = max(find(xpsin_S < spsin(is1)));%final FS on ray trajectory
0274         %
0275         xmask = min([ir0,ir1]):max([ir0,ir1]);
0276         %
0277         if spsin(is0) < spsin(is1),%ray propagating outward
0278             bdir = 1;
0279         else%ray propagating inward
0280             bdir = -1;
0281         end
0282         %
0283         xyprop(xmask,iy) = bdir;
0284         %
0285         mask_interp = [1,1 + find(diff(spsin(ismask)) ~= 0)];
0286         if ns_mask > 1,
0287             xys(xmask,iy) = interp1(spsin(ismask(mask_interp)),ss(ismask(mask_interp)),xpsin(xmask),method).';
0288         else
0289             xys(xmask,iy) = ray.ss(ismask);
0290         end
0291         %
0292         if isnan(xys(ir0,iy)),
0293             if ~isempty(ssiscross),
0294                 xys(ir0,iy) = (ss(is0) + ssiscross(1))/2;
0295             else
0296                 xys(ir0,iy) = (ss(is0) + ss(is1))/2;
0297             end
0298             %xyPinc_2piRp_lin(ir0,iy) = ray.sP_2piRp_lin(is0);
0299         end
0300         if isnan(xys(ir1,iy)),
0301             xys(ir1,iy) = (ssiscross(end) + ss(is1))/2;
0302             %xyPinc_2piRp_lin(ir1,iy) = ray.sP_2piRp_lin(is1);
0303         end
0304         %
0305         if isfield(ray,'sP_2piRp_lin')
0306               xyPinc_2piRp_lin(xmask,iy) = interp1(ss(ismask),ray.sP_2piRp_lin(ismask),xys(xmask,iy),method);
0307             %
0308             if isfield(ray,'sP_2piRp_coll'),
0309                 xyPinc_2piRp_coll(xmask,iy) = interp1(ss(ismask),ray.sP_2piRp_coll(ismask),xys(xmask,iy),method);
0310             else
0311                 xyPinc_2piRp_coll(xmask,iy) = ray.P0_2piRp*ones(size(xyPinc_2piRp_lin(xmask,iy)));%no collisional damping. P along ray is P0
0312             end
0313             %
0314         else
0315             xyPinc_2piRp_lin(xmask,iy) = ray.P0_2piRp;%no linear damping. P along ray is P0
0316             xyPinc_2piRp_coll(xmask,iy) = ray.P0_2piRp;%no collisional damping. P along ray is P0
0317         end
0318         %
0319         xyPinc_2piRp_noabs(xmask,iy) = ray.P0_2piRp;%P along ray is P0
0320     end
0321     if display_mode >= 1,info_dke_yp(2,['Preprocessing of launched ray #',int2str(ib),'/',int2str(nb),' done. Elapsed time : ',num2str(etime(clock,time0)),' (s)']);time0 = clock;end
0322 end
0323 %
0324 waveparam.model = 'RT';
0325 waveparam.biy_rf = yiy_rf;%Index of the ray the ray fragment belongs to
0326 waveparam.bir_rf = yir_rf;%Index of the flux surface the ray fragment belongs to
0327 %waveparam.bs_rf = ys_rf;%Ray path length coordinate at center of the ray fragment
0328 waveparam.bds_rf = yds_rf;%incremental length of the ray fragment
0329 waveparam.bomega_rf = yomega_rf;%RF wave frequency
0330 waveparam.bopt_rf = yopt_rf;%Options for Larmor radius effects and/or EBW calc.
0331 waveparam.bthetab_rf = ythetab_rf;%Poloidal location of RF beam (relevant for bounce_mode > 0 only)
0332 waveparam.bB_rf = yB_rf;%Magnetic field
0333 waveparam.bNpar_rf = yNpar_rf;%Gaussian N// spectrum: peak n//
0334 waveparam.bdNpar_rf = ydNpar_rf;%Gaussian N// spectrum: width in n//
0335 waveparam.bNperp_rf = yNperp_rf;%RF wave : n_perp
0336 waveparam.bepolp_rf = yepolp_rf;%RF wave : polarization plus
0337 waveparam.bepolm_rf = yepolm_rf;%RF wave : polarization minus
0338 waveparam.bepolz_rf = yepolz_rf;%RF wave : polarization parallel
0339 waveparam.bphinorm_rf = yphinorm_rf;%RF wave : amplitude of the normalized energy flow density
0340 %