fmoments_dke_yp

PURPOSE ^

SYNOPSIS ^

function [ZP0,Zcurr,Znorm,ZXXS,XxRR_fsav,xMRR_flux,xMRR_tau,xMRR_power_flux,xMRR_power_tau,xRRm_fsav,xRRp_fsav,xne_norm_out,xTe_norm_out,xft_fsav,xfteff_fsav,Zmom] =fmoments_dke_yp(dkepath,dkeparam,display,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi,waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,XXf0,XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp,gridindex_rf,ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,XXRintmask,XXILor)

DESCRIPTION ^

   Calculate moments of the electron distribution function.

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

 Calculation of the various moments of the distribution function (density, current, power, bremsstrahlung)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ZP0,Zcurr,Znorm,ZXXS,...
0002     XxRR_fsav,xMRR_flux,xMRR_tau,xMRR_power_flux,xMRR_power_tau,xRRm_fsav,xRRp_fsav,...
0003     xne_norm_out,xTe_norm_out,xft_fsav,xfteff_fsav,Zmom] = ...
0004     fmoments_dke_yp(dkepath,dkeparam,display,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi,waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,...
0005                     XXf0,XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,...
0006                     ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp,gridindex_rf,...
0007                     ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,XXRintmask,XXILor)
0008 %
0009 %   Calculate moments of the electron distribution function.
0010 %
0011 %   by Y. Peysson (CEA-DRFC) (yves.peysson@cea.fr) and J. Decker (CEA-DRFC)   (joan.decker@cea.fr)
0012 %
0013 % Calculation of the various moments of the distribution function (density, current, power, bremsstrahlung)
0014 %
0015 time0 = clock;
0016 dke_mode = dkeparam.dke_mode;
0017 ZP0.n_rf_list = waveparam.n_rf_list;
0018 %
0019 Xpnm = gridDKE.Xpnm;
0020 Xpnp = gridDKE.Xpnp;
0021 X1mmhu2 = gridDKE.X1mmhu2;
0022 Xdpnm = gridDKE.Xdpnm;
0023 Xdpn = gridDKE.Xdpn;
0024 Xdmhu = gridDKE.Xdmhu;
0025 %
0026 XXdpnm = gridDKE.XXdpnm;
0027 XXdpnp = gridDKE.XXdpnp;
0028 XXdmhum = gridDKE.XXdmhum;
0029 XXdmhup = gridDKE.XXdmhup;
0030 %
0031 Xgammam = Zmomcoef.Xgammam;
0032 Xgammap = Zmomcoef.Xgammap;
0033 %
0034 XXlambda = Zbouncecoef.XXlambda;
0035 xqtilde = Zbouncecoef.xqtilde;
0036 xqhat =  Zbouncecoef.xqhat;
0037 %
0038 if nargout > 1,
0039     %
0040     Bax = equilDKE.Bax;
0041     xx0 = equilDKE.xx0;
0042     xBT0 = equilDKE.xBT0;
0043     xB0 = equilDKE.xB0;
0044     %
0045     dpnm = gridDKE.dpnm;
0046     dpn = gridDKE.dpn;
0047     dpnp = gridDKE.dpnp;
0048     dmhu = gridDKE.dmhu;
0049     dmhum = gridDKE.dmhum;
0050     %
0051     X1mmhu2m = gridDKE.X1mmhu2m;
0052     X1mmhu2p = gridDKE.X1mmhu2p;
0053     Xpn = gridDKE.Xpn;
0054     Xmhu = gridDKE.Xmhu;
0055     %
0056     XXpnm = gridDKE.XXpnm;
0057     XXpnp = gridDKE.XXpnp;
0058     XXpn = gridDKE.XXpn;
0059     XX1mmhu2 = gridDKE.XX1mmhu2;
0060     XX1mmhu2m = gridDKE.XX1mmhu2m;
0061     XX1mmhu2p = gridDKE.XX1mmhu2p;
0062     %
0063     XXgammam = Zmomcoef.XXgammam;
0064     XXgammap = Zmomcoef.XXgammap;
0065     v2 = Zmomcoef.v2;
0066     %
0067     XXRlambda_m = Zbouncecoef.XXRlambda_m;
0068     %
0069 end    
0070 %
0071 XXdeltap_imj = Zdelta.XXdeltap_imj; 
0072 XXdeltap_ipj = Zdelta.XXdeltap_ipj;
0073 %
0074 XXdeltam_ijm = Zdelta.XXdeltam_ijm; 
0075 XXdeltam_ijp = Zdelta.XXdeltam_ijp;
0076 %
0077 npn = gridindex_rf.npn;
0078 nmhu = gridindex_rf.nmhu;
0079 nr = gridindex_rf.nr;
0080 nr_dke = gridindex_rf.nr_dke;
0081 ny = gridindex_rf.ny;
0082 nn_rf = gridindex_rf.nn_rf;
0083 rdke = gridDKE.rdke;
0084 %
0085 % Derivatives of the distibution function
0086 %
0087 indexi = cumsum(ones(npn,1));
0088 indexj = cumsum(ones(nmhu,1));
0089 %
0090 indexim = indexi;
0091 indexip = indexi;
0092 %
0093 indexjm = indexj;
0094 indexjp = indexj;
0095 %
0096 indexim(2:npn) = indexi(1:npn-1);
0097 %indexim(1) = indexi(1); since df/dp(p==0) = 0;
0098 %
0099 indexip(1:npn-1) = indexi(2:npn);
0100 %indexip(npn) = indexi(npn); f(p>pmax) = 0 enforced later;
0101 %
0102 indexjm(2:nmhu) = indexj(1:nmhu-1);
0103 %indexjm(1) = indexj(1); since df/dm(m==-1) = 0;
0104 %
0105 indexjp(1:nmhu-1) = indexj(2:nmhu);
0106 %indexjp(nmhu) = indexj(nmhu); since df/dm(m==1) = 0;
0107 %
0108 XXf0_ij = XXf0;
0109 %
0110 XXf0_immj = XXf0_ij(indexim,:,:);
0111 XXf0_ippj = XXf0_ij(indexip,:,:);XXf0_ippj(npn,:,:) = 0;
0112 %
0113 XXf0_imj = (1 - XXdeltap_imj).*XXf0_ij + XXdeltap_imj.*XXf0_immj;
0114 XXf0_ipj = (1 - XXdeltap_ipj).*XXf0_ippj + XXdeltap_ipj.*XXf0_ij;
0115 %
0116 XXf0_imjmm = XXf0_imj(:,indexjm,:);
0117 XXf0_imjpp = XXf0_imj(:,indexjp,:);
0118 %
0119 XXf0_ipjmm = XXf0_ipj(:,indexjm,:);
0120 XXf0_ipjpp = XXf0_ipj(:,indexjp,:);
0121 %
0122 XXdf0dpn_imj = (XXf0_ij - XXf0_immj)./XXdpnm;
0123 XXdf0dpn_ipj = (XXf0_ippj - XXf0_ij)./XXdpnp;
0124 %
0125 XXdf0dmhu_imj = (XXf0_imjpp - XXf0_imjmm)./(XXdmhum + XXdmhup);
0126 XXdf0dmhu_ipj = (XXf0_ipjpp - XXf0_ipjmm)./(XXdmhum + XXdmhup);
0127 %
0128 if dke_mode == 1,
0129     %
0130     XXf0_tp_ij = XXf0_tp;
0131     %
0132     XXf0_tp_immj = XXf0_tp_ij(indexim,:,:);
0133     XXf0_tp_ippj = XXf0_tp_ij(indexip,:,:);XXf0_tp_ippj(npn,:,:) = 0;
0134     %
0135     XXf0_tp_imj = (1 - XXdeltap_imj).*XXf0_tp_ij + XXdeltap_imj.*XXf0_tp_immj;
0136     XXf0_tp_ipj = (1 - XXdeltap_ipj).*XXf0_tp_ippj + XXdeltap_ipj.*XXf0_tp_ij;
0137     %
0138     XXf0_tp_imjmm = XXf0_tp_imj(:,indexjm,:);
0139     XXf0_tp_imjpp = XXf0_tp_imj(:,indexjp,:);
0140     %
0141     XXf0_tp_ipjmm = XXf0_tp_ipj(:,indexjm,:);
0142     XXf0_tp_ipjpp = XXf0_tp_ipj(:,indexjp,:);
0143     %
0144     XXdf0dpn_tp_imj = (XXf0_tp_ij - XXf0_tp_immj)./XXdpnm;
0145     XXdf0dpn_tp_ipj = (XXf0_tp_ippj - XXf0_tp_ij)./XXdpnp;
0146     %
0147     XXdf0dmhu_tp_imj = (XXf0_tp_imjpp - XXf0_tp_imjmm)./(XXdmhum + XXdmhup);
0148     XXdf0dmhu_tp_ipj = (XXf0_tp_ipjpp - XXf0_tp_ipjmm)./(XXdmhum + XXdmhup);
0149     %
0150     %
0151     XXf0_g_ij = XXf0_g;
0152     %
0153     XXf0_g_immj = XXf0_g_ij(indexim,:,:);
0154     XXf0_g_ippj = XXf0_g_ij(indexip,:,:);XXf0_g_ippj(npn,:,:) = 0;
0155     %
0156     XXf0_g_imj = (1 - XXdeltap_imj).*XXf0_g_ij + XXdeltap_imj.*XXf0_g_immj;
0157     XXf0_g_ipj = (1 - XXdeltap_ipj).*XXf0_g_ippj + XXdeltap_ipj.*XXf0_g_ij;
0158     %
0159     XXf0_g_imjmm = XXf0_g_imj(:,indexjm,:);
0160     XXf0_g_imjpp = XXf0_g_imj(:,indexjp,:);
0161     %
0162     XXf0_g_ipjmm = XXf0_g_ipj(:,indexjm,:);
0163     XXf0_g_ipjpp = XXf0_g_ipj(:,indexjp,:);
0164     %
0165     XXdf0dpn_g_imj = (XXf0_g_ij - XXf0_g_immj)./XXdpnm;
0166     XXdf0dpn_g_ipj = (XXf0_g_ippj - XXf0_g_ij)./XXdpnp;
0167     %
0168     XXdf0dmhu_g_imj = (XXf0_g_imjpp - XXf0_g_imjmm)./(XXdmhum + XXdmhup);
0169     XXdf0dmhu_g_ipj = (XXf0_g_ipjpp - XXf0_g_ipjmm)./(XXdmhum + XXdmhup);
0170     %
0171 end
0172 %
0173 if nargout > 1,
0174     %
0175     XXf0_ijmm = XXf0_ij(:,indexjm,:);
0176     XXf0_ijpp = XXf0_ij(:,indexjp,:);
0177     %
0178     XXf0_ijm = (1 - XXdeltam_ijm).*XXf0_ij + XXdeltam_ijm.*XXf0_ijmm;
0179     XXf0_ijp = (1 - XXdeltam_ijp).*XXf0_ijpp + XXdeltam_ijp.*XXf0_ij;
0180     %
0181     XXf0_immjm = XXf0_ijm(indexim,:,:);
0182     XXf0_ippjm = XXf0_ijm(indexip,:,:);XXf0_ippjm(npn,:,:) = 0;
0183     %
0184     XXf0_immjp = XXf0_ijp(indexim,:,:);
0185     XXf0_ippjp = XXf0_ijp(indexip,:,:);XXf0_ippjp(npn,:,:) = 0;
0186     %
0187     XXdf0dmhu_ijm = (XXf0_ij - XXf0_ijmm)./XXdmhum;
0188     XXdf0dmhu_ijp = (XXf0_ijpp - XXf0_ij)./XXdmhup;
0189     %
0190     XXdf0dpn_ijm = (XXf0_ippjm - XXf0_immjm)./(XXdpnm + XXdpnp);
0191     XXdf0dpn_ijp = (XXf0_ippjp - XXf0_immjp)./(XXdpnm + XXdpnp);
0192     %
0193     if dke_mode == 1,
0194         %
0195         XXf0_tp_ijmm = XXf0_tp_ij(:,indexjm,:);
0196         XXf0_tp_ijpp = XXf0_tp_ij(:,indexjp,:);
0197         %
0198         XXf0_tp_ijm = (1 - XXdeltam_ijm).*XXf0_tp_ij + XXdeltam_ijm.*XXf0_tp_ijmm;
0199         XXf0_tp_ijp = (1 - XXdeltam_ijp).*XXf0_tp_ippj + XXdeltam_ijp.*XXf0_tp_ij;
0200         %
0201         XXf0_tp_immjm = XXf0_tp_ijm(indexim,:,:);
0202         XXf0_tp_ippjm = XXf0_tp_ijm(indexip,:,:);XXf0_tp_ippjm(npn,:,:) = 0;
0203         %
0204         XXf0_tp_immjp = XXf0_tp_ijp(indexim,:,:);
0205         XXf0_tp_ippjp = XXf0_tp_ijp(indexip,:,:);XXf0_tp_ippjp(npn,:,:) = 0;
0206         %
0207         XXdf0dmhu_tp_ijm = (XXf0_tp_ij - XXf0_tp_ijmm)./XXdmhum;
0208         XXdf0dmhu_tp_ijp = (XXf0_tp_ijpp - XXf0_tp_ij)./XXdmhup;
0209         %
0210         XXdf0dpn_tp_ijm = (XXf0_tp_ippjm - XXf0_tp_immjm)./(XXdpnm + XXdpnp);
0211         XXdf0dpn_tp_ijp = (XXf0_tp_ippjp - XXf0_tp_immjp)./(XXdpnm + XXdpnp);
0212         %
0213         %
0214         XXf0_g_ijmm = XXf0_g_ij(:,indexjm,:);
0215         XXf0_g_ijpp = XXf0_g_ij(:,indexjp,:);
0216         %
0217         XXf0_g_ijm = (1 - XXdeltam_ijm).*XXf0_g_ij + XXdeltam_ijm.*XXf0_g_ijmm;
0218         XXf0_g_ijp = (1 - XXdeltam_ijp).*XXf0_g_ippj + XXdeltam_ijp.*XXf0_g_ij;
0219         %
0220         XXf0_g_immjm = XXf0_g_ijm(indexim,:,:);
0221         XXf0_g_ippjm = XXf0_g_ijm(indexip,:,:);XXf0_g_ippjm(npn,:,:) = 0;
0222         %
0223         XXf0_g_immjp = XXf0_g_ijp(indexim,:,:);
0224         XXf0_g_ippjp = XXf0_g_ijp(indexip,:,:);XXf0_g_ippjp(npn,:,:) = 0;
0225         %
0226         XXdf0dmhu_g_ijm = (XXf0_g_ij - XXf0_g_ijmm)./XXdmhum;
0227         XXdf0dmhu_g_ijp = (XXf0_g_ijpp - XXf0_g_ij)./XXdmhup;
0228         %
0229         XXdf0dpn_g_ijm = (XXf0_g_ippjm - XXf0_g_immjm)./(XXdpnm + XXdpnp);
0230         XXdf0dpn_g_ijp = (XXf0_g_ippjp - XXf0_g_immjp)./(XXdpnm + XXdpnp);
0231         %
0232     end
0233 end
0234 %
0235 % -------------------------------------------------------------------------
0236 %
0237 %                       RF Power and fluxes calculation
0238 %
0239 % -------------------------------------------------------------------------
0240 %
0241 % Vector construction for RF fluxes
0242 %
0243 iX = (gridindex_rf.mask_m - 1)*npn + gridindex_rf.mask_p;
0244 iXX = (gridindex_rf.mask_r - 1)*nmhu*npn + iX;
0245 ixy = (gridindex_rf.mask_y - 1)*nr + gridindex_rf.mask_r;
0246 ixyn = (gridindex_rf.mask_n - 1)*ny*nr + ixy;
0247 iyn = (gridindex_rf.mask_n - 1)*ny + gridindex_rf.mask_y;
0248 %
0249 gridindex_rf = rmfield(gridindex_rf,['mask_r';'mask_m';'mask_p';'mask_n';'mask_y']);
0250 %
0251 XYP0_2piRp = xyP0_2piRp(ixy);
0252 %
0253 ZXYD_rf.pp_ipj = ZXYD_rf.pp_ipj.*XYP0_2piRp;
0254 ZXYD_rf.pp_imj = ZXYD_rf.pp_imj.*XYP0_2piRp;
0255 ZXYD_rf.pm_ipj = ZXYD_rf.pm_ipj.*XYP0_2piRp;
0256 ZXYD_rf.pm_imj = ZXYD_rf.pm_imj.*XYP0_2piRp;
0257 %
0258 XYFp_rf_imj = 0;
0259 XYFp_rf_ipj = 0;
0260 %
0261 XYf0_imj = single(XXf0_imj(iXX));
0262 XYf0_ipj = single(XXf0_ipj(iXX));
0263 %
0264 XYdf0dpn_imj = single(XXdf0dpn_imj(iXX));
0265 XYdf0dpn_ipj = single(XXdf0dpn_ipj(iXX));
0266 %
0267 XYdf0dmhu_imj = single(XXdf0dmhu_imj(iXX));
0268 XYdf0dmhu_ipj = single(XXdf0dmhu_ipj(iXX));
0269 %
0270 XY1mmhu2 = single(X1mmhu2(iX));
0271 XYpnm = single(Xpnm(iX));
0272 XYpnp = single(Xpnp(iX));
0273 %
0274 if dke_mode == 1,
0275     %
0276     ZXYD_rf_tp.pp_ipj = ZXYD_rf_tp.pp_ipj.*XYP0_2piRp;
0277     ZXYD_rf_tp.pp_imj = ZXYD_rf_tp.pp_imj.*XYP0_2piRp;
0278     ZXYD_rf_tp.pm_ipj = ZXYD_rf_tp.pm_ipj.*XYP0_2piRp;
0279     ZXYD_rf_tp.pm_imj = ZXYD_rf_tp.pm_imj.*XYP0_2piRp;
0280     %
0281     ZXYF_rf_tp.p_ipj = ZXYF_rf_tp.p_ipj.*XYP0_2piRp;
0282     ZXYF_rf_tp.p_imj = ZXYF_rf_tp.p_imj.*XYP0_2piRp;
0283     %
0284     XYf0_tp_imj = single(XXf0_tp_imj(iXX));
0285     XYf0_tp_ipj = single(XXf0_tp_ipj(iXX));
0286     %
0287     XYdf0dpn_tp_imj = single(XXdf0dpn_tp_imj(iXX));
0288     XYdf0dpn_tp_ipj = single(XXdf0dpn_tp_ipj(iXX));
0289     %
0290     XYdf0dmhu_tp_imj = single(XXdf0dmhu_tp_imj(iXX));
0291     XYdf0dmhu_tp_ipj = single(XXdf0dmhu_tp_ipj(iXX));
0292     %
0293     XYf0_g_imj = single(XXf0_g_imj(iXX));
0294     XYf0_g_ipj = single(XXf0_g_ipj(iXX));
0295     %
0296     XYdf0dpn_g_imj = single(XXdf0dpn_g_imj(iXX));
0297     XYdf0dpn_g_ipj = single(XXdf0dpn_g_ipj(iXX));
0298     %
0299     XYdf0dmhu_g_imj = single(XXdf0dmhu_g_imj(iXX));
0300     XYdf0dmhu_g_ipj = single(XXdf0dmhu_g_ipj(iXX));
0301     %
0302 end
0303 %
0304 if nargout == 1,
0305     %
0306     %f0 contribution
0307     %
0308     [XYSp0_rf_imj,XYSp0_rf_ipj] = ...
0309         fluxes_dke_jd(XYdf0dpn_imj,XYdf0dpn_ipj,XYdf0dmhu_imj,XYdf0dmhu_ipj,XYf0_imj,XYf0_ipj,...
0310                       ZXYD_rf.pp_imj,ZXYD_rf.pp_ipj,ZXYD_rf.pm_imj,ZXYD_rf.pm_ipj,XYFp_rf_imj,XYFp_rf_ipj,XY1mmhu2,XYpnm,XYpnp);
0311     %
0312     if dke_mode == 1,
0313         %
0314         [XYSp0_rf_tp_imj,XYSp0_rf_tp_ipj] = ...
0315             fluxes_dke_jd(XYdf0dpn_tp_imj,XYdf0dpn_tp_ipj,XYdf0dmhu_tp_imj,XYdf0dmhu_tp_ipj,XYf0_tp_imj,XYf0_tp_ipj,...
0316                       ZXYD_rf_tp.pp_imj,ZXYD_rf_tp.pp_ipj,ZXYD_rf_tp.pm_imj,ZXYD_rf_tp.pm_ipj,ZXYF_rf_tp.p_imj,ZXYF_rf_tp.p_ipj,XY1mmhu2,XYpnm,XYpnp);
0317         %
0318         [XYSp0_rf_g_imj,XYSp0_rf_g_ipj] = ...
0319             fluxes_dke_jd(XYdf0dpn_g_imj,XYdf0dpn_g_ipj,XYdf0dmhu_g_imj,XYdf0dmhu_g_ipj,XYf0_g_imj,XYf0_g_ipj,...
0320                       ZXYD_rf.pp_imj,ZXYD_rf.pp_ipj,ZXYD_rf.pm_imj,ZXYD_rf.pm_ipj,XYFp_rf_imj,XYFp_rf_ipj,XY1mmhu2,XYpnm,XYpnp);
0321         %
0322     end
0323 else
0324     %
0325     ZXYD_rf.mm_ijm = ZXYD_rf.mm_ijm.*XYP0_2piRp;
0326     ZXYD_rf.mm_ijp = ZXYD_rf.mm_ijp.*XYP0_2piRp;
0327     ZXYD_rf.mp_ijm = ZXYD_rf.mp_ijm.*XYP0_2piRp;
0328     ZXYD_rf.mp_ijp = ZXYD_rf.mp_ijp.*XYP0_2piRp;
0329     %
0330     XYFm_rf_ijm = 0;
0331     XYFm_rf_ijp = 0;
0332     %
0333     XYf0_ijm = single(XXf0_ijm(iXX));
0334     XYf0_ijp = single(XXf0_ijp(iXX));
0335     %
0336     XYdf0dmhu_ijm = single(XXdf0dmhu_ijm(iXX));
0337     XYdf0dmhu_ijp = single(XXdf0dmhu_ijp(iXX));
0338     %
0339     XYdf0dpn_ijm = single(XXdf0dpn_ijm(iXX));
0340     XYdf0dpn_ijp = single(XXdf0dpn_ijp(iXX));
0341     %
0342     XY1mmhu2p = single(X1mmhu2p(iX));
0343     XY1mmhu2m = single(X1mmhu2m(iX));
0344     XYpn = single(Xpn(iX));
0345     %
0346     %f0 contribution
0347     %
0348     [XYSp0_rf_imj,XYSp0_rf_ipj,XYSm0_rf_ijm,XYSm0_rf_ijp] = ...
0349         fluxes_dke_jd(XYdf0dpn_imj,XYdf0dpn_ipj,XYdf0dmhu_imj,XYdf0dmhu_ipj,XYf0_imj,XYf0_ipj,...
0350                       ZXYD_rf.pp_imj,ZXYD_rf.pp_ipj,ZXYD_rf.pm_imj,ZXYD_rf.pm_ipj,XYFp_rf_imj,XYFp_rf_ipj,XY1mmhu2,XYpnm,XYpnp,...
0351                       XYdf0dmhu_ijm,XYdf0dmhu_ijp,XYdf0dpn_ijm,XYdf0dpn_ijp,XYf0_ijm,XYf0_ijp,...
0352                       ZXYD_rf.mm_ijm,ZXYD_rf.mm_ijp,ZXYD_rf.mp_ijm,ZXYD_rf.mp_ijp,XYFm_rf_ijm,XYFm_rf_ijp,XY1mmhu2m,XY1mmhu2p,XYpn);
0353     %
0354     if dke_mode == 1,
0355         %
0356         %f1 contribution
0357         %
0358         ZXYD_rf_tp.mm_ijm = ZXYD_rf_tp.mm_ijm.*XYP0_2piRp;
0359         ZXYD_rf_tp.mm_ijp = ZXYD_rf_tp.mm_ijp.*XYP0_2piRp;
0360         ZXYD_rf_tp.mp_ijm = ZXYD_rf_tp.mp_ijm.*XYP0_2piRp;
0361         ZXYD_rf_tp.mp_ijp = ZXYD_rf_tp.mp_ijp.*XYP0_2piRp;
0362         %
0363         ZXYF_rf_tp.m_ijm = ZXYF_rf_tp.m_ijm.*XYP0_2piRp;
0364         ZXYF_rf_tp.m_ijp = ZXYF_rf_tp.m_ijp.*XYP0_2piRp;
0365         %
0366         XYf0_tp_ijm = single(XXf0_tp_ijm(iXX));
0367         XYf0_tp_ijp = single(XXf0_tp_ijp(iXX));
0368         %
0369         XYdf0dmhu_tp_ijm = single(XXdf0dmhu_tp_ijm(iXX));
0370         XYdf0dmhu_tp_ijp = single(XXdf0dmhu_tp_ijp(iXX));
0371         %
0372         XYdf0dpn_tp_ijm = single(XXdf0dpn_tp_ijm(iXX));
0373         XYdf0dpn_tp_ijp = single(XXdf0dpn_tp_ijp(iXX));
0374         %
0375         XYf0_g_ijm = single(XXf0_g_ijm(iXX));
0376         XYf0_g_ijp = single(XXf0_g_ijp(iXX));
0377         %
0378         XYdf0dmhu_g_ijm = single(XXdf0dmhu_g_ijm(iXX));
0379         XYdf0dmhu_g_ijp = single(XXdf0dmhu_g_ijp(iXX));
0380         %
0381         XYdf0dpn_g_ijm = single(XXdf0dpn_g_ijm(iXX));
0382         XYdf0dpn_g_ijp = single(XXdf0dpn_g_ijp(iXX));
0383         %
0384         [XYSp0_rf_tp_imj,XYSp0_rf_tp_ipj,XYSm0_rf_tp_ijm,XYSm0_rf_tp_ijp] = ...
0385             fluxes_dke_jd(XYdf0dpn_tp_imj,XYdf0dpn_tp_ipj,XYdf0dmhu_tp_imj,XYdf0dmhu_tp_ipj,XYf0_tp_imj,XYf0_tp_ipj,...
0386                       ZXYD_rf_tp.pp_imj,ZXYD_rf_tp.pp_ipj,ZXYD_rf_tp.pm_imj,ZXYD_rf_tp.pm_ipj,ZXYF_rf_tp.p_imj,ZXYF_rf_tp.p_ipj,XY1mmhu2,XYpnm,XYpnp,...
0387                       XYdf0dmhu_tp_ijm,XYdf0dmhu_tp_ijp,XYdf0dpn_tp_ijm,XYdf0dpn_tp_ijp,XYf0_tp_ijm,XYf0_tp_ijp,...
0388                       ZXYD_rf_tp.mm_ijm,ZXYD_rf_tp.mm_ijp,ZXYD_rf_tp.mp_ijm,ZXYD_rf_tp.mp_ijp,ZXYF_rf_tp.m_ijm,ZXYF_rf_tp.m_ijp,XY1mmhu2m,XY1mmhu2p,XYpn);
0389         %
0390         [XYSp0_rf_g_imj,XYSp0_rf_g_ipj,XYSm0_rf_g_ijm,XYSm0_rf_g_ijp] = ...
0391             fluxes_dke_jd(XYdf0dpn_g_imj,XYdf0dpn_g_ipj,XYdf0dmhu_g_imj,XYdf0dmhu_g_ipj,XYf0_g_imj,XYf0_g_ipj,...
0392                       ZXYD_rf.pp_imj,ZXYD_rf.pp_ipj,ZXYD_rf.pm_imj,ZXYD_rf.pm_ipj,XYFp_rf_imj,XYFp_rf_ipj,XY1mmhu2,XYpnm,XYpnp,...
0393                       XYdf0dmhu_g_ijm,XYdf0dmhu_g_ijp,XYdf0dpn_g_ijm,XYdf0dpn_g_ijp,XYf0_g_ijm,XYf0_g_ijp,...
0394                       ZXYD_rf.mm_ijm,ZXYD_rf.mm_ijp,ZXYD_rf.mp_ijm,ZXYD_rf.mp_ijp,XYFm_rf_ijm,XYFm_rf_ijp,XY1mmhu2m,XY1mmhu2p,XYpn);
0395         %
0396         %  Wave-induced radial transport by adjoint technique (P. Helander model)
0397         %
0398         ZXYD_rf_tp.mm1_ijm = ZXYD_rf_tp.mm1_ijm.*XYP0_2piRp;
0399         ZXYD_rf_tp.mm1_ijp = ZXYD_rf_tp.mm1_ijp.*XYP0_2piRp;
0400         ZXYD_rf_tp.mp1_ijm = ZXYD_rf_tp.mp1_ijm.*XYP0_2piRp;
0401         ZXYD_rf_tp.mp1_ijp = ZXYD_rf_tp.mp1_ijp.*XYP0_2piRp;
0402         %
0403         %
0404         [XYSp0_rf_tp1_imj,XYSp0_rf_tp1_ipj,XYSm0_rf_tp1_ijm,XYSm0_rf_tp1_ijp] = ...
0405             fluxes_dke_jd(XYdf0dpn_imj,XYdf0dpn_ipj,XYdf0dmhu_imj,XYdf0dmhu_ipj,XYf0_imj,XYf0_ipj,...
0406                       ZXYD_rf_tp.pp_imj,ZXYD_rf_tp.pp_ipj,ZXYD_rf_tp.pm_imj,ZXYD_rf_tp.pm_ipj,XYFp_rf_imj,XYFp_rf_ipj,XY1mmhu2,XYpnm,XYpnp,...
0407                       XYdf0dmhu_ijm,XYdf0dmhu_ijp,XYdf0dpn_ijm,XYdf0dpn_ijp,XYf0_ijm,XYf0_ijp,...
0408                       ZXYD_rf_tp.mm1_ijm,ZXYD_rf_tp.mm1_ijp,ZXYD_rf_tp.mp1_ijm,ZXYD_rf_tp.mp1_ijp,XYFm_rf_ijm,XYFm_rf_ijp,XY1mmhu2m,XY1mmhu2p,XYpn);
0409     
0410         %
0411     end
0412     %
0413 end
0414 %
0415 clear XYP0_2piRp ZXYD_rf ZXYD_rf_tp ZXYF_rf_tp 
0416 %
0417 clear XYf0_imj XYf0_ipj XYdf0dpn_imj XYdf0dpn_ipj XYdf0dmhu_imj XYdf0dmhu_ipj
0418 clear XYf0_tp_imj XYf0_tp_ipj XYdf0dpn_tp_imj XYdf0dpn_tp_ipj XYdf0dmhu_tp_imj XYdf0dmhu_tp_ipj
0419 clear XYf0_g_imj XYf0_g_ipj XYdf0dpn_g_imj XYdf0dpn_g_ipj XYdf0dmhu_g_imj XYdf0dmhu_g_ipj
0420 %
0421 clear XYf0_ijm XYf0_ijp XYdf0dmhu_ijm XYdf0dmhu_ijp XYdf0dpn_ijm XYdf0dpn_ijp
0422 clear XYf0_tp_ijm XYf0_tp_ijp XYdf0dmhu_tp_ijm XYdf0dmhu_tp_ijp XYdf0dpn_tp_ijm XYdf0dpn_tp_ijp
0423 clear XYf0_g_ijm XYf0_g_ijp XYdf0dmhu_g_ijm XYdf0dmhu_g_ijp XYdf0dpn_g_ijm XYdf0dpn_g_ijp
0424 %
0425 clear XY1mmhu2 XY1mmhu2p 
0426 %
0427 XYdmhu    = single(Xdmhu(iX));     
0428 XYgammam  = single(Xgammam(iX));
0429 XYgammap  = single(Xgammap(iX));
0430 XYdpnm    = single(Xdpnm(iX));     
0431 XYdpn     = single(Xdpn(iX));     
0432 XYlambda  = single(XXlambda(iXX)); 
0433 %
0434 xqratio = xqtilde./xqhat;
0435 xynqratio = repmat(xqratio.',[1,ny,nn_rf]);
0436 %
0437 XYIntegrandm = XYpnm.*XYpnm.*XYpnm.*XYlambda.*XYSp0_rf_imj./XYgammam;
0438 XYIntegrandp = XYpnp.*XYpnp.*XYpnp.*XYlambda.*XYSp0_rf_ipj./XYgammap;
0439 %
0440 XYIntegrand = XYdmhu.*XYdpn.*(XYIntegrandm + XYIntegrandp)/2;
0441 %
0442 XYIntegrand = sparse(double(iX),double(ixyn),double(XYIntegrand),npn*nmhu,nr*ny*nn_rf);
0443 %
0444 ZP0.xyn_rf_fsav = 2*pi*xynqratio.*reshape(full(sum(XYIntegrand)),[nr,ny,nn_rf]);
0445 %
0446 ZP0.xy_rf_fsav = sum(ZP0.xyn_rf_fsav,3);%sum over harmonic numbers
0447 ZP0.x_rf_fsav = sum(ZP0.xy_rf_fsav,2);%sum over rays
0448 %
0449 if dke_mode == 1,
0450     %
0451     XYIntegrandm = XYpnm.*XYpnm.*XYpnm.*XYlambda.*XYSp0_rf_tp_imj./XYgammam;
0452     XYIntegrandp = XYpnp.*XYpnp.*XYpnp.*XYlambda.*XYSp0_rf_tp_ipj./XYgammap;
0453     %
0454     XYIntegrand = XYdmhu.*XYdpn.*(XYIntegrandm + XYIntegrandp)/2;
0455     %
0456     XYIntegrand = sparse(double(iX),double(ixyn),double(XYIntegrand),npn*nmhu,nr*ny*nn_rf);
0457     %
0458     ZP0.xyn_rf_tp_fsav = 2*pi*xynqratio.*reshape(full(sum(XYIntegrand)),[nr,ny,nn_rf]);
0459     %
0460     XYIntegrandm = XYpnm.*XYpnm.*XYpnm.*XYlambda.*XYSp0_rf_g_imj./XYgammam;
0461     XYIntegrandp = XYpnp.*XYpnp.*XYpnp.*XYlambda.*XYSp0_rf_g_ipj./XYgammap;
0462     %
0463     XYIntegrand = XYdmhu.*XYdpn.*(XYIntegrandm + XYIntegrandp)/2;
0464     %
0465     XYIntegrand = sparse(double(iX),double(ixyn),double(XYIntegrand),npn*nmhu,nr*ny*nn_rf);
0466     %
0467     ZP0.xyn_rf_g_fsav = 2*pi*xynqratio.*reshape(full(sum(XYIntegrand)),[nr,ny,nn_rf]);
0468     %
0469     ZP0.xy_rf_tp_fsav = sum(ZP0.xyn_rf_tp_fsav,3);%sum over harmonic numbers
0470     ZP0.xy_rf_g_fsav = sum(ZP0.xyn_rf_g_fsav,3);%sum over harmonic numbers
0471     %
0472     ZP0.x_rf_tp_fsav = sum(ZP0.xy_rf_tp_fsav,2);%sum over rays
0473     ZP0.x_rf_g_fsav = sum(ZP0.xy_rf_g_fsav,2);%sum over rays
0474     %
0475     if dke_mode == 1 & nargout > 1,
0476         %
0477         %  Wave-induced radial transport by adjoint technique (P. Helander model)
0478         %
0479         XYmhu = single(Xmhu(iX));     
0480         XYRlambda_m = single(Xdpn(iXX));
0481         %
0482         XYIntegrand = XYdmhu.*XYdpnm.*XYpnm.*XYpnm.*XYmhu.*XYlambda.*XYSp0_rf_tp1_imj;
0483         XYIntegrand = sparse(double(iX),double(ixyn),double(XYIntegrand),npn*nmhu,nr*ny*nn_rf);
0484         xggradpsi21 = 2*pi*xqratio.*(xR0.*xBT0./xB0).*sum(sum(reshape(full(sum(XYIntegrand)),[nr,ny,nn_rf]),3),2).';
0485         %
0486         XYIntegrand = XYdmhu.*XYdpn.*XYpn.*XYpn.*XY1mmhu2m.*XYlambda.*XYRlambda_m.*XYSm0_rf_tp1_ijm;
0487         XYIntegrand = sparse(double(iX),double(ixyn),double(XYIntegrand),npn*nmhu,nr*ny*nn_rf);
0488         xggradpsi1 = 2*pi*xynqratio.*reshape(full(sum(XYIntegrand)),[nr,ny,nn_rf]);
0489     end
0490     %
0491 end
0492 %
0493 clear iX ixy ixyn
0494 clear XYIntegrandm XYIntegrand XYIntegrandp
0495 clear XYdmhu XYgamma XYgammap XYdpnm XYpnm XYpnp XYlambda XYpn XY1mmhu2m XYmhu XYdpn XYRlambda_m
0496 %
0497 if nargout == 1,
0498     %
0499     return
0500     %
0501 end
0502 %
0503 % Sum over rays and harmonics for RF fluxes
0504 %
0505 ZXXS.p0_rf_imj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_imj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0506 ZXXS.p0_rf_ipj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_ipj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0507 ZXXS.m0_rf_ijm = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_ijm),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0508 ZXXS.m0_rf_ijp = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_ijp),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0509 %
0510 if dke_mode == 1,
0511     %
0512     ZXXS.p0_rf_tp_imj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_tp_imj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0513     ZXXS.p0_rf_tp_ipj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_tp_ipj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0514     ZXXS.m0_rf_tp_ijm = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_tp_ijm),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0515     ZXXS.m0_rf_tp_ijp = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_tp_ijp),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0516     %
0517     ZXXS.p0_rf_g_imj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_g_imj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0518     ZXXS.p0_rf_g_ipj = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSp0_rf_g_ipj),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0519     ZXXS.m0_rf_g_ijm = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_g_ijm),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0520     ZXXS.m0_rf_g_ijp = reshape(full(sum(sparse(double(iyn),double(iXX),double(XYSm0_rf_g_ijp),ny*nn_rf,npn*nmhu*nr),1)),[npn,nmhu,nr]);
0521     %
0522 end
0523 %
0524 clear iXX iyn
0525 clear XYSp0_rf_imj XYSp0_rf_ipj XYSm0_rf_ijm XYSm0_rf_ijp
0526 clear XYSp0_rf_tp_imj XYSp0_rf_tp_ipj XYSm0_rf_tp_ijm XYSm0_rf_tp_ijp
0527 clear XYSp0_rf_tp1_imj XYSp0_rf_tp1_ipj XYSm0_rf_tp1_ijm XYSm0_rf_tp1_ijp
0528 clear XYSp0_rf_g_imj XYSp0_rf_g_ipj XYSm0_rf_g_ijm XYSm0_rf_g_ijp
0529 %
0530 % -------------------------------------------------------------------------
0531 %
0532 %                       Other Power and Fluxes Calculations
0533 %
0534 % -------------------------------------------------------------------------
0535 %
0536 %   Collisional fluxes and power
0537 %
0538 XXIntegrand_im = XXpnm(:,:,rdke).*XXpnm(:,:,rdke).*XXpnm(:,:,rdke).*XXlambda(:,:,rdke)./XXgammam(:,:,rdke);
0539 XXIntegrand_ip = XXpnp(:,:,rdke).*XXpnp(:,:,rdke).*XXpnp(:,:,rdke).*XXlambda(:,:,rdke)./XXgammap(:,:,rdke);
0540 %
0541 % f0 contribution
0542 %
0543 [ZXXS.p0_c_imj,ZXXS.p0_c_ipj,ZXXS.m0_c_ijm,ZXXS.m0_c_ijp] = ...
0544     fluxes_dke_jd(XXdf0dpn_imj(:,:,rdke),XXdf0dpn_ipj(:,:,rdke),XXdf0dmhu_imj(:,:,rdke),XXdf0dmhu_ipj(:,:,rdke),XXf0_imj(:,:,rdke),XXf0_ipj(:,:,rdke),...
0545                   ZXXD_c.pp_imj(:,:,rdke),ZXXD_c.pp_ipj(:,:,rdke),ZXXD_c.pm_imj(:,:,rdke),ZXXD_c.pm_ipj(:,:,rdke),ZXXF_c.p_imj(:,:,rdke),ZXXF_c.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0546                   XXdf0dmhu_ijm(:,:,rdke),XXdf0dmhu_ijp(:,:,rdke),XXdf0dpn_ijm(:,:,rdke),XXdf0dpn_ijp(:,:,rdke),XXf0_ijm(:,:,rdke),XXf0_ijp(:,:,rdke),...
0547                   ZXXD_c.mm_ijm(:,:,rdke),ZXXD_c.mm_ijp(:,:,rdke),ZXXD_c.mp_ijm(:,:,rdke),ZXXD_c.mp_ijp(:,:,rdke),ZXXF_c.m_ijm(:,:,rdke),ZXXF_c.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0548 %
0549 ZP0.x_c_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_c_imj + XXIntegrand_ip.*ZXXS.p0_c_ipj)/2,1),1).';
0550 %
0551 if dke_mode >= 1
0552     %
0553     %f0_tp contribution
0554     %
0555     [ZXXS.p0_c_tp_imj,ZXXS.p0_c_tp_ipj,ZXXS.m0_c_tp_ijm,ZXXS.m0_c_tp_ijp] = ...
0556         fluxes_dke_jd(XXdf0dpn_tp_imj(:,:,rdke),XXdf0dpn_tp_ipj(:,:,rdke),XXdf0dmhu_tp_imj(:,:,rdke),XXdf0dmhu_tp_ipj(:,:,rdke),XXf0_tp_imj(:,:,rdke),XXf0_tp_ipj(:,:,rdke),...
0557                       ZXXD_c_tp.pp_imj(:,:,rdke),ZXXD_c_tp.pp_ipj(:,:,rdke),ZXXD_c_tp.pm_imj(:,:,rdke),ZXXD_c_tp.pm_ipj(:,:,rdke),ZXXF_c_tp.p_imj(:,:,rdke),ZXXF_c_tp.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0558                       XXdf0dmhu_tp_ijm(:,:,rdke),XXdf0dmhu_tp_ijp(:,:,rdke),XXdf0dpn_tp_ijm(:,:,rdke),XXdf0dpn_tp_ijp(:,:,rdke),XXf0_tp_ijm(:,:,rdke),XXf0_tp_ijp(:,:,rdke),...
0559                       ZXXD_c_tp.mm_ijm(:,:,rdke),ZXXD_c_tp.mm_ijp(:,:,rdke),ZXXD_c_tp.mp_ijm(:,:,rdke),ZXXD_c_tp.mp_ijp(:,:,rdke),ZXXF_c_tp.m_ijm(:,:,rdke),ZXXF_c_tp.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0560     %
0561     ZP0.x_c_tp_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_c_tp_imj + XXIntegrand_ip.*ZXXS.p0_c_tp_ipj)/2,1),1).';
0562     %
0563     %f0_g contribution
0564     %
0565     [ZXXS.p0_c_g_imj,ZXXS.p0_c_g_ipj,ZXXS.m0_c_g_ijm,ZXXS.m0_c_g_ijp] = ...
0566         fluxes_dke_jd(XXdf0dpn_g_imj(:,:,rdke),XXdf0dpn_g_ipj(:,:,rdke),XXdf0dmhu_g_imj(:,:,rdke),XXdf0dmhu_g_ipj(:,:,rdke),XXf0_g_imj(:,:,rdke),XXf0_g_ipj(:,:,rdke),...
0567                       ZXXD_c.pp_imj(:,:,rdke),ZXXD_c.pp_ipj(:,:,rdke),ZXXD_c.pm_imj(:,:,rdke),ZXXD_c.pm_ipj(:,:,rdke),ZXXF_c.p_imj(:,:,rdke),ZXXF_c.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0568                       XXdf0dmhu_g_ijm(:,:,rdke),XXdf0dmhu_g_ijp(:,:,rdke),XXdf0dpn_g_ijm(:,:,rdke),XXdf0dpn_g_ijp(:,:,rdke),XXf0_g_ijm(:,:,rdke),XXf0_g_ijp(:,:,rdke),...
0569                       ZXXD_c.mm_ijm(:,:,rdke),ZXXD_c.mm_ijp(:,:,rdke),ZXXD_c.mp_ijm(:,:,rdke),ZXXD_c.mp_ijp(:,:,rdke),ZXXF_c.m_ijm(:,:,rdke),ZXXF_c.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0570     %
0571     ZP0.x_c_g_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_c_g_imj + XXIntegrand_ip.*ZXXS.p0_c_g_ipj)/2,1),1).';
0572 end
0573 %
0574 %   Ohmic fluxes and power
0575 %
0576 %
0577 % f0 contribution
0578 %
0579 [ZXXS.p0_e_imj,ZXXS.p0_e_ipj,ZXXS.m0_e_ijm,ZXXS.m0_e_ijp] = ...
0580     fluxes_dke_jd(XXdf0dpn_imj(:,:,rdke),XXdf0dpn_ipj(:,:,rdke),XXdf0dmhu_imj(:,:,rdke),XXdf0dmhu_ipj(:,:,rdke),XXf0_imj(:,:,rdke),XXf0_ipj(:,:,rdke),...
0581                   ZXXD_e.pp_imj(:,:,rdke),ZXXD_e.pp_ipj(:,:,rdke),ZXXD_e.pm_imj(:,:,rdke),ZXXD_e.pm_ipj(:,:,rdke),ZXXF_e.p_imj(:,:,rdke),ZXXF_e.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0582                   XXdf0dmhu_ijm(:,:,rdke),XXdf0dmhu_ijp(:,:,rdke),XXdf0dpn_ijm(:,:,rdke),XXdf0dpn_ijp(:,:,rdke),XXf0_ijm(:,:,rdke),XXf0_ijp(:,:,rdke),...
0583                   ZXXD_e.mm_ijm(:,:,rdke),ZXXD_e.mm_ijp(:,:,rdke),ZXXD_e.mp_ijm(:,:,rdke),ZXXD_e.mp_ijp(:,:,rdke),ZXXF_e.m_ijm(:,:,rdke),ZXXF_e.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0584 %
0585 ZP0.x_e_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_e_imj + XXIntegrand_ip.*ZXXS.p0_e_ipj)/2,1),1).';
0586 %
0587 if dke_mode >= 1
0588     %
0589     %f0_tp contribution
0590     %
0591     [ZXXS.p0_e_tp_imj,ZXXS.p0_e_tp_ipj,ZXXS.m0_e_tp_ijm,ZXXS.m0_e_tp_ijp] = ...
0592         fluxes_dke_jd(XXdf0dpn_tp_imj(:,:,rdke),XXdf0dpn_tp_ipj(:,:,rdke),XXdf0dmhu_tp_imj(:,:,rdke),XXdf0dmhu_tp_ipj(:,:,rdke),XXf0_tp_imj(:,:,rdke),XXf0_tp_ipj(:,:,rdke),...
0593                       ZXXD_e_tp.pp_imj(:,:,rdke),ZXXD_e_tp.pp_ipj(:,:,rdke),ZXXD_e_tp.pm_imj(:,:,rdke),ZXXD_e_tp.pm_ipj(:,:,rdke),ZXXF_e_tp.p_imj(:,:,rdke),ZXXF_e_tp.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0594                       XXdf0dmhu_tp_ijm(:,:,rdke),XXdf0dmhu_tp_ijp(:,:,rdke),XXdf0dpn_tp_ijm(:,:,rdke),XXdf0dpn_tp_ijp(:,:,rdke),XXf0_tp_ijm(:,:,rdke),XXf0_tp_ijp(:,:,rdke),...
0595                       ZXXD_e_tp.mm_ijm(:,:,rdke),ZXXD_e_tp.mm_ijp(:,:,rdke),ZXXD_e_tp.mp_ijm(:,:,rdke),ZXXD_e_tp.mp_ijp(:,:,rdke),ZXXF_e_tp.m_ijm(:,:,rdke),ZXXF_e_tp.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0596     %
0597     ZP0.x_e_tp_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_e_tp_imj + XXIntegrand_ip.*ZXXS.p0_e_tp_ipj)/2,1),1).';
0598     %
0599     %f0_g contribution
0600     %
0601     [ZXXS.p0_e_g_imj,ZXXS.p0_e_g_ipj,ZXXS.m0_e_g_ijm,ZXXS.m0_e_g_ijp] = ...
0602         fluxes_dke_jd(XXdf0dpn_g_imj(:,:,rdke),XXdf0dpn_g_ipj(:,:,rdke),XXdf0dmhu_g_imj(:,:,rdke),XXdf0dmhu_g_ipj(:,:,rdke),XXf0_g_imj(:,:,rdke),XXf0_g_ipj(:,:,rdke),...
0603                       ZXXD_e.pp_imj(:,:,rdke),ZXXD_e.pp_ipj(:,:,rdke),ZXXD_e.pm_imj(:,:,rdke),ZXXD_e.pm_ipj(:,:,rdke),ZXXF_e.p_imj(:,:,rdke),ZXXF_e.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0604                       XXdf0dmhu_g_ijm(:,:,rdke),XXdf0dmhu_g_ijp(:,:,rdke),XXdf0dpn_g_ijm(:,:,rdke),XXdf0dpn_g_ijp(:,:,rdke),XXf0_g_ijm(:,:,rdke),XXf0_g_ijp(:,:,rdke),...
0605                       ZXXD_e.mm_ijm(:,:,rdke),ZXXD_e.mm_ijp(:,:,rdke),ZXXD_e.mp_ijm(:,:,rdke),ZXXD_e.mp_ijp(:,:,rdke),ZXXF_e.m_ijm(:,:,rdke),ZXXF_e.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0606     %
0607     ZP0.x_e_g_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_e_g_imj + XXIntegrand_ip.*ZXXS.p0_e_g_ipj)/2,1),1).';
0608 end
0609 %
0610 %   Synchrotron fluxes and power
0611 %
0612 %
0613 % f0 contribution
0614 %
0615 [ZXXS.p0_s_imj,ZXXS.p0_s_ipj,ZXXS.m0_s_ijm,ZXXS.m0_s_ijp] = ...
0616     fluxes_dke_jd(XXdf0dpn_imj(:,:,rdke),XXdf0dpn_ipj(:,:,rdke),XXdf0dmhu_imj(:,:,rdke),XXdf0dmhu_ipj(:,:,rdke),XXf0_imj(:,:,rdke),XXf0_ipj(:,:,rdke),...
0617                   ZXXD_s.pp_imj(:,:,rdke),ZXXD_s.pp_ipj(:,:,rdke),ZXXD_s.pm_imj(:,:,rdke),ZXXD_s.pm_ipj(:,:,rdke),ZXXF_s.p_imj(:,:,rdke),ZXXF_s.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0618                   XXdf0dmhu_ijm(:,:,rdke),XXdf0dmhu_ijp(:,:,rdke),XXdf0dpn_ijm(:,:,rdke),XXdf0dpn_ijp(:,:,rdke),XXf0_ijm(:,:,rdke),XXf0_ijp(:,:,rdke),...
0619                   ZXXD_s.mm_ijm(:,:,rdke),ZXXD_s.mm_ijp(:,:,rdke),ZXXD_s.mp_ijm(:,:,rdke),ZXXD_s.mp_ijp(:,:,rdke),ZXXF_s.m_ijm(:,:,rdke),ZXXF_s.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0620 %
0621 ZP0.x_s_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_s_imj + XXIntegrand_ip.*ZXXS.p0_s_ipj)/2,1),1).';
0622 %
0623 if dke_mode >= 1
0624     %
0625     %f0_tp contribution
0626     %
0627     [ZXXS.p0_s_tp_imj,ZXXS.p0_s_tp_ipj,ZXXS.m0_s_tp_ijm,ZXXS.m0_s_tp_ijp] = ...
0628         fluxes_dke_jd(XXdf0dpn_tp_imj(:,:,rdke),XXdf0dpn_tp_ipj(:,:,rdke),XXdf0dmhu_tp_imj(:,:,rdke),XXdf0dmhu_tp_ipj(:,:,rdke),XXf0_tp_imj(:,:,rdke),XXf0_tp_ipj(:,:,rdke),...
0629                       ZXXD_s_tp.pp_imj(:,:,rdke),ZXXD_s_tp.pp_ipj(:,:,rdke),ZXXD_s_tp.pm_imj(:,:,rdke),ZXXD_s_tp.pm_ipj(:,:,rdke),ZXXF_s_tp.p_imj(:,:,rdke),ZXXF_s_tp.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0630                       XXdf0dmhu_tp_ijm(:,:,rdke),XXdf0dmhu_tp_ijp(:,:,rdke),XXdf0dpn_tp_ijm(:,:,rdke),XXdf0dpn_tp_ijp(:,:,rdke),XXf0_tp_ijm(:,:,rdke),XXf0_tp_ijp(:,:,rdke),...
0631                       ZXXD_s_tp.mm_ijm(:,:,rdke),ZXXD_s_tp.mm_ijp(:,:,rdke),ZXXD_s_tp.mp_ijm(:,:,rdke),ZXXD_s_tp.mp_ijp(:,:,rdke),ZXXF_s_tp.m_ijm(:,:,rdke),ZXXF_s_tp.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0632     %
0633     ZP0.x_s_tp_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_s_tp_imj + XXIntegrand_ip.*ZXXS.p0_s_tp_ipj)/2,1),1).';
0634     %
0635     %f0_g contribution
0636     %
0637     [ZXXS.p0_s_g_imj,ZXXS.p0_s_g_ipj,ZXXS.m0_s_g_ijm,ZXXS.m0_s_g_ijp] = ...
0638         fluxes_dke_jd(XXdf0dpn_g_imj(:,:,rdke),XXdf0dpn_g_ipj(:,:,rdke),XXdf0dmhu_g_imj(:,:,rdke),XXdf0dmhu_g_ipj(:,:,rdke),XXf0_g_imj(:,:,rdke),XXf0_g_ipj(:,:,rdke),...
0639                       ZXXD_s.pp_imj(:,:,rdke),ZXXD_s.pp_ipj(:,:,rdke),ZXXD_s.pm_imj(:,:,rdke),ZXXD_s.pm_ipj(:,:,rdke),ZXXF_s.p_imj(:,:,rdke),ZXXF_s.p_ipj(:,:,rdke),XX1mmhu2(:,:,rdke),XXpnm(:,:,rdke),XXpnp(:,:,rdke),...
0640                       XXdf0dmhu_g_ijm(:,:,rdke),XXdf0dmhu_g_ijp(:,:,rdke),XXdf0dpn_g_ijm(:,:,rdke),XXdf0dpn_g_ijp(:,:,rdke),XXf0_g_ijm(:,:,rdke),XXf0_g_ijp(:,:,rdke),...
0641                       ZXXD_s.mm_ijm(:,:,rdke),ZXXD_s.mm_ijp(:,:,rdke),ZXXD_s.mp_ijm(:,:,rdke),ZXXD_s.mp_ijp(:,:,rdke),ZXXF_s.m_ijm(:,:,rdke),ZXXF_s.m_ijp(:,:,rdke),XX1mmhu2m(:,:,rdke),XX1mmhu2p(:,:,rdke),XXpn(:,:,rdke));
0642     %
0643     ZP0.x_s_g_fsav = 2*pi*xqratio(rdke).*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_s_g_imj + XXIntegrand_ip.*ZXXS.p0_s_g_ipj)/2,1),1).';
0644 end
0645 %
0646 XXSp0_imj = ZXXS.p0_c_imj + ZXXS.p0_e_imj + ZXXS.p0_s_imj + ZXXS.p0_rf_imj;
0647 XXSm0_ijm = ZXXS.m0_c_ijm + ZXXS.m0_e_ijm + ZXXS.m0_s_ijm + ZXXS.m0_rf_ijm;
0648 XXSp0_ipj = ZXXS.p0_c_ipj + ZXXS.p0_e_ipj + ZXXS.p0_s_ipj + ZXXS.p0_rf_ipj;
0649 XXSm0_ijp = ZXXS.m0_c_ijp + ZXXS.m0_e_ijp + ZXXS.m0_s_ijp + ZXXS.m0_rf_ijp;
0650 %
0651 % -------------------------------------------------------------------------
0652 %
0653 %                       Other Moments Calculation
0654 %
0655 % -------------------------------------------------------------------------
0656 %
0657 bounce_mode = dkeparam.bounce_mode;
0658 temppath = dkepath.temppath;
0659 %
0660 display_mode = display.display_mode;
0661 %
0662 nr = gridDKE.nr;
0663 npn = gridDKE.npn;
0664 nmhu = gridDKE.nmhu;
0665 nr_dke = gridDKE.nr_dke;
0666 rdke = gridDKE.rdke;
0667 mhu = gridDKE.mhu;
0668 Xdmhumm = gridDKE.Xdmhumm;
0669 Xdmhupp = gridDKE.Xdmhupp;
0670 Xdpnmm = gridDKE.Xdpnmm;
0671 Xdpn = gridDKE.Xdpn;
0672 Xdpnpp = gridDKE.Xdpnpp;
0673 Xmhumm = gridDKE.Xmhumm;
0674 Xmhum = gridDKE.Xmhum;
0675 Xmhu = gridDKE.Xmhu;
0676 Xmhup = gridDKE.Xmhup;
0677 Xmhupp = gridDKE.Xmhupp;
0678 Xpnmm = gridDKE.Xpnmm;
0679 Xpnpp = gridDKE.Xpnpp;
0680 Xpn2m = gridDKE.Xpn2m;
0681 Xpn2 = gridDKE.Xpn2;
0682 Xpn2p = gridDKE.Xpn2p;
0683 xpsin = gridDKE.xpsin;
0684 XXmhu = gridDKE.XXmhu;
0685 XXpn2 = gridDKE.XXpn2;
0686 XXpn = gridDKE.XXpn;
0687 XXpn2m = gridDKE.XXpn2m;
0688 %
0689 pn = gridDKE.pn;
0690 pn2 = gridDKE.pn2;
0691 pn3 = gridDKE.pn3;
0692 xpn = gridDKE.xpn;
0693 xpn2 = gridDKE.xpn2;
0694 xpn3 = gridDKE.xpn3;dpn = gridDKE.dpn;
0695 mhu = gridDKE.mhu;
0696 dmhu = gridDKE.dmhu;
0697 pn2 = gridDKE.pn2;
0698 masku = gridDKE.masku;
0699 maskl = gridDKE.maskl;
0700 ir_display = gridDKE.ir_display_out;
0701 %
0702 xpsin = gridDKE.xpsin;
0703 xrho = gridDKE.xrho;
0704 %
0705 z = Zmomcoef.z;
0706 v = Zmomcoef.v;
0707 z2 = Zmomcoef.z2;
0708 xz = Zmomcoef.xz;
0709 xz2 = Zmomcoef.xz2;
0710 gamma = Zmomcoef.gamma;
0711 xgamma = Zmomcoef.xgamma;
0712 sigma = Zmomcoef.sigma;
0713 xsigma = Zmomcoef.xsigma;
0714 J1 = Zmomcoef.J1;
0715 J2 = Zmomcoef.J2;
0716 J3 = Zmomcoef.J3;
0717 J4 = Zmomcoef.J4;
0718 xJ1 = Zmomcoef.xJ1;
0719 xJ2 = Zmomcoef.xJ2;
0720 xJ3 = Zmomcoef.xJ3;
0721 xJ4 = Zmomcoef.xJ4;
0722 XXgamma = Zmomcoef.XXgamma;
0723 %
0724 XXdeltam_ijp = Zdelta.XXdeltam_ijp;
0725 XXdeltam_ijm = Zdelta.XXdeltam_ijm;   
0726 XXdeltap_imjmm = Zdelta.XXdeltap_imjmm;
0727 XXdeltap_imjpp = Zdelta.XXdeltap_imjpp;
0728 XXdeltap_ipjmm = Zdelta.XXdeltap_ipjmm;
0729 XXdeltap_ipjpp = Zdelta.XXdeltap_ipjpp;
0730 XXdeltam_immjm = Zdelta.XXdeltam_immjm;
0731 XXdeltam_ippjm = Zdelta.XXdeltam_ippjm;
0732 XXdeltam_immjp = Zdelta.XXdeltam_immjp;
0733 XXdeltam_ippjp = Zdelta.XXdeltam_ippjp;
0734 %
0735 XXlossripple = Zmripple.XXlossripple;
0736 %
0737 XXlambda_p2m2p2 = Zbouncecoef.XXlambda_p2m2p2;
0738 XXlambda_b_p1m1 = Zbouncecoef.XXlambda_b_p1m1;
0739 XXlambda_b_p2 = Zbouncecoef.XXlambda_b_p2;
0740 XXheaviside = Zbouncecoef.XXheaviside;
0741 Zmask_f0_t = Zbouncecoef.Zmask_f0_t;
0742 ns_f = Zbouncecoef.ns_f; 
0743 xq = Zbouncecoef.xq;
0744 xqbar =  Zbouncecoef.xqbar;
0745 xqcronos =  Zbouncecoef.xqcronos;
0746 Zmask_g = Zbouncecoef.Zmask_g;
0747 ns_g = Zbouncecoef.ns_g;
0748 xftp_norm = mksa.xftp_norm_ref;
0749 %
0750 xB0 = equilDKE.xB0;        
0751 Rp = equilDKE.Rp;
0752 xx0 = equilDKE.xx0;
0753 xBT0 = equilDKE.xBT0;
0754 xZeff = equilDKE.xZeff;
0755 %
0756 betath_ref = mksa.betath_ref;
0757 %
0758 [xTe_norm,xne_norm,xzni_norm,xzTi_norm,xnloss_norm,xbetath,xlnc_e,xnhu,xrnhuth] = normcoef_dke_yp(mksa,equilDKE,gridDKE);   
0759 %
0760 %   Density
0761 %
0762 Znorm.x_0 = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0763 Znorm.x_tp = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0764 Znorm.x_g = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_g(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0765 %
0766 if dkeparam.coll_mode == 3,%Lorentz model (test mode only, see documentation)
0767     Znorm.x_L_tp = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_tp(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0768     Znorm.x_L_g = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_g(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0769 end 
0770 %
0771 Znorm.x_0_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,rdke).*XXf0(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0772 Znorm.x_tp_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda_b_p1m1(:,:,rdke).*XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0773 Znorm.x_g_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,rdke).*XXf0_g(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0774 %
0775 if dkeparam.coll_mode == 3,%Lorentz model (test mode only, see documentation)
0776     Znorm.x_L_tp_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda_b_p1m1(:,:,rdke).*XXf0_L_tp(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0777     Znorm.x_L_g_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,rdke).*XXf0_L_g(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0778 end 
0779 %
0780 % When radial transport is turned off, the distribution function is renormalized before calculating all other
0781 % moments. This makes physical sense since no electron is meant to be lost
0782 % or gained (! still need to include run-away contribution - Yves?)
0783 %
0784 if radial_mode ~= 1,
0785     xnorm_0_fsav = 2*pi.*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXf0.*XXpn2,1),1)';%Same as Znorm.x_0_fsav but at all radius (just for vector calculations in dke_mode=1)
0786     XXf0 = XXf0.*shiftdim(repmat(xne_norm'./xnorm_0_fsav',[1,npn,nmhu]),1);
0787 end
0788 %
0789 %
0790 xc_g = -(Znorm.x_tp_fsav + Znorm.x_g_fsav)./xne_norm(rdke);%Flux surface averaged density conservation
0791 XXxc_g = shiftdim(repmat(repmat(xc_g(:),[1,npn,1]),[1,1,nmhu]),1);
0792 XXf0_g(:,:,rdke) = XXf0_g(:,:,rdke) + XXxc_g.*XXf0(:,:,rdke);%Normalisation in order to maintain the flux surface averaged density at the right level
0793 %
0794 %   Current
0795 %
0796 Zcurr.x_0 = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0797 Zcurr.x_tp = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0798 Zcurr.x_g = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_g(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0799 Zcurr.x_tot = Zcurr.x_0 + Zcurr.x_tp + Zcurr.x_g;%Total current density
0800 %
0801 if dkeparam.coll_mode == 3,%Lorentz model (test mode only, see documentation)
0802     Zcurr.x_L_tp = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0803     Zcurr.x_L_g = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_g(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0804     Zcurr.x_L_tot = Zcurr.x_0 + Zcurr.x_L_tp + Zcurr.x_L_g;%Total current density
0805 end 
0806 %
0807 %   Flux surface averaged current (two definitions, fsav is LUKE, vfsav is CRONOS)
0808 %
0809 Zcurr.x_0_fsav = (xq(rdke)./xqbar(rdke)).*Zcurr.x_0;
0810 Zcurr.x_0_vfsav = (xqtilde(rdke)./xqhat(rdke)).*Zcurr.x_0;
0811 Zcurr.x_0_vcfsav = (xB0(rdke).*xqcronos(rdke)./xqhat(rdke)./abs(Bax)).*Zcurr.x_0;
0812 %
0813 Zcurr.x_tp_fsav = 2*pi*(xqtilde(rdke)./xqbar(rdke))./(1 + xx0(rdke)/Rp).*(xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke).*XXlambda_p2m2p2(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0814 %************* to be updated YP 14/06/07 *************
0815 Zcurr.x_tp_vfsav = 2*pi*(xqtilde(rdke)./xqhat(rdke))./(1 + xx0(rdke)/Rp).*(xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke).*XXlambda_p2m2p2(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0816 Zcurr.x_tp_vcfsav = 2*pi*(xB0(rdke).*xqcronos(rdke)./xqhat(rdke)./Bax)./(1 + xx0(rdke)/Rp).*(xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke).*XXlambda_p2m2p2(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0817 %*************
0818 Zcurr.x_g_fsav = 2*pi*(xq(rdke)./xqbar(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_g(:,:,rdke).*XXheaviside(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0819 %************* to be updated YP 14/06/07 *************
0820 Zcurr.x_g_vfsav = 2*pi*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_g(:,:,rdke).*XXheaviside(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0821 Zcurr.x_g_vcfsav = 2*pi*(xB0(rdke).*xqcronos(rdke)./xqhat(rdke)./Bax).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_g(:,:,rdke).*XXheaviside(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0822 %*************
0823 Zcurr.x_tot_fsav = Zcurr.x_0_fsav + Zcurr.x_tp_fsav + Zcurr.x_g_fsav;%Total flux surface averaged current density (LUKE flow FSAV definition)
0824 Zcurr.x_tot_vfsav = Zcurr.x_0_vfsav + Zcurr.x_tp_vfsav + Zcurr.x_g_vfsav;%Total flux surface averaged current density (LUKE volumic FSAV definition)
0825 Zcurr.x_tot_vcfsav = Zcurr.x_0_vcfsav + Zcurr.x_tp_vcfsav + Zcurr.x_g_vcfsav;%Total flux surface averaged current density (CRONOS volumic FSAV definition)
0826 %
0827 if dkeparam.coll_mode == 3,%Lorentz model (test mode only, see documentation)
0828     Zcurr.x_L_tp_fsav = 2*pi*(xqtilde(rdke)./xqbar(rdke))./(1 + xx0(rdke)/Rp).*(xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_tp(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke).*XXlambda_p2m2p2(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0829     Zcurr.x_L_g_fsav = 2*pi*(xq(rdke)./xqbar(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0_L_g(:,:,rdke).*XXheaviside(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0830     Zcurr.x_L_tot = Zcurr.x_0_fsav + Zcurr.x_L_tp_fsav + Zcurr.x_L_g_fsav;%Total current density
0831 end
0832 %
0833 % Momentum transfer rate
0834 %
0835 % flux-suface averaging to be done
0836 %
0837 XXIntegrand_im = XXpnm(:,:,rdke).*XXpnm(:,:,rdke).*XXmhu(:,:,rdke);
0838 XXIntegrand_ip = XXpnp(:,:,rdke).*XXpnp(:,:,rdke).*XXmhu(:,:,rdke);
0839 %
0840 XXIntegrand_jm = XXpn(:,:,rdke).*XXpn(:,:,rdke).*sqrt(XX1mmhu2m);
0841 XXIntegrand_jp = XXpn(:,:,rdke).*XXpn(:,:,rdke).*sqrt(XX1mmhu2p);
0842 %
0843 Zmom.x_rf = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_im.*ZXXS.p0_rf_imj + XXIntegrand_ip.*ZXXS.p0_rf_ipj)/2,1),1).'...
0844            -2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXIntegrand_jm.*ZXXS.m0_rf_ijm + XXIntegrand_jp.*ZXXS.m0_rf_ijp)/2,1),1).';
0845 %
0846 %  Wave-induced radial transport by adjoint technique (P. Helander model)
0847 %
0848 xX = ones(nr,1)*pn/sqrt(2);
0849 xXG = (erf(xX) - 2*xX.*exp(-xX.^2)/sqrt(pi))./2./xX./xX;
0850 XXnhuDei = permute(repmat(((xrnhuth.*xZeff)'*ones(1,npn)).*((erf(xX) - xXG)./xX./xX./xX),[1,1,nmhu]),[2,3,1]);
0851 %
0852 if dke_mode == 1,
0853     xggradpsi22 = 2*pi*(xqtilde(rdke)./xqhat(rdke)).*(xR0(rdke).*xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,rdke).*XXnhuDei(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke).*XXlambda_b_p1m1(:,:,rdke),1),1)';
0854 end
0855 %
0856 %   Wave-induced radial transport by adjoint technique (P. Helander model)
0857 %
0858 if dke_mode == 1,
0859     if ny > 0,
0860          xggradpsi = xggradpsi1 - xggradpsi21 - xggradpsi22;
0861          xVr = xggradpsi.*xR0(rdke).*xBp0(rdke)./xne_norm(rdke);
0862          ZWiTransp_fsav{1} = xggradpsi;%For all rays
0863          ZWiTransp_fsav{2} = xggradpsi1;%For all rays
0864          ZWiTransp_fsav{3} = xggradpsi21;%For all rays
0865          ZWiTransp_fsav{4} = xggradpsi22;%For all rays
0866     else
0867          xggradpsi = zeros(length(rdke),1);
0868          xVr = zeros(length(rdke),1);
0869          ZWiTransp_fsav{1} = zeros(length(rdke),1);%For all rays
0870          ZWiTransp_fsav{2} = zeros(length(rdke),1);%For all rays
0871          ZWiTransp_fsav{3} = zeros(length(rdke),1);%For all rays
0872          ZWiTransp_fsav{4} = zeros(length(rdke),1);%For all rays
0873     end    
0874 end
0875 %
0876 % Run-away loss rate (Primary and secondary (avalanche) generation)
0877 %
0878 xepsi_rdke = xepsi(rdke);
0879 Xxqtilde = repmat(xqtilde,npn,1);
0880 Xxqhat = repmat(xqtilde,npn,1);
0881 %
0882 XXmaskm = XXmhu < 0;
0883 XXmaskp = XXmhu > 0;
0884 %
0885 XxRRm_fsav = 2*pi*(Xxqtilde(:,rdke)./Xxqhat(:,rdke)).*integral_dke_jd(dmhu,XXmaskm.*XXSp0_imj.*XXpn2m(:,:,rdke).*XXlambda(:,:,rdke),2);%Primary generation run-away rate (Dreicer field effect)
0886 XxRRp_fsav = 2*pi*(Xxqtilde(:,rdke)./Xxqhat(:,rdke)).*integral_dke_jd(dmhu,XXmaskp.*XXSp0_imj.*XXpn2m(:,:,rdke).*XXlambda(:,:,rdke),2);%Primary generation run-away rate (Dreicer field effect)
0887 XxRR_fsav = XxRRm_fsav + XxRRp_fsav;
0888 %
0889 xRRm_primary_fsav = XxRRm_fsav(end,:);%Radial profile of primary runaway electron generation rate
0890 xRRp_primary_fsav = XxRRp_fsav(end,:);%Radial profile of primary runaway electron generation rate
0891 %
0892 mask_xepsi = find(xepsi_rdke~=0);
0893 if ~isempty(mask_xepsi)
0894     Zcurr.x_B_fsav = Zcurr.x_tot_fsav - 0.5*XxRR_fsav(npn-1,mask_xepsi)*max(v2)./xepsi_rdke(mask_xepsi);%Bulk current as defined by Karney (useful for non-zero large dc electric field only)
0895 else
0896     Zcurr.x_B_fsav = Zcurr.x_tot_fsav;
0897 end 
0898 %
0899 % Knock-on runaway generation
0900 %
0901 %pnminRR = 2;%
0902 %zminRR = betath_ref*pnminRR;
0903 %zminRR2 = zminRR.*zminRR;
0904 %gammaminRR = sqrt(1 + zminRR2);
0905 %xRR_secondary_fsav = 1/(gammaminRR-1)*pi*xnloss(rdke).*xrnhuth(rdke)./xlnc_e(rdke);
0906 %
0907 xRRm_fsav = xRRm_primary_fsav;% + xRRm_secondary_fsav;%Runaway rate profile (primary + secondary contribution)
0908 xRRp_fsav = xRRp_primary_fsav;% + xRRp_secondary_fsav;%Runaway rate profile (primary + secondary contribution)
0909 %
0910 % Runaway fraction and current through internal boundary
0911 %
0912 Znorm.x_Rfrac_0 = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,rdke).*XXRintmask(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0913 Znorm.x_Rfrac_0_fsav = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,rdke).*XXf0(:,:,rdke).*XXRintmask(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0914 %
0915 Zcurr.x_Rfrac_0 = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,rdke).*XXRintmask(:,:,rdke).*XXpn2(:,:,rdke).*XXpn(:,:,rdke).*XXmhu(:,:,rdke)./XXgamma(:,:,rdke),1),1)';
0916 Zcurr.x_Rfrac_0_fsav = (xq(rdke)./xqbar(rdke)).*Zcurr.x_Rfrac_0;
0917 %
0918 % Equivalent output temperature
0919 %
0920 Xxf0 = 0.5*integral_dke_jd(dmhu,XXf0(:,:,rdke),2);
0921 XxfM = squeeze(XXfM(:,1,rdke));
0922 Xxpn = repmat(pn',[1,nr_dke]);
0923 Xxv = repmat(v',[1,nr_dke]);
0924 %
0925 xvth20 = 4*pi*integral_dke_jd(dpn,Xxpn.*Xxpn.*Xxv.*Xxv.*Xxf0,1)'/3;
0926 xvth2M = 4*pi*integral_dke_jd(dpn,Xxpn.*Xxpn.*Xxv.*Xxv.*XxfM,1)'/3;
0927 xTe_norm_out = (xvth20./xvth2M).*xTe_norm(rdke);%New bulk temperature (bulk heating)
0928 %
0929 %Exact trapped fraction (not the effective trapped fraction that rely on simplified plasma modeling with Maxwellian background)
0930 %
0931 if bounce_mode == 0,
0932     xft_fsav= zeros(1,nr_dke);%Fraction of trapped electrons
0933 elseif bounce_mode == 1,
0934     xne_trap_fsav = integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXf0(:,:,rdke) + XXf0_g(:,:,rdke)).*~XXheaviside(:,:,rdke).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke),1),1)';%trapped region only
0935     xne_tot_fsav = integral_dke_jd(dmhu,integral_dke_jd(dpn,(XXf0(:,:,rdke) + XXf0_tp(:,:,rdke) + XXf0_g(:,:,rdke)).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke),1),1)';%
0936     xft_fsav = xne_trap_fsav./xne_tot_fsav;
0937 end
0938 %
0939 %Effective trapped fraction
0940 %
0941 if bounce_mode == 0,
0942     xfteff_fsav= zeros(1,nr_dke);%Effective fraction of trapped electrons
0943 elseif bounce_mode == 1,
0944     xfteff_fsav = 1.5*(xBT0(rdke)./xB0(rdke)).*((xqtilde(rdke)./xqbar(rdke)).*(xBT0(rdke)./xB0(rdke)).*integral_dke_jd(dmhu,XXmhu(1,:,rdke).*XXmhu(1,:,rdke).*XXlambda_p2m2p2(1,:,rdke),2));    
0945     xfteff_fsav = xfteff_fsav - 1.5*(xBT0(rdke)./xB0(rdke)).*((xq(rdke)./xqbar(rdke)).*(1 + xx0(rdke)/Rp).*integral_dke_jd(dmhu,abs(XXmhu(1,:,rdke)).*XXheaviside(1,:,rdke).*abs(XXILor(1,:,rdke)),2));
0946 end
0947 %
0948 %   Magnetic ripple loss rate calculation (for Tore Supra or WEST tokamak only)
0949 %   Volume or surface loss calculations must match if no additional source term
0950 %   is added (XXlossripple = XXsource). Otherwise, only the surface loss must be considered
0951 %
0952 if sum(sum(sum(XXlossripple))) ~= 0 & (~isempty(findstr(equilDKE.tokname,'TS')) || ~isempty(findstr(equilDKE.tokname,'WEST'))),
0953     %
0954     XXEc = sqrt(XXpn2*mksa.betath_ref^2 + 1) - 1;
0955     XXEcm = sqrt(XXpn2m*mksa.betath_ref^2 + 1) - 1;
0956     %
0957     xMRR_tau = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlossripple(:,:,rdke).*XXf0(:,:,rdke).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke),1),1)';%Flux of electrons leaving the plasma through the loss cone
0958     xMRR_power_tau = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlossripple(:,:,rdke).*XXf0(:,:,rdke).*XXpn2(:,:,rdke).*XXEc(:,:,rdke).*XXlambda(:,:,rdke),1),1)';%Flux of electrons leaving the plasma through the loss cone
0959     %
0960     XXmaskMRR_im = [zeros(1,nmhu,nr);abs(diff(XXlossripple~=0))];%Mask in p (only nhuloss ~=0 is considered)
0961     XXmaskMRR_jm = [zeros(npn,1,nr),abs(permute(diff(permute(XXlossripple~=0,[2,1,3])),[2,1,3]))];%Mask in mhu (mhu < 0 only because of symmetry) (only nhuloss ~=0 is considered)
0962     XXmaskMRR_jm(:,mhu>0,:)=0;%mhu < 0 only because of symmetry)
0963     XxMRR_flux_p = 2*pi*integral_dke_jd(dmhum,XXmaskMRR_im(:,:,rdke).*XXSp0_imj.*XXpn2m(:,:,rdke).*XXlambda(:,:,rdke).*XXRlambda_m(:,:,rdke),2);
0964     XxMRR_flux_mhu = 2*2*pi*integral_dke_jd(dpnm,XXmaskMRR_jm(:,:,rdke).*XXSm0_ijm.*XXpnm(:,:,rdke).*sqrt(XX1mmhu2m(:,:,rdke)).*XXlambda(:,:,rdke).*XXRlambda_m(:,:,rdke),1);
0965     xMRR_flux_p = XxMRR_flux_p(XxMRR_flux_p~=0)';
0966     xMRR_flux_mhu = XxMRR_flux_mhu(XxMRR_flux_mhu~=0)';
0967     xMRR_flux = xMRR_flux_p - xMRR_flux_mhu;%Reduce number of points since only nhuloss ~=0 is considered
0968     xMRR_flux =[xMRR_tau(find(xMRR_tau==0)),xMRR_flux];%Restore the points for nhuloss =0
0969 
0970     XxMRR_power_flux_p = 2*pi*integral_dke_jd(dmhum,XXmaskMRR_im(:,:,rdke).*XXSp0_imj.*XXpn2m(:,:,rdke).*XXEcm(:,:,rdke).*XXlambda(:,:,rdke).*XXRlambda_m(:,:,rdke),2);
0971     XxMRR_power_flux_mhu = 2*2*pi*integral_dke_jd(dpnm,XXmaskMRR_jm(:,:,rdke).*XXSm0_ijm.*XXpnm(:,:,rdke).*XXEcm(:,:,rdke).*sqrt(XX1mmhu2m(:,:,rdke)).*XXlambda(:,:,rdke).*XXRlambda_m(:,:,rdke),1);
0972     xMRR_power_flux_p = XxMRR_power_flux_p(XxMRR_power_flux_p~=0)';
0973     xMRR_power_flux_mhu = XxMRR_power_flux_mhu(XxMRR_power_flux_mhu~=0)';
0974     xMRR_power_flux = xMRR_power_flux_p - xMRR_power_flux_mhu;%Reduce number of points since only nhuloss ~=0 is considered
0975     xMRR_power_flux =[xMRR_power_tau(find(xMRR_power_tau==0)),xMRR_power_flux];%Restore the points for nhuloss =0
0976 else
0977     xMRR_tau = NaN*ones(1,nr_dke);
0978     xMRR_flux = NaN*ones(1,nr_dke);
0979     xMRR_power_flux = NaN*ones(1,nr_dke);
0980     xMRR_power_tau = NaN*ones(1,nr_dke);
0981 end
0982 %
0983 xne_norm_out = 2*pi.*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXsinksource(:,:,rdke).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke),1),1)';%transfer rate trought sink/source term
0984 %
0985 Znorm.xRRAtot = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXSavalanches(:,:,rdke).*XXpn2(:,:,rdke),1),1)';
0986 Znorm.xRRAint = 2*pi*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXSavalanches(:,:,rdke).*XXpn2(:,:,rdke).*XXRintmask(:,:,rdke),1),1)';
0987 %
0988 Znorm.xRRAtot_fsav = 2*pi*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXSavalanches(:,:,rdke).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke),1),1)';
0989 Znorm.xRRAint_fsav = 2*pi*(xqtilde(rdke)./xqhat(rdke)).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXSavalanches(:,:,rdke).*XXpn2(:,:,rdke).*XXlambda(:,:,rdke).*XXRintmask(:,:,rdke),1),1)';
0990 %
0991 % 0-D moments of the distribution
0992 %
0993 Zcurr.I_tot = sum(Zcurr.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref);
0994 Zcurr.I_Rfrac_tot = sum(Zcurr.x_Rfrac_0_fsav.*equilDKE.xdA_dke*mksa.j_ref);
0995 ZP0.P_rf_2piRp = sum(ZP0.x_rf_fsav.'.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0996 ZP0.P_ohm_2piRp = sum(ZP0.x_e_fsav.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0997 %
0998 if display_mode >= 1,info_dke_yp(2,['Moments of the distribution function calculated in ',num2str(etime(clock,time0)),' (s)']);end
0999 %
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.