rfdiff_dke_jd

 Calculates the RF quasilinear diffusion coefficient and the diffusion
 tensor elements on a 3-D Grid: 2-D (p,mhu) in momentum; 1-D (rho) in space. 
 The coefficients associated with each ray are saved in the file
 temp_DQL.mat while the resulting coefficients are returned in the
 variables ZYXXD_rf,ZYXXD_rf_tp and ZYXXF_rf_tp.
 Note that since the coefficients are saved, only the missing elements are
 calculated in rfdiff_dke_jd while others are retrived. The coefficients
 are finally weighted by the magnitude of the diffusion coefficient.

 INPUTS:

       - waveparam:        RF wave parameters (See detailed description below in code)
                           waveparam is build by the function rfwave_dke_jd.m
       - nmax_rf:          maximum harmonic number to be included
       - save_id:          identification string, based on running date and time
       - temppath:         path to temporary folder
       - theta...xqtilde:  equilibrium properties
       - Xmhup...Xdmhupp:  pitch angle grid 
       - Xpnpp...Xdpnpp:   momentum grid 
       - Xgamma...:        relativistic factor
       - Xxlambda...:      normalized bounce time
       - betath_ref:       sqrt(Te_ref/mc2)
       - bounce_mode       bounce averaged calculation: (0) OFF, (1) ON
       - dke_mode          equation solved: (0) FPE, (1) DKE
       - rdke              list of points for DKE calculations
       - display_mode      display mode parameter (full: 2, partial: 1, no: = 0)

 OUTPUTS:
   
       - ZYXXD_rf:         QL diffusion tensor elements for f0 
       - ZYXXD_rf_tp:      QL diffusion tensor elements for ftilde 
       - ZYXXF_rf_tp:      QL convection vector elements for ftilde 
       - gridindex_rf:     List of grid indices in resonant point vectors

 NOTE:

   the first characters of a variable name are related to its structure:
       - x_    1-D radial (rho)
       - X_    2-D momentum (p,mhu) 
       - Xx_   1-D momentum + 1-D radial (mhu,rho)
       - xy_   1-D radial + 1-D ray number (rho,b)
       - XX_   1-D radial + 2-D momentum (rho,p,mhu)

   refer to documentation for detailed description    

 by J. Decker (MIT/RLE) <jodecker@mit.edu> and Y. Peysson CEA/DRFC <yves.peysson@cea.fr>

0001 function [ZXYD_rf,ZXYF_rf,ZXYD_rf_tp,ZXYF_rf_tp,gridindex_rf] = rfdiff_dke_jd(dkeparam,display,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,mksa,radialDKE,waveparam);
0002 %
0003 % Calculates the RF quasilinear diffusion coefficient and the diffusion
0004 % tensor elements on a 3-D Grid: 2-D (p,mhu) in momentum; 1-D (rho) in space.
0005 % The coefficients associated with each ray are saved in the file
0006 % temp_DQL.mat while the resulting coefficients are returned in the
0007 % variables ZYXXD_rf,ZYXXD_rf_tp and ZYXXF_rf_tp.
0008 % Note that since the coefficients are saved, only the missing elements are
0009 % calculated in rfdiff_dke_jd while others are retrived. The coefficients
0010 % are finally weighted by the magnitude of the diffusion coefficient.
0011 %
0012 % INPUTS:
0013 %
0014 %       - waveparam:        RF wave parameters (See detailed description below in code)
0015 %                           waveparam is build by the function rfwave_dke_jd.m
0016 %       - nmax_rf:          maximum harmonic number to be included
0017 %       - save_id:          identification string, based on running date and time
0018 %       - temppath:         path to temporary folder
0019 %       - theta...xqtilde:  equilibrium properties
0020 %       - Xmhup...Xdmhupp:  pitch angle grid
0021 %       - Xpnpp...Xdpnpp:   momentum grid
0022 %       - Xgamma...:        relativistic factor
0023 %       - Xxlambda...:      normalized bounce time
0024 %       - betath_ref:       sqrt(Te_ref/mc2)
0025 %       - bounce_mode       bounce averaged calculation: (0) OFF, (1) ON
0026 %       - dke_mode          equation solved: (0) FPE, (1) DKE
0027 %       - rdke              list of points for DKE calculations
0028 %       - display_mode      display mode parameter (full: 2, partial: 1, no: = 0)
0029 %
0030 % OUTPUTS:
0031 %
0032 %       - ZYXXD_rf:         QL diffusion tensor elements for f0
0033 %       - ZYXXD_rf_tp:      QL diffusion tensor elements for ftilde
0034 %       - ZYXXF_rf_tp:      QL convection vector elements for ftilde
0035 %       - gridindex_rf:     List of grid indices in resonant point vectors
0036 %
0037 % NOTE:
0038 %
0039 %   the first characters of a variable name are related to its structure:
0040 %       - x_    1-D radial (rho)
0041 %       - X_    2-D momentum (p,mhu)
0042 %       - Xx_   1-D momentum + 1-D radial (mhu,rho)
0043 %       - xy_   1-D radial + 1-D ray number (rho,b)
0044 %       - XX_   1-D radial + 2-D momentum (rho,p,mhu)
0045 %
0046 %   refer to documentation for detailed description
0047 %
0048 % by J. Decker (MIT/RLE) <jodecker@mit.edu> and Y. Peysson CEA/DRFC <yves.peysson@cea.fr>
0049 %
0050 etime_rf = 0;
0051 time0 = clock;
0052 %
0053 [qe,me,mp,mn,e0] = pc_dke_yp;%Universal physics constant
0054 %
0055 rf_tol = eps('single');% threshold for calculation mask on deltares
0056 %
0057 n_rf_list = single(waveparam.n_rf_list);
0058 bounce_mode = dkeparam.bounce_mode;
0059 dke_mode = dkeparam.dke_mode;
0060 %
0061 display_mode = display.display_mode;
0062 %
0063 theta = equilDKE.mtheta;
0064 Xx = equilDKE.Xx;
0065 Xy = equilDKE.Xy;
0066 XBx = equilDKE.XBx;
0067 XBy = equilDKE.XBy;
0068 XBphi = equilDKE.XBphi;
0069 xB0 = equilDKE.xB0;
0070 Rp = equilDKE.Rp;
0071 Zp = equilDKE.Zp;
0072 ap = equilDKE.ap;
0073 xdV_2piRp_dke = equilDKE.xdV_2piRp_dke;
0074 %
0075 Xmhup = gridDKE.Xmhup;
0076 Xmhu = gridDKE.Xmhu;
0077 Xmhum = gridDKE.Xmhum;
0078 X1mmhu2p = gridDKE.X1mmhu2p;
0079 X1mmhu2 = gridDKE.X1mmhu2;
0080 X1mmhu2m = gridDKE.X1mmhu2m;
0081 Xpnpp = gridDKE.Xpnpp;
0082 Xpnp = gridDKE.Xpnp;
0083 Xpn = gridDKE.Xpn;
0084 Xpnm = gridDKE.Xpnm;
0085 Xpnmm = gridDKE.Xpnmm;
0086 %
0087 Xgammap = Zmomcoef.Xgammap;
0088 Xgamma = Zmomcoef.Xgamma;
0089 Xgammam = Zmomcoef.Xgammam;
0090 %
0091 xqtilde = Zbouncecoef.xqtilde;
0092 xqhat = Zbouncecoef.xqhat;
0093 Xxlambdap = Zbouncecoef.Xxlambda_p;
0094 Xxlambda = Zbouncecoef.Xxlambda;
0095 Xxlambdam = Zbouncecoef.Xxlambda_m;
0096 XXheavisidep = Zbouncecoef.XXheaviside_p;
0097 XXheaviside = Zbouncecoef.XXheaviside;
0098 XXheavisidem = Zbouncecoef.XXheaviside_m;
0099 %
0100 betath_ref = mksa.betath_ref;
0101 ne_ref = mksa.ne_ref;
0102 lnc_e_ref = mksa.lnc_e_ref;
0103 %
0104 rdke = radialDKE.r_dke;
0105 %
0106 Xsigmap = sign(Xmhup);
0107 Xsigma = sign(Xmhu);
0108 Xsigmam = sign(Xmhum);
0109 %
0110 pnmax = max(Xpn(:,1));
0111 %
0112 % RF wave parameters
0113 %
0114 biy_rf = uint32(waveparam.biy_rf);% subray the fragment belongs to
0115 bir_rf = uint32(waveparam.bir_rf);% FS the fragment belongs to
0116 bds_rf = single(waveparam.bds_rf);% ray fragement path length
0117 bomega_rf = single(waveparam.bomega_rf);% RF wave frequency
0118 bopt_rf = uint32(waveparam.bopt_rf);% options for Larmor radius effects and/or LH calc.
0119 %                               - 0:exact,
0120 %                               - 1:approx (kperp*rhoe<<1)
0121 %                               - 2: kperp*rhoe<<1 and old LH calc. (Karney, Shoucri)
0122 bthetab_rf = single(waveparam.bthetab_rf);% poloidal location along the ray (bounce_mode>0 only)
0123 bB_rf = single(waveparam.bB_rf);% B field along the ray
0124 bnpar_rf = single(waveparam.bNpar_rf);% n_par along the ray
0125 bdnpar_rf = single(waveparam.bdNpar_rf);% dn_par along the ray
0126 bnperp_rf = single(waveparam.bNperp_rf);%  n_perp along the ray
0127 bepolp_rf = single(waveparam.bepolp_rf);%  polarization plus along the ray
0128 bepolm_rf = single(waveparam.bepolm_rf);%  polarization minus along the ray
0129 bepolz_rf = single(waveparam.bepolz_rf);%  polarization parallel along the ray
0130 bphinorm_rf = single(waveparam.bphinorm_rf);%  normalized energy flow density along the ray
0131 %
0132 clear dkeparam display equilDKE gridDKE Zmomcoef Zbouncecoef mksa radialDKE waveparam
0133 %
0134 [npn,nmhu] = size(Xpn);
0135 nr = length(xB0);
0136 nr_dke = length(rdke);
0137 nb = length(biy_rf);
0138 ny = double(max([biy_rf,0]));
0139 nn_rf = length(n_rf_list);
0140 %
0141 gridindex_rf.npn = npn;
0142 gridindex_rf.nmhu = nmhu;
0143 gridindex_rf.nr = nr;
0144 gridindex_rf.nr_dke = nr_dke;
0145 gridindex_rf.ny = ny;
0146 gridindex_rf.nn_rf = nn_rf;
0147 %
0148 % Case where no RF wave is to be accounted for
0149 %
0150 if nb == 0,
0151     ZXYD_rf = [];
0152     ZXYD_rf_tp = [];
0153     ZXYF_rf_tp = [];
0154     %
0155     gridindex_rf.mask_p = zeros(0,0,'uint32');%mask of indices in the momentum grid
0156     gridindex_rf.mask_m = zeros(0,0,'uint32');%mask of indices in the pitch angle grid
0157     gridindex_rf.mask_r = zeros(0,0,'uint32');%mask of indices in the radial grid
0158     gridindex_rf.mask_y = zeros(0,0,'uint32');%mask of indices in the ray fragments list
0159     gridindex_rf.mask_n = zeros(0,0,'uint32');%mask of indices in the harmonic number list
0160     gridindex_rf.mask_dke = zeros(0,0,'uint32');%mask of indices in the DKE grid list
0161     %
0162     ZXYD_rf.parpar_ij = [];
0163     ZXYD_rf.ij = []; 
0164 %
0165     % Terms from Sp at l+1/2 (Notations: _imj -> (i,j+1+2);_ipj -> (i+1,j+1/2)
0166     %
0167     ZXYD_rf.pp_ipj = [];%momentum dynamics only
0168     ZXYD_rf.pp_imj = [];%momentum dynamics only
0169     ZXYD_rf.pm_ipj = [];%momentum dynamics only
0170     ZXYD_rf.pm_imj = [];%momentum dynamics only
0171     ZXYD_rf.pr_ipj = [];
0172     ZXYD_rf.pr_imj = [];   
0173     ZXYF_rf.p_ipj = [];%momentum dynamics only
0174     ZXYF_rf.p_imj = [];%momentum dynamics only
0175     %
0176     % Terms from Sksi at l+1/2 (Notations: _ijp -> (i+1/2,j+1);_ijm -> (i+1/2,j)
0177     %
0178     ZXYD_rf.mm_ijp = [];%momentum dynamics only
0179     ZXYD_rf.mm_ijm = [];%momentum dynamics only
0180     ZXYD_rf.mp_ijp = [];%momentum dynamics only
0181     ZXYD_rf.mp_ijm = [];%momentum dynamics only
0182     ZXYD_rf.mr_ijp = [];
0183     ZXYD_rf.mr_ijm = [];
0184     ZXYF_rf.m_ijp = [];%momentum dynamics only
0185     ZXYF_rf.m_ijm = [];%momentum dynamics only
0186     %
0187     % Terms from Spsi at (i+1/2,j+1/2) (Notations: _lm -> (l);_lp -> (l+1);
0188     %
0189     ZXYD_rf.rr_lm = [];%radial dynamics only
0190     ZXYD_rf.rr_lp = [];%radial dynamics only
0191     ZXYD_rf.rp_lm = [];
0192     ZXYD_rf.rp_lp = [];
0193     ZXYD_rf.rm_lm = [];
0194     ZXYD_rf.rm_lp = [];
0195     ZXYF_rf.r_lm = [];%radial dynamics only
0196     ZXYF_rf.r_lp = [];%radial dynamics only
0197     %
0198     % Terms used by the old scheme only for momentum cross-derivatives (M.
0199     % Shoucri and I. Shkarofsky,Comp. Phys. Comm., 82 (1994) 287)
0200     %
0201     ZXYD_rf.pp_ij = [];
0202     ZXYD_rf.pm_ij = [];
0203     ZXYD_rf.mp_ij = [];
0204     ZXYF_rf.p_ij = [];     
0205     %
0206     if dke_mode == 1,
0207         %
0208         ZXYD_rf_tp.pp_ippj = [];
0209         ZXYD_rf_tp.pp_ipj =  [];
0210         ZXYD_rf_tp.pp_ij =   [];
0211         ZXYD_rf_tp.pp_imj =  [];
0212         ZXYD_rf_tp.pp_immj = [];
0213         ZXYD_rf_tp.pm_ipj =  [];
0214         ZXYD_rf_tp.pm_ij =   [];
0215         ZXYD_rf_tp.pm_imj =  [];
0216         ZXYD_rf_tp.mp_ijp =  [];
0217         ZXYD_rf_tp.mp_ij =   [];
0218         ZXYD_rf_tp.mp_ijm =  [];
0219         ZXYD_rf_tp.mm_ijp =  [];
0220         ZXYD_rf_tp.mm_ijm =  [];
0221         %
0222         % ---------------- Terms for wave-induced radial transport by adjoint technique (P. Helander model) ----------------
0223         %
0224         ZXYD_rf_tp.mp1_ijp = [];
0225         ZXYD_rf_tp.mp1_ij =  [];
0226         ZXYD_rf_tp.mp1_ijm = [];
0227         ZXYD_rf_tp.mm1_ijp = [];
0228         ZXYD_rf_tp.mm1_ijm = [];
0229         %
0230         ZXYF_rf_tp.p_ippj =  [];
0231         ZXYF_rf_tp.p_ipj =   [];
0232         ZXYF_rf_tp.p_ij =    [];
0233         ZXYF_rf_tp.p_imj =   [];
0234         ZXYF_rf_tp.p_immj =  [];
0235         ZXYF_rf_tp.m_ijp =   [];
0236         ZXYF_rf_tp.m_ijm =   [];
0237     else
0238         ZXYD_rf_tp = NaN;
0239         ZXYF_rf_tp = NaN;
0240     end            
0241     %
0242     return
0243 end
0244 %
0245 %
0246 %
0247 [biy_rf_sort,ib_sort] = sort(biy_rf);
0248 bomega_rf_sort = bomega_rf(ib_sort);
0249 iy_sort = 1+[0,find(diff(biy_rf_sort) ~= 0)];
0250 if iy_sort ~= (1:ny),
0251     error('subray list is inconsistent')
0252 end
0253 %
0254 yomega_rf = bomega_rf_sort(iy_sort);
0255 %
0256 if length(theta) > 1,
0257     xbxb_rf = interp1(theta,Xx.',bthetab_rf);
0258     xbyb_rf = interp1(theta,Xy.',bthetab_rf);
0259     xbBxb_rf = interp1(theta,XBx.',bthetab_rf);
0260     xbByb_rf = interp1(theta,XBy.',bthetab_rf);
0261     xbBphib_rf = interp1(theta,XBphi.',bthetab_rf);
0262 else
0263     xbxb_rf = Xx;
0264     xbyb_rf = Xy;
0265     xbBxb_rf = XBx;
0266     xbByb_rf = XBy;
0267     xbBphib_rf = XBphi;
0268 end    
0269 %
0270 if nr > 1,
0271     xbxb_rf = xbxb_rf.';
0272     xbyb_rf = xbyb_rf.';
0273     xbBxb_rf = xbBxb_rf.';
0274     xbByb_rf = xbByb_rf.';
0275     xbBphib_rf = xbBphib_rf.';    
0276 end
0277 %
0278 xbrb_rf = sqrt(xbxb_rf.^2 + xbyb_rf.^2);% minor radius at theta = thetab_rf at the center of any FS
0279 xbBPb_rf = sqrt(xbBxb_rf.^2 + xbByb_rf.^2);% poloidal magnetic field at theta = thetab_rf at the center of any FS
0280 xbBb_rf = sqrt(xbBPb_rf.^2 + xbBphib_rf.^2);% magnetic field at theta = thetab_rf at the center of any FS
0281 xbcosb_rf = abs(xbBxb_rf.*xbyb_rf - xbByb_rf.*xbxb_rf)./xbrb_rf./xbBPb_rf;
0282 %
0283 xbmask_rf = uint32((0:nb-1)*nr) + bir_rf;
0284 bxb_rf = xbxb_rf(xbmask_rf);% horizontal poloidal position at theta = thetab_rf at the center of the FS
0285 bBPb_rf = xbBPb_rf(xbmask_rf);% poloidal magnetic field at theta = thetab_rf at the center of the FS
0286 bBb_rf = xbBb_rf(xbmask_rf);% magnetic field at theta = thetab_rf at the center of the FS
0287 brb_rf = xbrb_rf(xbmask_rf);% minor radius at theta = thetab_rf at the center of the FS
0288 bcosb_rf = xbcosb_rf(xbmask_rf);% cos alpha at theta = thetab_rf at the center of the FS
0289 %
0290 %if bounce_mode > 0,
0291     bB0 = single(xB0(bir_rf));
0292     %
0293     bpsib = bBb_rf./bB0; %Ratio of magnetic fields at theta = thetab_rf at the center of the FS
0294 %else
0295 %    bpsib = ones(1,nb); %Ratio of magnetic fields at theta = thetab_rf at the center of the FS
0296 %end
0297 %
0298 bomegaratio = qe*bB_rf./bomega_rf/me; % cyclotron frequency along the ray path
0299 %
0300 bkperprho = bnperp_rf*betath_ref./bomegaratio;
0301 bZbfac = bkperprho.*sqrt(bpsib);% momentum-independent factor in XBZb
0302 %
0303 % Spectral properties
0304 %
0305 bnpar_acc = imag(bnpar_rf);%accessibility condition
0306 bnpar_rf = real(bnpar_rf);%actual npar value
0307 %
0308 bnparmin_rf = NaN*ones(1,nb,'single');
0309 bnparmax_rf = NaN*ones(1,nb,'single');
0310 %
0311 bmask_gauss = (imag(bdnpar_rf) == 0);
0312 %
0313 bnparmin_rf(~bmask_gauss) = bnpar_rf(~bmask_gauss) - imag(bdnpar_rf(~bmask_gauss))/2;
0314 bnparmax_rf(~bmask_gauss) = bnpar_rf(~bmask_gauss) + imag(bdnpar_rf(~bmask_gauss))/2;
0315 %
0316 bnparmin_rf(bmask_gauss) = bnpar_rf(bmask_gauss) - bdnpar_rf(bmask_gauss).*sqrt(-log(rf_tol*sqrt(pi)*bdnpar_rf(bmask_gauss)));
0317 bnparmax_rf(bmask_gauss) = bnpar_rf(bmask_gauss) + bdnpar_rf(bmask_gauss).*sqrt(-log(rf_tol*sqrt(pi)*bdnpar_rf(bmask_gauss)));
0318 %
0319 bdnpar_rf = abs(bdnpar_rf);
0320 %
0321 % Diffusion coefficient (normalized to the incident power per unit length along the major radius in W/m)
0322 %
0323 wpe_ref = sqrt(ne_ref*qe^2/e0/me);
0324 %
0325 bphinorm_rf(bphinorm_rf == 0) = eps('single');%avoid singularities
0326 %
0327 bqtilde = single(xqtilde(bir_rf));
0328 bdV_2piRp_dke = single(xdV_2piRp_dke(bir_rf));
0329 %
0330 bFtilde = brb_rf.*bBb_rf./(ap*bBPb_rf.*bqtilde.*bcosb_rf);
0331 %
0332 if ~isempty(bds_rf),
0333     bD0_rf = 4*pi^2*bds_rf.*bFtilde./(me*lnc_e_ref*wpe_ref^2.*bomega_rf.*bphinorm_rf.*bdV_2piRp_dke);
0334 else
0335     %
0336     % note that in this case bphinorm_rf should be replaced by bphi_perp_rf
0337     %
0338     bqhat = single(xqhat(bir_rf));
0339     bFhat = brb_rf.*bB0./(ap*bBPb_rf.*bqhat.*bcosb_rf*(1 + bxb_rf/Rp));    
0340     bD0_rf = 4*pi^2.*bFhat.*bFtilde./(me*lnc_e_ref*wpe_ref^2.*bomega_rf.*bphinorm_rf*2*pi.*brb_rf);
0341 end
0342 %
0343 bD0_rf = bD0_rf./(1 + bxb_rf/Rp);
0344 %
0345 % Determination of (np*nmhu*nb*nr*n_rf)-space points to be accounted for
0346 %
0347 %
0348 % sum over harmonic numbers and rays
0349 %
0350 if display_mode > 1, h = waitbar(0);end
0351 %
0352 % mask calculation of point that must be calculated
0353 %
0354 Nmat = npn*nmhu*min([nb,100]);
0355 %
0356 mask_p = zeros(Nmat,1,'uint32');%mask of indices in the momentum grid
0357 mask_m = zeros(Nmat,1,'uint32');%mask of indices in the pitch angle grid
0358 mask_b = zeros(Nmat,1,'uint32');%mask of indices in the ray fragments
0359 mask_n = zeros(Nmat,1,'uint32');%mask of indices in the harmonic number list
0360 %
0361 iX = 1;
0362 iMat = 1;
0363 %
0364 for in = 1:nn_rf,
0365     %
0366     n_rf = n_rf_list(in);
0367     %
0368     % Preselection based on the momentum space localization of the resonance curves
0369     %
0370     % See J. Decker PoP 13, 112503 (2006) Eq. 32
0371     %
0372     byn_rf = n_rf*bomegaratio;
0373     %
0374     bmask_n1 = (bnpar_rf == 1);
0375     bmask_yn = (byn_rf >= sqrt(1 - min([bnpar_rf.^2;ones(1,nb)])));
0376     bpn_res = abs((abs(bnpar_rf).*byn_rf - sqrt(byn_rf.^2 + bnpar_rf.^2 - 1))./(1 - bnpar_rf.^2))/betath_ref;
0377     bpn_res_n1 = NaN(1,nb);
0378     bpn_res_n1(byn_rf ~= 0) = abs(1 - byn_rf(byn_rf ~= 0).^2)./(2*byn_rf(byn_rf ~= 0)*betath_ref);
0379     bpn_res(bmask_n1) = bpn_res_n1(bmask_n1);
0380     %
0381     bmask_n1_min = (bnparmin_rf == 1);
0382     bmask_yn_min = (byn_rf >= sqrt(1 - min([bnparmin_rf.^2;ones(1,nb)])));
0383     bpn_res_min = abs((abs(bnparmin_rf).*byn_rf - sqrt(byn_rf.^2 + bnparmin_rf.^2 - 1))./(1 - bnparmin_rf.^2))/betath_ref;
0384     bpn_res_n1_min = NaN(1,nb);
0385     bpn_res_n1_min(byn_rf ~= 0) = abs(1 - byn_rf(byn_rf ~= 0).^2)./(2*byn_rf(byn_rf ~= 0)*betath_ref);
0386     bpn_res_min(bmask_n1_min) = bpn_res_n1_min(bmask_n1_min);
0387     %
0388     bmask_n1_max = (bnparmax_rf == 1);
0389     bmask_yn_max = (byn_rf >= sqrt(1 - min([bnparmax_rf.^2;ones(1,nb)])));
0390     bpn_res_max = abs((abs(bnparmax_rf).*byn_rf - sqrt(byn_rf.^2 + bnparmax_rf.^2 - 1))./(1 - bnparmax_rf.^2))/betath_ref;
0391     bpn_res_n1_max = NaN(1,nb);
0392     bpn_res_n1_max(byn_rf ~= 0) = abs(1 - byn_rf(byn_rf ~= 0).^2)./(2*byn_rf(byn_rf ~= 0)*betath_ref);
0393     bpn_res_max(bmask_n1_max) = bpn_res_n1_max(bmask_n1_max);
0394     %
0395     bmask = find((bmask_yn & bpn_res < pnmax) | (bmask_yn_min & bpn_res_min < pnmax) | (bmask_yn_max & bpn_res_max < pnmax));
0396     %
0397     for ib = bmask,   
0398         %
0399         Xmaskb = (bpsib(ib).*X1mmhu2 <= 1); % condition thetab_rf < thetaT
0400         %
0401         Xmhub = Xmaskb.*Xsigma.*sqrt(1 - bpsib(ib).*X1mmhu2);%mhu at theta = thetab_rf at the center of the FS
0402         Xmhub(Xmhub == 0) = eps('single');
0403         %
0404         % Resonance Condition: resonant Npar
0405         %
0406         Xnparres_ij = (Xgamma - n_rf.*bomegaratio(ib))./Xpn./Xmhub/betath_ref;
0407         %
0408         % Resonance Condition: delta function (normalized to: int(deltares d_npar) = 1)
0409         %
0410         if bmask_gauss(ib),
0411             Xdeltares_ij = exp(-(Xnparres_ij - bnpar_rf(ib)).^2./bdnpar_rf(ib).^2)./bdnpar_rf(ib)/sqrt(pi);
0412         else
0413             Xdeltares_ij = (Xnparres_ij > bnparmin_rf(ib)).*(bnparmax_rf(ib) > Xnparres_ij)./bdnpar_rf(ib);
0414         end
0415         %
0416         %Xdeltares_ij(abs(Xnparres_ij) < bnpar_acc(ib)) = 0;%accessibility condition
0417         %
0418         X_index_ij = (Xdeltares_ij >= rf_tol) | (~XXheaviside(:,:,bir_rf(ib)) & fliplr(Xdeltares_ij) >= rf_tol);
0419         %
0420         % neighboring points are also accounted for
0421         %
0422         X_index_immj = X_index_ij;X_index_immj(2:npn,:) = X_index_ij(1:npn-1,:);
0423         X_index_ippj = X_index_ij;X_index_ippj(1:npn-1,:) = X_index_ij(2:npn,:);
0424         X_index_ijmm = X_index_ij;X_index_ijmm(:,2:nmhu) = X_index_ij(:,1:nmhu-1);
0425         X_index_ijpp = X_index_ij;X_index_ijpp(:,1:nmhu-1) = X_index_ij(:,2:nmhu);      
0426         %
0427         X_index = find(X_index_ij | X_index_immj | X_index_ippj | X_index_ijmm | X_index_ijpp);
0428         %
0429         m_index = ceil(X_index/npn);
0430         p_index = X_index - (m_index - 1)*npn;
0431         b_index = ib*ones(size(X_index));
0432         n_index = in*ones(size(X_index));
0433         %
0434         N_index = length(m_index);
0435         %
0436         mask_index = iX:iX + N_index - 1;
0437         %
0438         while iX + N_index - 1 > iMat*Nmat,
0439             %
0440             iMat = iMat + 1; 
0441             %
0442             mask_p = [mask_p;zeros(Nmat,1,'uint32')];%mask of indices in the momentum grid
0443             mask_m = [mask_m;zeros(Nmat,1,'uint32')];%mask of indices in the pitch angle grid
0444             mask_b = [mask_b;zeros(Nmat,1,'uint32')];%mask of indices in the ray fragments
0445             mask_n = [mask_n;zeros(Nmat,1,'uint32')];%mask of indices in the harmonic number list
0446             %
0447         end
0448         %
0449         mask_p(mask_index) = p_index;
0450         mask_m(mask_index) = m_index;
0451         mask_b(mask_index) = b_index;
0452         mask_n(mask_index) = n_index;
0453         %
0454         iX = iX + N_index;
0455     %
0456         if display_mode > 1, waitbar((ib + nb*(in - 1))/(nb*nn_rf),h);end
0457         %
0458     end
0459     %
0460 end
0461 %
0462 if display_mode >1, close(h); end
0463 %
0464 mask_0 = 1:uint32(iX - 1);
0465 %
0466 mask_p = mask_p(mask_0);%mask of indices in the momentum grid
0467 mask_m = mask_m(mask_0);%mask of indices in the pitch angle grid
0468 mask_b = mask_b(mask_0);%mask of indices in the ray fragments
0469 mask_n = mask_n(mask_0);%mask of indices in the harmonic number list
0470 %
0471 clear XXsigma XX1mmhu2 XXpn XXgamma XXpsib XXomegac XXmaskb XXmhub 
0472 clear XXomega_rf XXnparmin_rf XXnparmax_rf XXnpar_rf XXdnpar_rf
0473 clear XXmask_square  XXmask_gauss XXnparres_ij XXdeltares_ij 
0474 clear XX_index XX_index_ij XX_index_immj XX_index_ippj XX_index_ijmm XX_index_ijpp
0475 clear p_index m_index r_index b_index n_index
0476 clear mask_index mask_0 N_index N iX imat
0477 %
0478 Nmax = 2000000;
0479 %
0480 %SXB = size(mask_p);
0481 nXB = length(mask_p);
0482 %
0483 XYDrf.ipj = spalloc(nr*ny,npn*nmhu*nn_rf,0);
0484 XYDrf.ij =  spalloc(nr*ny,npn*nmhu*nn_rf,0);
0485 XYDrf.imj = spalloc(nr*ny,npn*nmhu*nn_rf,0);
0486 XYDrf.ijp = spalloc(nr*ny,npn*nmhu*nn_rf,0);
0487 XYDrf.ijm = spalloc(nr*ny,npn*nmhu*nn_rf,0);
0488 %
0489 % ispalloc = 1;
0490 %
0491 % XBDrf.ipj = zeros(SXB,'single');
0492 % XBDrf.ij = zeros(SXB,'single');
0493 % XBDrf.imj = zeros(SXB,'single');
0494 % XBDrf.ijp = zeros(SXB,'single');
0495 % XBDrf.ijm = zeros(SXB,'single');
0496 %
0497 if display_mode > 0,
0498     disp(['Number of DQL elements : ',num2str(nXB)])
0499 end
0500 %
0501 for ibloc = 1:ceil(nXB/Nmax),
0502     %
0503     bloc_mask = (ibloc - 1)*Nmax + 1:min([ibloc*Nmax,nXB]);
0504     %
0505     mask_p_bloc = mask_p(bloc_mask);
0506     mask_m_bloc = mask_m(bloc_mask);
0507     mask_b_bloc = mask_b(bloc_mask);
0508     mask_n_bloc = mask_n(bloc_mask);
0509     %
0510     SXB_bloc = size(mask_p_bloc);
0511     %
0512     Xmask = (mask_m_bloc - 1)*npn + mask_p_bloc;%mask in momentum 2-D matrices
0513     %
0514     % Vector construction
0515     %
0516     XBn_rf = n_rf_list(mask_n_bloc);
0517     if nn_rf > 1,
0518         XBn_rf = XBn_rf.';
0519     end
0520     %
0521     if nb > 1, 
0522         XByn = XBn_rf.*bomegaratio(mask_b_bloc).';
0523         XBpsib = bpsib(mask_b_bloc).';   
0524     else
0525         XByn = XBn_rf.*bomegaratio(mask_b_bloc);
0526         XBpsib = bpsib(mask_b_bloc);   
0527     end
0528     %
0529     XBpnp = single(Xpnp(Xmask));
0530     XBpn = single(Xpn(Xmask));
0531     XBpnm = single(Xpnm(Xmask));
0532     %
0533     XBgammap = single(Xgammap(Xmask));
0534     XBgamma = single(Xgamma(Xmask));
0535     XBgammam = single(Xgammam(Xmask));
0536     %
0537     XBsigmap = single(Xsigmap(Xmask));
0538     XBsigma = single(Xsigma(Xmask));
0539     XBsigmam = single(Xsigmam(Xmask));
0540     %
0541     XB1mmhu2p = single(X1mmhu2p(Xmask));
0542     XB1mmhu2 = single(X1mmhu2(Xmask));
0543     XB1mmhu2m = single(X1mmhu2m(Xmask));
0544     %
0545     XB1mmhu2p(XB1mmhu2p == 0) = eps('single');
0546     XB1mmhu2(XB1mmhu2 == 0) =  eps('single');
0547     XB1mmhu2m(XB1mmhu2m == 0) =  eps('single');
0548     %
0549     XBmaskb = (XBpsib.*XB1mmhu2p <= 1) & (XBpsib.*XB1mmhu2 <= 1) & (XBpsib.*XB1mmhu2m <= 1);
0550     %
0551     XBmhubp = XBmaskb.*XBsigmap.*sqrt(1 - XBpsib.*XB1mmhu2p);
0552     XBmhub = XBmaskb.*XBsigma.*sqrt(1 - XBpsib.*XB1mmhu2);
0553     XBmhubm = XBmaskb.*XBsigmam.*sqrt(1 - XBpsib.*XB1mmhu2m);
0554     %
0555     if dke_mode == 0,
0556         %
0557         clear XBsigmap XBsigma XBsigmam
0558         %
0559     end     
0560     %
0561     XBmhubp(XBmhubp == 0) = eps('single');
0562     XBmhub(XBmhub == 0) =  eps('single');
0563     XBmhubm(XBmhubm == 0) =  eps('single');
0564     %
0565     % Resonant Npar
0566     %
0567     XBnparres.ipj = (XBgammap - XByn)./XBpnp./XBmhub/betath_ref;
0568     XBnparres.ij = (XBgamma - XByn)./XBpn./XBmhub/betath_ref;
0569     XBnparres.imj = (XBgammam - XByn)./XBpnm./XBmhub/betath_ref;
0570     XBnparres.ijp = (XBgamma - XByn)./XBpn./XBmhubp/betath_ref;
0571     XBnparres.ijm = (XBgamma - XByn)./XBpn./XBmhubm/betath_ref;
0572     %
0573     clear XByn
0574     %
0575     % Resonance Condition: delta function (normalized to: int(deltares d_npar) = 1)
0576     %
0577     if nb > 1,
0578         XBmask_gauss = bmask_gauss(mask_b_bloc).';
0579         XBmask_square = ~bmask_gauss(mask_b_bloc).';
0580         %
0581         XBnparmin_rf = bnparmin_rf(mask_b_bloc).';
0582         XBnparmax_rf = bnparmax_rf(mask_b_bloc).';
0583         XBnpar_rf = bnpar_rf(mask_b_bloc).';
0584         XBdnpar_rf = bdnpar_rf(mask_b_bloc).';
0585         XBnpar_acc = bnpar_acc(mask_b_bloc).';
0586         %
0587         XBZbfac = bZbfac(mask_b_bloc).';
0588     else
0589         XBmask_gauss = bmask_gauss(mask_b_bloc);
0590         XBmask_square = ~bmask_gauss(mask_b_bloc);
0591         %
0592         XBnparmin_rf = bnparmin_rf(mask_b_bloc);
0593         XBnparmax_rf = bnparmax_rf(mask_b_bloc);
0594         XBnpar_rf = bnpar_rf(mask_b_bloc);
0595         XBdnpar_rf = bdnpar_rf(mask_b_bloc);
0596         XBnpar_acc = bnpar_acc(mask_b_bloc);
0597         %
0598         XBZbfac = bZbfac(mask_b_bloc);
0599     end    
0600     %
0601     XBdeltares.ipj = zeros(SXB_bloc,'single');
0602     XBdeltares.ij = zeros(SXB_bloc,'single');
0603     XBdeltares.imj = zeros(SXB_bloc,'single');
0604     XBdeltares.ijp = zeros(SXB_bloc,'single');
0605     XBdeltares.ijm = zeros(SXB_bloc,'single');
0606     %
0607     XBnpar_acc(XBn_rf > 0) = -Inf;%Accessibility Condition Only for Landau Damping
0608     XBnpar_acc(XBnpar_acc == 0) = -Inf;%Accessibility Condition not given
0609     %
0610     XBnorm = erfc((XBnpar_acc - abs(XBnpar_rf))./XBdnpar_rf)/2;
0611     XBnorm(XBnorm == 0) = rf_tol;
0612     %
0613     XBdeltares.ipj(XBmask_square) = (XBnparres.ipj(XBmask_square) > XBnparmin_rf(XBmask_square)).*(XBnparmax_rf(XBmask_square) > XBnparres.ipj(XBmask_square))./XBdnpar_rf(XBmask_square);
0614     XBdeltares.ipj(XBmask_gauss) = exp(-(XBnparres.ipj(XBmask_gauss) - XBnpar_rf(XBmask_gauss)).^2./XBdnpar_rf(XBmask_gauss).^2)./XBdnpar_rf(XBmask_gauss)./XBnorm(XBmask_gauss)/sqrt(pi);
0615     XBdeltares.ipj(abs(XBnparres.ipj) < XBnpar_acc) = 0;
0616     %
0617     XBdeltares.ij(XBmask_square) = (XBnparres.ij(XBmask_square) > XBnparmin_rf(XBmask_square)).*(XBnparmax_rf(XBmask_square) > XBnparres.ij(XBmask_square))./XBdnpar_rf(XBmask_square);
0618     XBdeltares.ij(XBmask_gauss) = exp(-(XBnparres.ij(XBmask_gauss) - XBnpar_rf(XBmask_gauss)).^2./XBdnpar_rf(XBmask_gauss).^2)./XBdnpar_rf(XBmask_gauss)./XBnorm(XBmask_gauss)/sqrt(pi);
0619     XBdeltares.ij(abs(XBnparres.ij) < XBnpar_acc) = 0;
0620     %
0621     XBdeltares.imj(XBmask_square) = (XBnparres.imj(XBmask_square) > XBnparmin_rf(XBmask_square)).*(XBnparmax_rf(XBmask_square) > XBnparres.imj(XBmask_square))./XBdnpar_rf(XBmask_square);
0622     XBdeltares.imj(XBmask_gauss) = exp(-(XBnparres.imj(XBmask_gauss) - XBnpar_rf(XBmask_gauss)).^2./XBdnpar_rf(XBmask_gauss).^2)./XBdnpar_rf(XBmask_gauss)./XBnorm(XBmask_gauss)/sqrt(pi);
0623     XBdeltares.imj(abs(XBnparres.imj) < XBnpar_acc) = 0;
0624     %
0625     XBdeltares.ijp(XBmask_square) = (XBnparres.ijp(XBmask_square) > XBnparmin_rf(XBmask_square)).*(XBnparmax_rf(XBmask_square) > XBnparres.ijp(XBmask_square))./XBdnpar_rf(XBmask_square);
0626     XBdeltares.ijp(XBmask_gauss) = exp(-(XBnparres.ijp(XBmask_gauss) - XBnpar_rf(XBmask_gauss)).^2./XBdnpar_rf(XBmask_gauss).^2)./XBdnpar_rf(XBmask_gauss)./XBnorm(XBmask_gauss)/sqrt(pi);
0627     XBdeltares.ijp(abs(XBnparres.ijp) < XBnpar_acc) = 0;
0628     %
0629     XBdeltares.ijm(XBmask_square) = (XBnparres.ijm(XBmask_square) > XBnparmin_rf(XBmask_square)).*(XBnparmax_rf(XBmask_square) > XBnparres.ijm(XBmask_square))./XBdnpar_rf(XBmask_square);
0630     XBdeltares.ijm(XBmask_gauss) = exp(-(XBnparres.ijm(XBmask_gauss) - XBnpar_rf(XBmask_gauss)).^2./XBdnpar_rf(XBmask_gauss).^2)./XBdnpar_rf(XBmask_gauss)./XBnorm(XBmask_gauss)/sqrt(pi);
0631     XBdeltares.ijm(abs(XBnparres.ijm) < XBnpar_acc) = 0;
0632     %
0633     clear XBnparres XBmask_square XBmask_gauss XBnparmin_rf XBnparmax_rf XBnpar_rf XBdnpar_rf XBnpar_acc
0634     %
0635     % FLR factor k_perp*rho_perp
0636     %
0637     XBZb.ipj = -XBZbfac.*XBpnp.*sqrt(XB1mmhu2);
0638     XBZb.ij = -XBZbfac.*XBpn.*sqrt(XB1mmhu2);
0639     XBZb.imj = -XBZbfac.*XBpnm.*sqrt(XB1mmhu2);
0640     XBZb.ijp = -XBZbfac.*XBpn.*sqrt(XB1mmhu2p);
0641     XBZb.ijm = -XBZbfac.*XBpn.*sqrt(XB1mmhu2m);
0642     %
0643     clear XBZbfac
0644     %
0645     % Polarization term
0646     %
0647     % Polarization vector is normalized to 1.
0648     %
0649     if nb > 1,
0650         XBepolp_rf = bepolp_rf(mask_b_bloc).'; 
0651         XBepolm_rf = bepolm_rf(mask_b_bloc).'; 
0652         XBepolz_rf = bepolz_rf(mask_b_bloc).'; 
0653         %
0654         XBlre = (bopt_rf(mask_b_bloc).' == 0); %(0): exact calculation for the Bessel functions, otherwise, limit k_perp*rho_perp << 1
0655     else
0656         XBepolp_rf = bepolp_rf(mask_b_bloc); 
0657         XBepolm_rf = bepolm_rf(mask_b_bloc); 
0658         XBepolz_rf = bepolz_rf(mask_b_bloc); 
0659         %
0660         XBlre = (bopt_rf(mask_b_bloc) == 0); %(0): exact calculation for the Bessel functions, otherwise, limit k_perp*rho_perp << 1
0661     end
0662     %
0663     XBmhuratiop = XBmhubp./sqrt(XBpsib.*XB1mmhu2p);
0664     XBmhuratio = XBmhub./sqrt(XBpsib.*XB1mmhu2);
0665     XBmhuratiom = XBmhubm./sqrt(XBpsib.*XB1mmhu2m);
0666     %
0667     clear XB1mmhu2p XB1mmhu2 XB1mmhu2m
0668     %
0669     XBpolb.ipj = XBepolp_rf.*besselj_dke_jd(XBn_rf + 1,XBZb.ipj,XBlre)/sqrt(2) + XBepolm_rf.*besselj_dke_jd(XBn_rf - 1,XBZb.ipj,XBlre)/sqrt(2)...
0670        - XBmhuratio.*XBepolz_rf.*besselj_dke_jd(XBn_rf,XBZb.ipj,XBlre);
0671     XBpolb.ij = XBepolp_rf.*besselj_dke_jd(XBn_rf + 1,XBZb.ij,XBlre)/sqrt(2) + XBepolm_rf.*besselj_dke_jd(XBn_rf - 1,XBZb.ij,XBlre)/sqrt(2)...
0672        - XBmhuratio.*XBepolz_rf.*besselj_dke_jd(XBn_rf,XBZb.ij,XBlre);
0673     XBpolb.imj = XBepolp_rf.*besselj_dke_jd(XBn_rf + 1,XBZb.imj,XBlre)/sqrt(2) + XBepolm_rf.*besselj_dke_jd(XBn_rf - 1,XBZb.imj,XBlre)/sqrt(2)...
0674        - XBmhuratio.*XBepolz_rf.*besselj_dke_jd(XBn_rf,XBZb.imj,XBlre);
0675     XBpolb.ijp = XBepolp_rf.*besselj_dke_jd(XBn_rf + 1,XBZb.ijp,XBlre)/sqrt(2) + XBepolm_rf.*besselj_dke_jd(XBn_rf - 1,XBZb.ijp,XBlre)/sqrt(2)...
0676        - XBmhuratiop.*XBepolz_rf.*besselj_dke_jd(XBn_rf,XBZb.ijp,XBlre);
0677     XBpolb.ijm = XBepolp_rf.*besselj_dke_jd(XBn_rf + 1,XBZb.ijm,XBlre)/sqrt(2) + XBepolm_rf.*besselj_dke_jd(XBn_rf - 1,XBZb.ijm,XBlre)/sqrt(2)...
0678        - XBmhuratiom.*XBepolz_rf.*besselj_dke_jd(XBn_rf,XBZb.ijm,XBlre);    
0679     %
0680     clear XBZb XBn_rf XBepolp_rf XBepolm_rf XBepolz_rf XBlre XBmhuratiop XBmhuratio XBmhuratiom 
0681     %
0682     XBrespol.ipj = XBmaskb.*XBdeltares.ipj.*abs(XBpolb.ipj).^2;
0683     XBrespol.ij = XBmaskb.*XBdeltares.ij.*abs(XBpolb.ij).^2;
0684     XBrespol.imj = XBmaskb.*XBdeltares.imj.*abs(XBpolb.imj).^2;
0685     XBrespol.ijp = XBmaskb.*XBdeltares.ijp.*abs(XBpolb.ijp).^2;
0686     XBrespol.ijm = XBmaskb.*XBdeltares.ijm.*abs(XBpolb.ijm).^2;
0687     %
0688     clear XBdeltares XBpolb XBmaskb
0689     %
0690     %Sum over sigma for trapped particles - case FP
0691     %
0692     if dke_mode == 1,
0693         %
0694         XBrespol_tp = XBrespol;
0695         %
0696     end
0697     %
0698     Yxmask = ((mask_n_bloc - 1)*nb + (mask_b_bloc - 1))*npn + mask_p_bloc;%mask in pnp-y-rho-n 4-D matrices
0699     %
0700     XBrespol_mhuflip.ipj = flipud(sparse(double(mask_m_bloc),double(Yxmask),double(XBrespol.ipj),nmhu,npn*nb*nn_rf));
0701     XBrespol_mhuflip.ij = flipud(sparse(double(mask_m_bloc),double(Yxmask),double(XBrespol.ij),nmhu,npn*nb*nn_rf));
0702     XBrespol_mhuflip.imj = flipud(sparse(double(mask_m_bloc),double(Yxmask),double(XBrespol.imj),nmhu,npn*nb*nn_rf));
0703     XBrespol_mhuflip.ijp = flipud(sparse(double(mask_m_bloc),double(Yxmask),double(XBrespol.ijm),nmhu,npn*nb*nn_rf));
0704     XBrespol_mhuflip.ijm = flipud(sparse(double(mask_m_bloc),double(Yxmask),double(XBrespol.ijp),nmhu,npn*nb*nn_rf));
0705     %
0706     XBmask = (Yxmask - 1)*nmhu + mask_m_bloc;
0707     %
0708     XBrespol_mhuflip.ipj = single(full(XBrespol_mhuflip.ipj(XBmask)));
0709     XBrespol_mhuflip.ij = single(full(XBrespol_mhuflip.ij(XBmask)));
0710     XBrespol_mhuflip.imj = single(full(XBrespol_mhuflip.imj(XBmask)));
0711     XBrespol_mhuflip.ijp = single(full(XBrespol_mhuflip.ijp(XBmask)));
0712     XBrespol_mhuflip.ijm = single(full(XBrespol_mhuflip.ijm(XBmask))); 
0713     %
0714     if dke_mode == 0,
0715         %
0716         clear Yxmask XBmask
0717         %
0718     end
0719     %
0720     if nb > 1,
0721         mask_r_bloc = bir_rf(mask_b_bloc).';
0722     else
0723         mask_r_bloc = bir_rf(mask_b_bloc);
0724     end    
0725     %
0726     XXmask = ((mask_r_bloc - 1)*nmhu + (mask_m_bloc - 1))*npn + mask_p_bloc;%mask in rho-momentum 3-D matrices
0727     %
0728     XBtrapped_p = ~XXheavisidep(XXmask);
0729     XBtrapped = ~XXheaviside(XXmask);
0730     XBtrapped_m = ~XXheavisidem(XXmask);     
0731     %
0732     clear XXmask
0733     %
0734     XBrespol.ipj(XBtrapped) = (XBrespol.ipj(XBtrapped) + XBrespol_mhuflip.ipj(XBtrapped))/2;
0735     XBrespol.ij(XBtrapped) = (XBrespol.ij(XBtrapped) + XBrespol_mhuflip.ij(XBtrapped))/2;
0736     XBrespol.imj(XBtrapped) = (XBrespol.imj(XBtrapped) + XBrespol_mhuflip.imj(XBtrapped))/2;
0737     XBrespol.ijp(XBtrapped_p) = (XBrespol.ijp(XBtrapped_p) + XBrespol_mhuflip.ijp(XBtrapped_p))/2;
0738     XBrespol.ijm(XBtrapped_m) = (XBrespol.ijm(XBtrapped_m) + XBrespol_mhuflip.ijm(XBtrapped_m))/2;
0739     %
0740     clear XBrespol_mhuflip
0741     %
0742     if dke_mode == 0,
0743         %
0744         clear XBtrapped_p XBtrapped XBtrapped_m
0745         %
0746     end
0747     %
0748     Xxmask = (mask_r_bloc - 1)*nmhu + mask_m_bloc;%mask in mhu-rho 2-D matrices
0749     %
0750     XBlambdap = Xxlambdap(Xxmask);
0751     XBlambda = Xxlambda(Xxmask);        
0752     XBlambdam = Xxlambdam(Xxmask);        
0753     %
0754     clear Xxmask 
0755     %
0756     XBlambdap(XBlambdap == 0) = eps('single');
0757     XBlambda(XBlambda == 0) =  eps('single');
0758     XBlambdam(XBlambdam == 0) =  eps('single');
0759     %
0760     XBmhup = single(Xmhup(Xmask));
0761     XBmhu = single(Xmhu(Xmask));
0762     XBmhum = single(Xmhum(Xmask));
0763     %
0764     clear Xmask 
0765     %
0766     if nb > 1,
0767         XBD0_rf = bD0_rf(mask_b_bloc).';
0768     else
0769         XBD0_rf = bD0_rf(mask_b_bloc);
0770     end
0771     %
0772     XBDfacp = XBD0_rf.*XBmhup.*XBmhup.*XBpsib./XBlambdap./XBmhubp./XBmhubp;
0773     XBDfac = XBD0_rf.*XBmhu.*XBmhu.*XBpsib./XBlambda./XBmhub./XBmhub;
0774     XBDfacm = XBD0_rf.*XBmhum.*XBmhum.*XBpsib./XBlambdam./XBmhubm./XBmhubm;
0775     %
0776     clear XBD0_rf
0777     %
0778     if dke_mode == 0,
0779         %
0780         clear XBpsib XBmhubp XBmhub XBmhubm XBlambdap XBlambda XBlambdam XBqratio
0781         %
0782     end     
0783     %
0784     %QL Diffusion Coefficient
0785     %
0786     if nb > 1,
0787         XBold = (bopt_rf(mask_b_bloc).' == 2); %(2): no 1/vpar dependence (old method, benchmark Karney Shoucri)
0788     else
0789         XBold = (bopt_rf(mask_b_bloc) == 2); %(2): no 1/vpar dependence (old method, benchmark Karney Shoucri)
0790     end
0791     %
0792     if max(XBold(:)) > 0
0793         if display_mode >= 1,info_dke_yp(1,'WARNING: Old technique for calculating the LH matrix coefficients, for comparaison with previous versions (Karney, Shoucri)');end
0794     end
0795     %
0796     XBDrf_bloc.ipj = zeros(SXB_bloc,'single');
0797     XBDrf_bloc.ij = zeros(SXB_bloc,'single');
0798     XBDrf_bloc.imj = zeros(SXB_bloc,'single');
0799     XBDrf_bloc.ijp = zeros(SXB_bloc,'single');
0800     XBDrf_bloc.ijm = zeros(SXB_bloc,'single');
0801     %
0802     XBDrf_bloc.ipj(XBold) = XBrespol.ipj(XBold).*XBDfac(XBold);
0803     XBDrf_bloc.ij(XBold) = XBrespol.ij(XBold).*XBDfac(XBold);
0804     XBDrf_bloc.imj(XBold) = XBrespol.imj(XBold).*XBDfac(XBold);
0805     XBDrf_bloc.ijp(XBold) = XBrespol.ijp(XBold).*XBDfacp(XBold);
0806     XBDrf_bloc.ijm(XBold) = XBrespol.ijm(XBold).*XBDfacm(XBold);
0807     %
0808     XBDrf_bloc.ipj(~XBold) = XBrespol.ipj(~XBold).*XBDfac(~XBold).*XBgammap(~XBold)./XBpnp(~XBold)./abs(XBmhu(~XBold));
0809     XBDrf_bloc.ij(~XBold) = XBrespol.ij(~XBold).*XBDfac(~XBold).*XBgamma(~XBold)./XBpn(~XBold)./abs(XBmhu(~XBold));
0810     XBDrf_bloc.imj(~XBold) = XBrespol.imj(~XBold).*XBDfac(~XBold).*XBgammam(~XBold)./XBpnm(~XBold)./abs(XBmhu(~XBold));
0811     XBDrf_bloc.ijp(~XBold) = XBrespol.ijp(~XBold).*XBDfacp(~XBold).*XBgamma(~XBold)./XBpn(~XBold)./abs(XBmhup(~XBold));
0812     XBDrf_bloc.ijm(~XBold) = XBrespol.ijm(~XBold).*XBDfacm(~XBold).*XBgamma(~XBold)./XBpn(~XBold)./abs(XBmhum(~XBold));
0813     %
0814     clear XBrespol XBDfacp XBDfac XBDfacm 
0815     clear XBgammap XBgamma XBgammam   
0816     clear XBmhup XBmhu XBmhum
0817     %
0818     if dke_mode == 0,
0819         %
0820         clear XBpnpp XBpnp XBpn XBpnm XBpnmm XBold
0821         %
0822     end
0823     %
0824 %     XBDrf.ipj(bloc_mask) = XBDrf_bloc.ipj;
0825 %     XBDrf.ij(bloc_mask) = XBDrf_bloc.ij;
0826 %     XBDrf.imj(bloc_mask) = XBDrf_bloc.imj;
0827 %     XBDrf.ijp(bloc_mask) = XBDrf_bloc.ijp;
0828 %     XBDrf.ijm(bloc_mask) = XBDrf_bloc.ijm;
0829     %
0830     %
0831     % Sum over all fragments within a flux surface
0832     %
0833     %
0834     if nb > 1,
0835         mask_y_bloc = biy_rf(mask_b_bloc).';
0836     else
0837         mask_y_bloc = biy_rf(mask_b_bloc);
0838         %
0839     end
0840     %
0841     xymask_bloc = double((mask_y_bloc - 1)*nr + mask_r_bloc);
0842     Xnmask_bloc = double(((mask_n_bloc - 1)*nmhu + mask_m_bloc - 1)*npn + mask_p_bloc);%mask in pnp-y-rho-n 4-D matrices
0843     %
0844     XYDrf.ipj = XYDrf.ipj + sparse(xymask_bloc,Xnmask_bloc,double(XBDrf_bloc.ipj),nr*ny,npn*nmhu*nn_rf);
0845     XYDrf.ij = XYDrf.ij + sparse(xymask_bloc,Xnmask_bloc,double(XBDrf_bloc.ij),nr*ny,npn*nmhu*nn_rf);
0846     XYDrf.imj = XYDrf.imj + sparse(xymask_bloc,Xnmask_bloc,double(XBDrf_bloc.imj),nr*ny,npn*nmhu*nn_rf);
0847     XYDrf.ijp = XYDrf.ijp + sparse(xymask_bloc,Xnmask_bloc,double(XBDrf_bloc.ijp),nr*ny,npn*nmhu*nn_rf);
0848     XYDrf.ijm = XYDrf.ijm + sparse(xymask_bloc,Xnmask_bloc,double(XBDrf_bloc.ijm),nr*ny,npn*nmhu*nn_rf);
0849     %
0850     clear Xnmask_bloc xymask_bloc% XBDrf
0851     %
0852 %     if length(find(XYDrf.ij)) > Nmax*ispalloc/2,
0853 %     end
0854     %
0855 end
0856 %
0857 clear XBDrf_bloc  XXheavisidep XXheaviside XXheavisidem mask_b
0858 %
0859 XYmask = find(XYDrf.ipj | XYDrf.ij | XYDrf.imj | XYDrf.ijp | XYDrf.ijm);
0860 %
0861 if nr*ny == 1,
0862     XYmask = XYmask.';
0863 end
0864 %
0865 mask_n = ceil(XYmask/(nmhu*npn*ny*nr));
0866 mask_m = ceil((XYmask - (mask_n - 1)*nmhu*npn*ny*nr)/(npn*ny*nr));
0867 mask_p = ceil((XYmask - ((mask_n - 1)*nmhu + mask_m - 1)*npn*ny*nr)/(ny*nr));
0868 mask_y = ceil((XYmask - (((mask_n - 1)*nmhu + mask_m - 1)*npn + mask_p - 1)*ny*nr)/nr);
0869 mask_r = XYmask - ((((mask_n - 1)*nmhu + mask_m - 1)*npn + mask_p - 1)*ny + mask_y - 1)*nr;
0870 %
0871 XYmask = uint32(XYmask);
0872 mask_n = uint32(mask_n);
0873 mask_m = uint32(mask_m);
0874 mask_p = uint32(mask_p);
0875 mask_y = uint32(mask_y);
0876 mask_r = uint32(mask_r);
0877 %
0878 XYDrf.ipj = single(full(XYDrf.ipj(XYmask)));
0879 XYDrf.ij = single(full(XYDrf.ij(XYmask)));
0880 XYDrf.imj = single(full(XYDrf.imj(XYmask)));
0881 XYDrf.ijp = single(full(XYDrf.ijp(XYmask)));
0882 XYDrf.ijm = single(full(XYDrf.ijm(XYmask))); 
0883 %
0884 if nr*ny == 1,
0885     XYDrf.ipj = XYDrf.ipj.';
0886     XYDrf.ij = XYDrf.ij.';
0887     XYDrf.imj = XYDrf.imj.';
0888     XYDrf.ijp = XYDrf.ijp.';
0889     XYDrf.ijm = XYDrf.ijm.';
0890 end
0891 %
0892 clear XYmask
0893 %
0894 % QL diffusion tensor elements
0895 %
0896 Xmask = (mask_m - 1)*npn + mask_p;
0897 %
0898 XY1mmhu2p = single(X1mmhu2p(Xmask));
0899 XY1mmhu2 = single(X1mmhu2(Xmask));
0900 XY1mmhu2m = single(X1mmhu2m(Xmask));
0901 %
0902 XY1mmhu2p(XY1mmhu2p == 0) = eps('single');
0903 XY1mmhu2(XY1mmhu2 == 0) =  eps('single');
0904 XY1mmhu2m(XY1mmhu2m == 0) =  eps('single');
0905 %
0906 XYgammap = single(Xgammap(Xmask));
0907 XYgamma = single(Xgamma(Xmask));
0908 XYgammam = single(Xgammam(Xmask));
0909 %
0910 XYmhup = single(Xmhup(Xmask));
0911 XYmhu = single(Xmhu(Xmask));
0912 XYmhum = single(Xmhum(Xmask));
0913 %
0914 XYomega_rf = yomega_rf(mask_y);
0915 if ny > 1,
0916     XYomega_rf = XYomega_rf.';
0917 end    
0918 XYB0 = single(xB0(mask_r));
0919 if nr > 1,
0920     XYB0 = XYB0.';
0921 end    
0922 XYomegaratio0 = qe*XYB0./XYomega_rf/me; % cyclotron frequency along the ray path
0923 %
0924 XYn = n_rf_list(mask_n);
0925 if nn_rf > 1,
0926     XYn = XYn.';
0927 end
0928 %
0929 XYyn0 = XYomegaratio0.*XYn;
0930 %
0931 clear XYomega_rf XYB0 XYomegaratio0
0932 %
0933 ZXYD_rf.pp_ipj = XY1mmhu2.*XYDrf.ipj;   
0934 ZXYD_rf.pp_imj = XY1mmhu2.*XYDrf.imj;
0935 ZXYD_rf.pm_ipj = - sqrt(XY1mmhu2).*(XY1mmhu2 - XYyn0./XYgammap).*XYDrf.ipj./XYmhu;
0936 ZXYD_rf.pm_imj = - sqrt(XY1mmhu2).*(XY1mmhu2 - XYyn0./XYgammam).*XYDrf.imj./XYmhu;
0937 ZXYD_rf.pr_ipj = zeros(size(XYmhu));
0938 ZXYD_rf.pr_imj = zeros(size(XYmhu));
0939 ZXYF_rf.p_ipj = zeros(size(XYmhu));
0940 ZXYF_rf.p_imj = zeros(size(XYmhu));
0941 %
0942 ZXYD_rf.mm_ijp = (XY1mmhu2p - XYyn0./XYgamma).^2.*XYDrf.ijp./XYmhup./XYmhup;
0943 ZXYD_rf.mm_ijm = (XY1mmhu2m - XYyn0./XYgamma).^2.*XYDrf.ijm./XYmhum./XYmhum;
0944 ZXYD_rf.mp_ijp = - sqrt(XY1mmhu2p).*(XY1mmhu2p - XYyn0./XYgamma).*XYDrf.ijp./XYmhup;
0945 ZXYD_rf.mp_ijm = - sqrt(XY1mmhu2m).*(XY1mmhu2m - XYyn0./XYgamma).*XYDrf.ijm./XYmhum;
0946 ZXYD_rf.mr_ijp = zeros(size(XYmhu));
0947 ZXYD_rf.mr_ijm = zeros(size(XYmhu));
0948 ZXYF_rf.m_ijp = zeros(size(XYmhu));
0949 ZXYF_rf.m_ijm = zeros(size(XYmhu));
0950 %
0951 ZXYD_rf.rr_lp = zeros(size(XYmhu));
0952 ZXYD_rf.rr_lm = zeros(size(XYmhu));
0953 ZXYD_rf.rp_lp = zeros(size(XYmhu));
0954 ZXYD_rf.rp_lm = zeros(size(XYmhu));
0955 ZXYD_rf.rm_lp = zeros(size(XYmhu));
0956 ZXYD_rf.rm_lm = zeros(size(XYmhu));
0957 ZXYF_rf.r_lp = zeros(size(XYmhu));
0958 ZXYF_rf.r_lm = zeros(size(XYmhu));
0959 %
0960 ZXYD_rf.pp_ij = XY1mmhu2.*XYDrf.ij;                                      
0961 ZXYD_rf.pm_ij = - sqrt(XY1mmhu2).*(XY1mmhu2 - XYyn0./XYgamma).*XYDrf.ij./XYmhu;
0962 ZXYD_rf.mp_ij = - sqrt(XY1mmhu2).*(XY1mmhu2 - XYyn0./XYgamma).*XYDrf.ij./XYmhu;
0963 %
0964 ZXYD_rf.parpar_ij = XY1mmhu2.*(1 - XYyn0./XYgamma).^2.*XYDrf.ij./XYmhu./XYmhu; 
0965 ZXYD_rf.ij = XYDrf.ij; 
0966 %
0967 clear XYDrf 
0968 %
0969 if dke_mode == 1,
0970     %
0971     % Sum over sigma for trapped particles - case DKE
0972     %
0973     XYrespol_tp_mhuflip.ippj = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.ippj.*XYsigma),nmhu,npn*nb*nr));
0974     XYrespol_tp_mhuflip.ipj = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.ipj.*XYsigma),nmhu,npn*nb*nr));
0975     XYrespol_tp_mhuflip.ij = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.ij.*XYsigma),nmhu,npn*nb*nr));
0976     XYrespol_tp_mhuflip.imj = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.imj.*XYsigma),nmhu,npn*nb*nr));
0977     XYrespol_tp_mhuflip.immj = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.immj.*XYsigma),nmhu,npn*nb*nr));
0978     XYrespol_tp_mhuflip.ijp = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.ijp.*XYsigmap),nmhu,npn*nb*nr));
0979     XYrespol_tp_mhuflip.ijm = flipud(sparse(double(mask_m),double(Yxmask),double(XYrespol_tp.ijm.*XYsigmam),nmhu,npn*nb*nr));
0980     %
0981     XYrespol_tp_mhuflip.ippj = XYsigma.*single(full(XYrespol_tp_mhuflip.ippj(XYmask)));
0982     XYrespol_tp_mhuflip.ipj = XYsigma.*single(full(XYrespol_tp_mhuflip.ipj(XYmask)));
0983     XYrespol_tp_mhuflip.ij = XYsigma.*single(full(XYrespol_tp_mhuflip.ij(XYmask)));
0984     XYrespol_tp_mhuflip.imj = XYsigma.*single(full(XYrespol_tp_mhuflip.imj(XYmask)));
0985     XYrespol_tp_mhuflip.immj = XYsigma.*single(full(XYrespol_tp_mhuflip.immj(XYmask)));
0986     XYrespol_tp_mhuflip.ijp = XYsigmap.*single(full(XYrespol_tp_mhuflip.ijp(XYmask)));
0987     XYrespol_tp_mhuflip.ijm = XYsigmam.*single(full(XYrespol_tp_mhuflip.ijm(XYmask))); 
0988     %
0989     clear Yxmask XYmask
0990     %
0991     XYrespol_tp.ippj(XYtrapped) = (XYrespol_tp.ippj(XYtrapped) + XYrespol_tp_mhuflip.ippj(XYtrapped))/2;
0992     XYrespol_tp.ipj(XYtrapped) = (XYrespol_tp.ipj(XYtrapped) + XYrespol_tp_mhuflip.ipj(XYtrapped))/2;
0993     XYrespol_tp.ij(XYtrapped) = (XYrespol_tp.ij(XYtrapped) + XYrespol_tp_mhuflip.ij(XYtrapped))/2;
0994     XYrespol_tp.imj(XYtrapped) = (XYrespol_tp.imj(XYtrapped) + XYrespol_tp_mhuflip.imj(XYtrapped))/2;
0995     XYrespol_tp.immj(XYtrapped) = (XYrespol_tp.immj(XYtrapped) + XYrespol_tp_mhuflip.immj(XYtrapped))/2;
0996     XYrespol_tp.ijp(XYtrapped_p) = (XYrespol_tp.ijp(XYtrapped_p) + XYrespol_tp_mhuflip.ijp(XYtrapped_p))/2;
0997     XYrespol_tp.ijm(XYtrapped_m) = (XYrespol_tp.ijm(XYtrapped_m) + XYrespol_tp_mhuflip.ijm(XYtrapped_m))/2;
0998     %
0999     clear XYtrapped_p XYtrapped XYtrapped_m XYrespol_tp_mhuflip
1000     %
1001     XYDfacp_D_tp = XYsigmap.*abs(XYmhup)./XYlambdap./XYmhubp;
1002     XYDfac_D_tp = XYsigma.*abs(XYmhu)./XYlambda./XYmhub;
1003     XYDfacm_D_tp = XYsigmam.*abs(XYmhum)./XYlambdam./XYmhubm;
1004     %
1005     XYDfacp_F_tp = (1 - XYpsib).*XYmhup.*XYmhup.*XYmhup./XYlambdap./XYmhubp./XYmhubp./XYmhubp;
1006     XYDfac_F_tp = (1 - XYpsib).*XYmhu.*XYmhu.*XYmhu./XYlambda./XYmhub./XYmhub./XYmhub;
1007     XYDfacm_F_tp = (1 - XYpsib).*XYmhum.*XYmhum.*XYmhum./XYlambdam./XYmhubm./XYmhubm./XYmhubm;
1008     %
1009     clear XYpsib XYsigmap XYsigma XYsigmam XYlambdap XYlambda XYlambdam XYqratio XYmhubp XYmhub XYmhubm 
1010     %
1011     %QL Diffusion Coefficient
1012     %
1013     XYDrf_D_tp.ippj = zeros(SXY,'single');
1014     XYDrf_D_tp.ipj = zeros(SXY,'single');
1015     XYDrf_D_tp.ij = zeros(SXY,'single');
1016     XYDrf_D_tp.imj = zeros(SXY,'single');
1017     XYDrf_D_tp.immj = zeros(SXY,'single');
1018     XYDrf_D_tp.ijp = zeros(SXY,'single');
1019     XYDrf_D_tp.ijm = zeros(SXY,'single');
1020     %
1021     XYDrf_D_tp.ippj(XYold) = XYrespol_tp.ippj(XYold).*XYDfac_D_tp(XYold);
1022     XYDrf_D_tp.ipj(XYold) = XYrespol_tp.ipj(XYold).*XYDfac_D_tp(XYold);
1023     XYDrf_D_tp.ij(XYold) = XYrespol_tp.ij(XYold).*XYDfac_D_tp(XYold);
1024     XYDrf_D_tp.imj(XYold) = XYrespol_tp.imj(XYold).*XYDfac_D_tp(XYold);
1025     XYDrf_D_tp.immj(XYold) = XYrespol_tp.immj(XYold).*XYDfac_D_tp(XYold);
1026     XYDrf_D_tp.ijp(XYold) = XYrespol_tp.ijp(XYold).*XYDfacp_D_tp(XYold);
1027     XYDrf_D_tp.ijm(XYold) = XYrespol_tp.ijm(XYold).*XYDfacm_D_tp(XYold);   
1028     %
1029     XYDrf_D_tp.ippj(~XYold) = XYrespol_tp.ippj(~XYold).*XYDfac_D_tp(~XYold).*XYgammapp(~XYold)./XYpnpp(~XYold)./abs(XYmhu(~XYold));
1030     XYDrf_D_tp.ipj(~XYold) = XYrespol_tp.ipj(~XYold).*XYDfac_D_tp(~XYold).*XYgammap(~XYold)./XYpnp(~XYold)./abs(XYmhu(~XYold));
1031     XYDrf_D_tp.ij(~XYold) = XYrespol_tp.ij(~XYold).*XYDfac_D_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhu(~XYold));
1032     XYDrf_D_tp.imj(~XYold) = XYrespol_tp.imj(~XYold).*XYDfac_D_tp(~XYold).*XYgammam(~XYold)./XYpnm(~XYold)./abs(XYmhu(~XYold));
1033     XYDrf_D_tp.immj(~XYold) = XYrespol_tp.immj(~XYold).*XYDfac_D_tp(~XYold).*XYgammamm(~XYold)./XYpnmm(~XYold)./abs(XYmhu(~XYold));
1034     XYDrf_D_tp.ijp(~XYold) = XYrespol_tp.ijp(~XYold).*XYDfacp_D_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhup(~XYold));
1035     XYDrf_D_tp.ijm(~XYold) = XYrespol_tp.ijm(~XYold).*XYDfacm_D_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhum(~XYold));   
1036     %
1037     clear XYDfacp_D_tp XYDfac_D_tp XYDfacm_D_tp
1038     %
1039     XYDrf_F_tp.ippj = zeros(SXY,'single');
1040     XYDrf_F_tp.ipj = zeros(SXY,'single');
1041     XYDrf_F_tp.ij = zeros(SXY,'single');
1042     XYDrf_F_tp.imj = zeros(SXY,'single');
1043     XYDrf_F_tp.immj = zeros(SXY,'single');
1044     XYDrf_F_tp.ijp = zeros(SXY,'single');
1045     XYDrf_F_tp.ijm = zeros(SXY,'single');
1046     %
1047     XYDrf_F_tp.ippj(XYold) = XYrespol_tp.ippj(XYold).*XYDfac_F_tp(XYold);
1048     XYDrf_F_tp.ipj(XYold) = XYrespol_tp.ipj(XYold).*XYDfac_F_tp(XYold);
1049     XYDrf_F_tp.ij(XYold) = XYrespol_tp.ij(XYold).*XYDfac_F_tp(XYold);
1050     XYDrf_F_tp.imj(XYold) = XYrespol_tp.imj(XYold).*XYDfac_F_tp(XYold);
1051     XYDrf_F_tp.immj(XYold) = XYrespol_tp.immj(XYold).*XYDfac_F_tp(XYold);
1052     XYDrf_F_tp.ijp(XYold) = XYrespol_tp.ijp(XYold).*XYDfacp_F_tp(XYold);
1053     XYDrf_F_tp.ijm(XYold) = XYrespol_tp.ijm(XYold).*XYDfacm_F_tp(XYold);       
1054     %
1055     XYDrf_F_tp.ippj(~XYold) = XYrespol_tp.ippj(~XYold).*XYDfac_F_tp(~XYold).*XYgammapp(~XYold)./XYpnpp(~XYold)./abs(XYmhu(~XYold));
1056     XYDrf_F_tp.ipj(~XYold) = XYrespol_tp.ipj(~XYold).*XYDfac_F_tp(~XYold).*XYgammap(~XYold)./XYpnp(~XYold)./abs(XYmhu(~XYold));
1057     XYDrf_F_tp.ij(~XYold) = XYrespol_tp.ij(~XYold).*XYDfac_F_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhu(~XYold));
1058     XYDrf_F_tp.imj(~XYold) = XYrespol_tp.imj(~XYold).*XYDfac_F_tp(~XYold).*XYgammam(~XYold)./XYpnm(~XYold)./abs(XYmhu(~XYold));
1059     XYDrf_F_tp.immj(~XYold) = XYrespol_tp.immj(~XYold).*XYDfac_F_tp(~XYold).*XYgammamm(~XYold)./XYpnmm(~XYold)./abs(XYmhu(~XYold));
1060     XYDrf_F_tp.ijp(~XYold) = XYrespol_tp.ijp(~XYold).*XYDfacp_F_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhup(~XYold));
1061     XYDrf_F_tp.ijm(~XYold) = XYrespol_tp.ijm(~XYold).*XYDfacm_F_tp(~XYold).*XYgamma(~XYold)./XYpn(~XYold)./abs(XYmhum(~XYold));       
1062     %
1063     % ---------------- Terms for wave-induced radial transport by adjoint technique (P. Helander model) ----------------
1064     %
1065     XYDrf_F1_tp.ijp = XYrespol_tp.ijp.*XYDfacp_F_tp.*XYgamma./XYpn;
1066     XYDrf_F1_tp.ij = XYrespol_tp.ij.*XYDfac_F_tp.*XYgamma./XYpn;
1067     XYDrf_F1_tp.ijm = XYrespol_tp.ijm.*XYDfacm_F_tp.*XYgamma./XYpn;       
1068     %
1069     clear XYrespol_tp XYDfacp_F_tp XYDfac_F_tp XYDfacm_F_tp
1070     clear XYgammapp XYgammap XYgamma XYgammam XYgammamm
1071     clear XYpnpp XYpnp XYpn XYpnm XYpnmm
1072     %
1073     % Mask for DKE points
1074     %
1075     mask_dke = any((1./rdke.')*(mask_r.') == 1).';%selects DKE points
1076     %
1077     %QL diffusion tensor elements
1078     %
1079     ZXYD_rf_tp.pp_ippj = XY1mmhu2(mask_dke).*XYDrf_D_tp.ippj(mask_dke);
1080     ZXYD_rf_tp.pp_ipj = XY1mmhu2(mask_dke).*XYDrf_D_tp.ipj(mask_dke);
1081     ZXYD_rf_tp.pp_ij = XY1mmhu2(mask_dke).*XYDrf_D_tp.ij(mask_dke);
1082     ZXYD_rf_tp.pp_imj = XY1mmhu2(mask_dke).*XYDrf_D_tp.imj(mask_dke);
1083     ZXYD_rf_tp.pp_immj = XY1mmhu2(mask_dke).*XYDrf_D_tp.immj(mask_dke);
1084     %
1085     ZXYD_rf_tp.pm_ipj = - sqrt(XY1mmhu2(mask_dke)).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_D_tp.ipj(mask_dke)./XYmhu(mask_dke);
1086     ZXYD_rf_tp.pm_ij = - sqrt(XY1mmhu2(mask_dke)).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_D_tp.ij(mask_dke)./XYmhu(mask_dke);
1087     ZXYD_rf_tp.pm_imj = - sqrt(XY1mmhu2(mask_dke)).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_D_tp.imj(mask_dke)./XYmhu(mask_dke);
1088     %
1089     ZXYD_rf_tp.mp_ijp = - sqrt(XY1mmhu2p(mask_dke)).*(XY1mmhu2p(mask_dke) - XYyn(mask_dke)).*XYDrf_D_tp.ijp(mask_dke)./XYmhup(mask_dke);
1090     ZXYD_rf_tp.mp_ij = - sqrt(XY1mmhu2(mask_dke)).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_D_tp.ij(mask_dke)./XYmhu(mask_dke);
1091     ZXYD_rf_tp.mp_ijm = - sqrt(XY1mmhu2m(mask_dke)).*(XY1mmhu2m (mask_dke)- XYyn(mask_dke)).*XYDrf_D_tp.ijm(mask_dke)./XYmhum(mask_dke);
1092     %
1093     ZXYD_rf_tp.mm_ijp = (XY1mmhu2p(mask_dke) - XYyn(mask_dke)).^2.*XYDrf_D_tp.ijp(mask_dke)./XYmhu2p(mask_dke);
1094     ZXYD_rf_tp.mm_ijm = (XY1mmhu2m(mask_dke) - XYyn(mask_dke)).^2.*XYDrf_D_tp.ijm(mask_dke)./XYmhu2m(mask_dke);
1095     %
1096     % ---------------- Terms for wave-induced radial transport by adjoint technique (P. Helander model) ----------------
1097     %
1098     ZXYD_rf_tp.mp1_ijp = - sqrt(XY1mmhu2p(mask_dke)).*(XY1mmhu2p(mask_dke) - XYyn(mask_dke)).*(XYDrf_D_tp.ijp(mask_dke) + XYDrf_F1_tp.ijp(mask_dke))./XYmhup(mask_dke);%XYDmp_rf_ijp
1099     ZXYD_rf_tp.mp1_ij = - sqrt(XY1mmhu2(mask_dke)).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*(XYDrf_D_tp.ij(mask_dke) + XYDrf_F1_tp.ij(mask_dke))./XYmhu(mask_dke);%XYDmp_rf_ij
1100     ZXYD_rf_tp.mp1_ijm = - sqrt(XY1mmhu2m(mask_dke)).*(XY1mmhu2m (mask_dke)- XYyn(mask_dke)).*(XYDrf_D_tp.ijm(mask_dke) + XYDrf_F1_tp.ijm(mask_dke))./XYmhum(mask_dke);%XYDmp_rf_ijm
1101     %
1102     ZXYD_rf_tp.mm1_ijp = (XY1mmhu2p(mask_dke) - XYyn(mask_dke)).^2.*(XYDrf_D_tp.ijp(mask_dke) + XYDrf_F1_tp.ijp(mask_dke))./XYmhu2p(mask_dke);%XYDmm_rf_ijp
1103     ZXYD_rf_tp.mm1_ijm = (XY1mmhu2m(mask_dke) - XYyn(mask_dke)).^2.*(XYDrf_D_tp.ijm(mask_dke) + XYDrf_F1_tp.ijm(mask_dke))./XYmhu2m(mask_dke);%XYDmm_rf_ijm
1104     %
1105     clear XYDrf_F1_tp XYDrf_D_tp
1106     %
1107     % ------------------------------------------------------------------------------------------------------------------
1108     %
1109     XYmhu5p = XYmhu2p./XYmhu2p./XYmhup;
1110     XYmhu4 = XYmhu2./XYmhu2;
1111     XYmhu5m = XYmhu2m./XYmhu2m./XYmhum;
1112     %
1113     ZXYF_rf_tp.p_ippj = - XY1mmhu2(mask_dke).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_F_tp.ippj(mask_dke)./XYmhu4(mask_dke)./XYpnpp(mask_dke); 
1114     ZXYF_rf_tp.p_ipj = - XY1mmhu2(mask_dke).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_F_tp.ipj(mask_dke)./XYmhu4(mask_dke)./XYpnp(mask_dke);
1115     ZXYF_rf_tp.p_ij = - XY1mmhu2(mask_dke).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_F_tp.ij(mask_dke)./XYmhu4(mask_dke)./XYpn(mask_dke); 
1116     ZXYF_rf_tp.p_imj = - XY1mmhu2(mask_dke).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_F_tp.imj(mask_dke)./XYmhu4(mask_dke)./XYpnm(mask_dke);
1117     ZXYF_rf_tp.p_immj = - XY1mmhu2(mask_dke).*(XY1mmhu2(mask_dke) - XYyn(mask_dke)).*XYDrf_F_tp.immj(mask_dke)./XYmhu4(mask_dke)./XYpnmm(mask_dke);
1118     %
1119     ZXYF_rf_tp.m_ijp = sqrt(XY1mmhu2p(mask_dke)).*(XY1mmhu2p(mask_dke) - XYyn(mask_dke)).^2.*XYDrf_F_tp.ijp(mask_dke)./XYmhu5p(mask_dke)./XYpn(mask_dke);
1120     ZXYF_rf_tp.m_ijm = sqrt(XY1mmhu2m(mask_dke)).*(XY1mmhu2m(mask_dke) - XYyn(mask_dke)).^2.*XYDrf_F_tp.ijm(mask_dke)./XYmhu5m(mask_dke)./XYpn(mask_dke);
1121     %
1122     clear XYDrf_F_tp XYmhu5p XYmhu4 XYmhu5m 
1123     %
1124     gridindex_rf.mask_dke = mask_dke;
1125     %
1126     clear mask_dke
1127     %
1128     %iXX = double(((mask_r(mask_dke) - 1)*nmhu + (mask_m(mask_dke) - 1))*npn + mask_p(mask_dke));
1129     %iyn = double((mask_n(mask_dke) - 1)*nb + mask_y(mask_dke));
1130     %
1131     %ZXXD_rf_tp.pp_ippj =   sparse(iyn,iXX,double(ZXYD_rf_tp.pp_ippj),nb*nn_rf,npn*nmhu*nr_dke);
1132     %ZXXD_rf_tp.pp_ipj =    sparse(iyn,iXX,double(ZXYD_rf_tp.pp_ipj ),nb*nn_rf,npn*nmhu*nr_dke);
1133     %ZXXD_rf_tp.pp_ij =     sparse(iyn,iXX,double(ZXYD_rf_tp.pp_ij  ),nb*nn_rf,npn*nmhu*nr_dke);
1134     %ZXXD_rf_tp.pp_imj =    sparse(iyn,iXX,double(ZXYD_rf_tp.pp_imj ),nb*nn_rf,npn*nmhu*nr_dke);
1135     %ZXXD_rf_tp.pp_immj =   sparse(iyn,iXX,double(ZXYD_rf_tp.pp_immj),nb*nn_rf,npn*nmhu*nr_dke);
1136     %ZXXD_rf_tp.pm_ipj =    sparse(iyn,iXX,double(ZXYD_rf_tp.pm_ipj ),nb*nn_rf,npn*nmhu*nr_dke);
1137     %ZXXD_rf_tp.pm_ij =     sparse(iyn,iXX,double(ZXYD_rf_tp.pm_ij  ),nb*nn_rf,npn*nmhu*nr_dke);
1138     %ZXXD_rf_tp.pm_imj =    sparse(iyn,iXX,double(ZXYD_rf_tp.pm_imj ),nb*nn_rf,npn*nmhu*nr_dke);
1139     %ZXXD_rf_tp.mp_ijp =    sparse(iyn,iXX,double(ZXYD_rf_tp.mp_ijp ),nb*nn_rf,npn*nmhu*nr_dke);
1140     %ZXXD_rf_tp.mp_ij =     sparse(iyn,iXX,double(ZXYD_rf_tp.mp_ij  ),nb*nn_rf,npn*nmhu*nr_dke);
1141     %ZXXD_rf_tp.mp_ijm =    sparse(iyn,iXX,double(ZXYD_rf_tp.mp_ijm ),nb*nn_rf,npn*nmhu*nr_dke);
1142     %ZXXD_rf_tp.mm_ijp =    sparse(iyn,iXX,double(ZXYD_rf_tp.mm_ijp ),nb*nn_rf,npn*nmhu*nr_dke);
1143     %ZXXD_rf_tp.mm_ijm =    sparse(iyn,iXX,double(ZXYD_rf_tp.mm_ijm ),nb*nn_rf,npn*nmhu*nr_dke);
1144     %
1145     %ZXXD_rf_tp.mp1_ijp =   sparse(iyn,iXX,double(ZXYD_rf_tp.mp1_ijp),nb*nn_rf,npn*nmhu*nr_dke);
1146     %ZXXD_rf_tp.mp1_ij =    sparse(iyn,iXX,double(ZXYD_rf_tp.mp1_ij ),nb*nn_rf,npn*nmhu*nr_dke);
1147     %ZXXD_rf_tp.mp1_ijm =   sparse(iyn,iXX,double(ZXYD_rf_tp.mp1_ijm),nb*nn_rf,npn*nmhu*nr_dke);
1148     %ZXXD_rf_tp.mm1_ijp =   sparse(iyn,iXX,double(ZXYD_rf_tp.mm1_ijp),nb*nn_rf,npn*nmhu*nr_dke);
1149     %ZXXD_rf_tp.mm1_ijm =   sparse(iyn,iXX,double(ZXYD_rf_tp.mm1_ijm),nb*nn_rf,npn*nmhu*nr_dke);
1150     %
1151     %ZXXF_rf_tp.p_ippj =    sparse(iyn,iXX,double(ZXYF_rf_tp.p_ippj ),nb*nn_rf,npn*nmhu*nr_dke);
1152     %ZXXF_rf_tp.p_ipj =     sparse(iyn,iXX,double(ZXYF_rf_tp.p_ipj  ),nb*nn_rf,npn*nmhu*nr_dke);
1153     %ZXXF_rf_tp.p_ij =      sparse(iyn,iXX,double(ZXYF_rf_tp.p_ij   ),nb*nn_rf,npn*nmhu*nr_dke);
1154     %ZXXF_rf_tp.p_imj =     sparse(iyn,iXX,double(ZXYF_rf_tp.p_imj  ),nb*nn_rf,npn*nmhu*nr_dke);
1155     %ZXXF_rf_tp.p_immj =    sparse(iyn,iXX,double(ZXYF_rf_tp.p_immj ),nb*nn_rf,npn*nmhu*nr_dke);
1156     %ZXXF_rf_tp.m_ijp =     sparse(iyn,iXX,double(ZXYF_rf_tp.m_ijp  ),nb*nn_rf,npn*nmhu*nr_dke);
1157     %ZXXF_rf_tp.m_ijm =     sparse(iyn,iXX,double(ZXYD_rf_tp.m_ijm  ),nb*nn_rf,npn*nmhu*nr_dke);
1158     %
1159 else
1160     ZXYD_rf_tp = NaN;
1161     ZXYF_rf_tp = NaN;
1162 end            
1163 %
1164 clear XYgammap XYgamma XYgammam
1165 clear XYmhup XYmhu XYmhum XY1mmhu2p XY1mmhu2 XY1mmhu2m XYyn0 XYold
1166 %
1167 gridindex_rf.mask_p = mask_p;
1168 gridindex_rf.mask_m = mask_m;
1169 gridindex_rf.mask_y = mask_y;
1170 gridindex_rf.mask_r = mask_r;
1171 gridindex_rf.mask_n = mask_n;
1172 %
1173 %iXX = double(((mask_r - 1)*nmhu + (mask_m - 1))*npn + mask_p);
1174 %iyn = double((mask_n - 1)*nb + mask_y);
1175 %
1176 %clear mask_r mask_m mask_p mask_n mask_y
1177 %
1178 %ZXXD_rf.pp_ippj =  sparse(iyn,iXX,double(ZXYD_rf.pp_ippj   ),nb*nn_rf,npn*nmhu*nr);
1179 %ZXXD_rf.pp_ipj =   sparse(iyn,iXX,double(ZXYD_rf.pp_ipj    ),nb*nn_rf,npn*nmhu*nr);
1180 %ZXXD_rf.pp_ij =    sparse(iyn,iXX,double(ZXYD_rf.pp_ij     ),nb*nn_rf,npn*nmhu*nr);
1181 %ZXXD_rf.pp_imj =   sparse(iyn,iXX,double(ZXYD_rf.pp_imj    ),nb*nn_rf,npn*nmhu*nr);
1182 %ZXXD_rf.pp_immj =  sparse(iyn,iXX,double(ZXYD_rf.pp_immj   ),nb*nn_rf,npn*nmhu*nr);
1183 %ZXXD_rf.pm_ipj =   sparse(iyn,iXX,double(ZXYD_rf.pm_ipj    ),nb*nn_rf,npn*nmhu*nr);
1184 %ZXXD_rf.pm_ij =    sparse(iyn,iXX,double(ZXYD_rf.pm_ij     ),nb*nn_rf,npn*nmhu*nr);
1185 %ZXXD_rf.pm_imj =   sparse(iyn,iXX,double(ZXYD_rf.pm_imj    ),nb*nn_rf,npn*nmhu*nr);
1186 %ZXXD_rf.mp_ijp =   sparse(iyn,iXX,double(ZXYD_rf.mp_ijp    ),nb*nn_rf,npn*nmhu*nr);
1187 %ZXXD_rf.mp_ij =    sparse(iyn,iXX,double(ZXYD_rf.mp_ij     ),nb*nn_rf,npn*nmhu*nr);
1188 %ZXXD_rf.mp_ijm =   sparse(iyn,iXX,double(ZXYD_rf.mp_ijm    ),nb*nn_rf,npn*nmhu*nr);
1189 %ZXXD_rf.mm_ijp =   sparse(iyn,iXX,double(ZXYD_rf.mm_ijp    ),nb*nn_rf,npn*nmhu*nr);
1190 %ZXXD_rf.mm_ijm =   sparse(iyn,iXX,double(ZXYD_rf.mm_ijm    ),nb*nn_rf,npn*nmhu*nr);
1191 %
1192 %
1193 etime_rf = etime_rf + etime(clock,time0);
1194 if display_mode >= 1,info_dke_yp(2,['RF flux coefficients calculations done in ',num2str(etime_rf),' (s)']);end