matrixcalc_dke_yp

   Calculate the super matrix for the 3D electron relativistic drift kinetic solver.

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

0001 function [Zdelta,MMXT_f_t,MMXP_f_t,MMXR_f_t,Xmask_r_edge,MMXP_tp,MMXP_g_t,MMH_tp,Zpn2_t,sm_f,sm_g,sm_tp] = matrixcalc_dke_yp(dkeparam,display,gridDKE,equilDKE,Zbouncecoef,mksa,Zmripple,Zmomcoef,transpfaste,...
0002                                                                        ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,XXILor,...
0003                                                                        ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,...
0004                                                                        ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,...
0005                                                                        ZXXD_rf,ZXXF_rf,ZXXD_rf_tp,ZXXF_rf_tp,...
0006                                                                        ZXXD_a,ZXXF_a)
0007 %
0008 %   Calculate the super matrix for the 3D electron relativistic drift kinetic solver.
0009 %
0010 %   by Y. Peysson (CEA-IRFM) (yves.peysson@cea.fr) and J. Decker (CEA-IRFM) (joan.decker@cea.fr)
0011 %
0012 time0 = clock;
0013 etime_fp = 0;
0014 %
0015 boundary_mode_f = dkeparam.boundary_mode_f;
0016 dke_mode = dkeparam.dke_mode;
0017 bounce_mode = dkeparam.bounce_mode;
0018 dtn = dkeparam.dtn;
0019 %
0020 display_mode = display.display_mode;
0021 %
0022 nr = gridDKE.nr;
0023 npn = gridDKE.npn;
0024 nmhu = gridDKE.nmhu;
0025 nr_dke = gridDKE.nr_dke;
0026 mhu = gridDKE.mhu;
0027 Xdmhumm = gridDKE.Xdmhumm;
0028 Xdmhum = gridDKE.Xdmhum;
0029 Xdmhu = gridDKE.Xdmhu;
0030 Xdmhup = gridDKE.Xdmhup;
0031 Xdmhupp = gridDKE.Xdmhupp;
0032 Xdpnmm = gridDKE.Xdpnmm;
0033 Xdpnm = gridDKE.Xdpnm;
0034 Xdpn = gridDKE.Xdpn;
0035 Xdpnp = gridDKE.Xdpnp;
0036 Xdpnpp = gridDKE.Xdpnpp;
0037 Xmhumm = gridDKE.Xmhumm;
0038 Xmhum = gridDKE.Xmhum;
0039 Xmhu = gridDKE.Xmhu;
0040 Xmhu2 = gridDKE.Xmhu2;
0041 Xmhup = gridDKE.Xmhup;
0042 Xmhupp = gridDKE.Xmhupp;
0043 X1mmhu2m = gridDKE.X1mmhu2m;
0044 X1mmhu2 = gridDKE.X1mmhu2;
0045 X1mmhu2p = gridDKE.X1mmhu2p;
0046 pn = gridDKE.pn;
0047 Xpnmm = gridDKE.Xpnmm;
0048 Xpnm = gridDKE.Xpnm;
0049 Xpn = gridDKE.Xpn;
0050 Xpnp = gridDKE.Xpnp;
0051 Xpnpp = gridDKE.Xpnpp;
0052 Xpn2m = gridDKE.Xpn2m;
0053 Xpn2 = gridDKE.Xpn2;
0054 Xpn2p = gridDKE.Xpn2p;
0055 XXpn2 = gridDKE.XXpn2;
0056 xpsin = gridDKE.xpsin;
0057 xdpsin = gridDKE.xdpsin_dke;
0058 xdpsinm = gridDKE.xdpsinm_dke;
0059 xdpsinp = gridDKE.xdpsinp_dke;
0060 xmhubounce2 = gridDKE.xmhubounce2;
0061 rdke = gridDKE.rdke;
0062 xrho = gridDKE.xrho;
0063 %
0064 gamma = Zmomcoef.gamma;
0065 Xvpar = Zmomcoef.Xvpar;
0066 %
0067 djmhut = Zbouncecoef.djmhut;
0068 Zmask_f0_P = Zbouncecoef.Zmask_f0_P;
0069 Zbouncecoef.Zmask_f0_P = Zmask_f0_P;
0070 Zmask_f0_Tm = Zbouncecoef.Zmask_f0_Tm;
0071 Zmask_f0_Tp = Zbouncecoef.Zmask_f0_Tp;
0072 jmhut_min = Zbouncecoef.jmhut_min;
0073 jmhut_max = Zbouncecoef.jmhut_max;
0074 Zmask_f0_t = Zbouncecoef.Zmask_f0_t;
0075 Zmask_f0_t_P = Zbouncecoef.Zmask_f0_t_P;
0076 Zmask_f0_t_Tp = Zbouncecoef.Zmask_f0_t_Tp;
0077 ns_f = Zbouncecoef.ns_f;    
0078 XXlambda = Zbouncecoef.XXlambda;
0079 XXRlambda_m = Zbouncecoef.XXRlambda_m;
0080 XXRlambda_p = Zbouncecoef.XXRlambda_p;
0081 sm_f = Zbouncecoef.sm_f;
0082 sm_T_f = Zbouncecoef.sm_T_f;
0083 blocksize_f_t = Zbouncecoef.blocksize_f_t;
0084 blocksize_T_f_t = Zbouncecoef.blocksize_T_f_t;
0085 xqtilde = Zbouncecoef.xqtilde;
0086 %
0087 xx0 = equilDKE.xx0;
0088 xBp0 = equilDKE.xBp0;
0089 xB0 = equilDKE.xB0;
0090 ap = equilDKE.ap;
0091 Rp = equilDKE.Rp;
0092 psia_apRp = equilDKE.psia_apRp;
0093 %
0094 betath_ref = mksa.betath_ref;
0095 nhu_ref = mksa.nhu_ref;
0096 %
0097 bounce_mode = dkeparam.bounce_mode;    
0098 %
0099 %if dke_mode == 1,
0100 %Zmask_g = Zbouncecoef.Zmask_g;
0101 %ns_g = Zbouncecoef.ns_g;
0102 %xftp_norm = mksa.xftp_norm_ref;
0103 %boundary_mode_g = dkeparam.boundary_mode_g;
0104 sm_tp = Zbouncecoef.sm_tp;
0105 sm_g = Zbouncecoef.sm_g;
0106 %blocksize_tp = Zbouncecoef.blocksize_tp;
0107 %blocksize_g_t = Zbouncecoef.blocksize_g_t;
0108 %
0109 %end
0110 %
0111 MMXT_f_t = [];
0112 MMXP_f_t = [];
0113 MMXR_f_t = [];
0114 %
0115 Xmask_r_edge = [];
0116 %
0117 nMMXP_f_t = max(cumsum(Zbouncecoef.blocksize_f_t));
0118 %
0119 MMXP_tp = [];
0120 MMXP_g_t = [];
0121 MMH_tp = [];
0122 %
0123 spparms('default');%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions only
0124 spparms('supernd',4);%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions only
0125 spparms('rreduce',4);%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions onlyXXdpsin
0126 %
0127 % Case where radial transport is set to be proportional to RF diffusion coefficient
0128 %
0129 XXD_rf_norm = ZXXD_rf.ij./max(ZXXD_rf.ij(:));
0130 %
0131 if isfield(transpfaste,'Dr_model') && imag(transpfaste.Dr_model) == 1,
0132     fnames = fieldnames(ZXXD_a);
0133     for ifield = 1:length(fnames),
0134         ZXXD_a.(fnames{ifield}) = ZXXD_a.(fnames{ifield}).*XXD_rf_norm;
0135     end
0136 end
0137 %
0138 if isfield(transpfaste,'Vr_model') && imag(transpfaste.Vr_model) == 1,
0139     fnames = fieldnames(ZXXF_a);
0140     for ifield = 1:length(fnames),
0141         ZXXF_a.(fnames{ifield}) = ZXXF_a.(fnames{ifield}).*XXD_rf_norm;
0142     end
0143 end   
0144 %
0145 %TODO: remove transpfaste dependence when collisional transport is set
0146 %
0147 if ~isempty(transpfaste),%boundary conditions
0148     Xmask_r_edge = (ZXXD_a.rr_lp(:,:,nr) > 0);
0149 end
0150 %
0151 % Build momentum supermatrix MMXP_f_t
0152 %
0153 for ir = 1:nr,
0154     %
0155     XDpp_ipj = ZXXD_c.pp_ipj(:,:,ir) + ZXXD_e.pp_ipj(:,:,ir) + ZXXD_s.pp_ipj(:,:,ir) + ZXXD_rf.pp_ipj(:,:,ir) + ZXXD_a.pp_ipj(:,:,ir);
0156     XDpp_imj = ZXXD_c.pp_imj(:,:,ir) + ZXXD_e.pp_imj(:,:,ir) + ZXXD_s.pp_imj(:,:,ir) + ZXXD_rf.pp_imj(:,:,ir) + ZXXD_a.pp_imj(:,:,ir);
0157     XDpm_ipj = ZXXD_c.pm_ipj(:,:,ir) + ZXXD_e.pm_ipj(:,:,ir) + ZXXD_s.pm_ipj(:,:,ir) + ZXXD_rf.pm_ipj(:,:,ir) + ZXXD_a.pm_ipj(:,:,ir);
0158     XDpm_imj = ZXXD_c.pm_imj(:,:,ir) + ZXXD_e.pm_imj(:,:,ir) + ZXXD_s.pm_imj(:,:,ir) + ZXXD_rf.pm_imj(:,:,ir) + ZXXD_a.pm_imj(:,:,ir);
0159     XFp_ipj = ZXXF_c.p_ipj(:,:,ir) + ZXXF_e.p_ipj(:,:,ir) + ZXXF_s.p_ipj(:,:,ir) + ZXXF_rf.p_ipj(:,:,ir) + ZXXF_a.p_ipj(:,:,ir);
0160     XFp_imj = ZXXF_c.p_imj(:,:,ir) + ZXXF_e.p_imj(:,:,ir) + ZXXF_s.p_imj(:,:,ir) + ZXXF_rf.p_imj(:,:,ir) + ZXXF_a.p_imj(:,:,ir);
0161     %
0162     XDmm_ijp = ZXXD_c.mm_ijp(:,:,ir) + ZXXD_e.mm_ijp(:,:,ir) + ZXXD_s.mm_ijp(:,:,ir) + ZXXD_rf.mm_ijp(:,:,ir) + ZXXD_a.mm_ijp(:,:,ir);
0163     XDmm_ijm = ZXXD_c.mm_ijm(:,:,ir) + ZXXD_e.mm_ijm(:,:,ir) + ZXXD_s.mm_ijm(:,:,ir) + ZXXD_rf.mm_ijm(:,:,ir) + ZXXD_a.mm_ijm(:,:,ir);
0164     XDmp_ijp = ZXXD_c.mp_ijp(:,:,ir) + ZXXD_e.mp_ijp(:,:,ir) + ZXXD_s.mp_ijp(:,:,ir) + ZXXD_rf.mp_ijp(:,:,ir) + ZXXD_a.mp_ijp(:,:,ir);
0165     XDmp_ijm = ZXXD_c.mp_ijm(:,:,ir) + ZXXD_e.mp_ijm(:,:,ir) + ZXXD_s.mp_ijm(:,:,ir) + ZXXD_rf.mp_ijm(:,:,ir) + ZXXD_a.mp_ijm(:,:,ir);
0166     XFm_ijp = ZXXF_c.m_ijp(:,:,ir) + ZXXF_e.m_ijp(:,:,ir) + ZXXF_s.m_ijp(:,:,ir) + ZXXF_rf.m_ijp(:,:,ir) + ZXXF_a.m_ijp(:,:,ir);
0167     XFm_ijm = ZXXF_c.m_ijm(:,:,ir) + ZXXF_e.m_ijm(:,:,ir) + ZXXF_s.m_ijm(:,:,ir) + ZXXF_rf.m_ijm(:,:,ir) + ZXXF_a.m_ijm(:,:,ir);
0168     %
0169     XDpp_ij = ZXXD_c.pp_ij(:,:,ir) + ZXXD_e.pp_ij(:,:,ir) + ZXXD_s.pp_ij(:,:,ir) + ZXXD_rf.pp_ij(:,:,ir) + ZXXD_a.pp_ij(:,:,ir);
0170     XDpm_ij = ZXXD_c.pm_ij(:,:,ir) + ZXXD_e.pm_ij(:,:,ir) + ZXXD_s.pm_ij(:,:,ir) + ZXXD_rf.pm_ij(:,:,ir) + ZXXD_a.pm_ij(:,:,ir);
0171     XDmp_ij = ZXXD_c.mp_ij(:,:,ir) + ZXXD_e.mp_ij(:,:,ir) + ZXXD_s.mp_ij(:,:,ir) + ZXXD_rf.mp_ij(:,:,ir) + ZXXD_a.mp_ij(:,:,ir);
0172     XFp_ij = ZXXF_c.p_ij(:,:,ir) + ZXXF_e.p_ij(:,:,ir) + ZXXF_s.p_ij(:,:,ir) + ZXXF_rf.p_ij(:,:,ir) + ZXXF_a.p_ij(:,:,ir);
0173 %
0174 %Grid weights calculations (for f0 and ftp also)
0175 %
0176     [Xdeltap_imj,Xdeltap_ipj,...
0177          Xdeltam_ijm,Xdeltam_ijp,...
0178          Xdeltap_imjmm,Xdeltap_imjpp,Xdeltap_ipjmm,Xdeltap_ipjpp,...
0179          Xdeltam_immjm,Xdeltam_ippjm,Xdeltam_immjp,Xdeltam_ippjp,...
0180          xdeltap_npj,xdeltap_npjmm,xdeltap_npjpp] = ...
0181          fppgridweights_dke_yp(Xdpnmm,Xdpnm,Xdpn,Xdpnp,Xdpnpp,...
0182          Xdmhumm,Xdmhum,Xdmhu,Xdmhup,Xdmhupp,...     
0183          XDpp_ipj,XDpp_imj,XFp_ipj,XFp_imj,...
0184          XDpp_ij,XFp_ij);
0185 %
0186     XXdeltap_imj(:,:,ir) = Xdeltap_imj; 
0187     XXdeltap_ipj(:,:,ir) = Xdeltap_ipj;
0188     XXdeltam_ijp(:,:,ir) = Xdeltam_ijp;
0189     XXdeltam_ijm(:,:,ir) = Xdeltam_ijm;   
0190     XXdeltap_imjmm(:,:,ir) = Xdeltap_imjmm;
0191     XXdeltap_imjpp(:,:,ir) = Xdeltap_imjpp;
0192     XXdeltap_ipjmm(:,:,ir) = Xdeltap_ipjmm;
0193     XXdeltap_ipjpp(:,:,ir) = Xdeltap_ipjpp;
0194     XXdeltam_immjm(:,:,ir) = Xdeltam_immjm;
0195     XXdeltam_ippjm(:,:,ir) = Xdeltam_ippjm;
0196     XXdeltam_immjp(:,:,ir) = Xdeltam_immjp;
0197     XXdeltam_ippjp(:,:,ir) = Xdeltam_ippjp;
0198 %
0199 %Momentum dynamics matrix coefficients calculations
0200 %
0201     [XXa(:,:,ir),XXb(:,:,ir),XXc(:,:,ir),XXd(:,:,ir),XXe(:,:,ir),XXf(:,:,ir),XXg(:,:,ir),XXh(:,:,ir),XXi(:,:,ir)] = ...
0202     fpengine_dke_yp(dkeparam,0,...
0203     XXdeltap_imj(:,:,ir),XXdeltap_ipj(:,:,ir),XXdeltam_ijm(:,:,ir),XXdeltam_ijp(:,:,ir),...               
0204     XXdeltap_imjmm(:,:,ir),XXdeltap_imjpp(:,:,ir),XXdeltap_ipjmm(:,:,ir),XXdeltap_ipjpp(:,:,ir),...
0205     XXdeltam_immjm(:,:,ir),XXdeltam_ippjm(:,:,ir),XXdeltam_immjp(:,:,ir),XXdeltam_ippjp(:,:,ir),...
0206     XXRlambda_m(:,:,ir),XXRlambda_p(:,:,ir),...
0207     Xpnmm,Xpnm,Xpn,Xpnp,Xpnpp,Xpn2m,Xpn2p,Xdpnmm,Xdpnm,Xdpn,Xdpnp,Xdpnpp,...
0208     Xmhumm,Xmhum,Xmhu,Xmhup,Xmhupp,X1mmhu2,X1mmhu2m,X1mmhu2p,Xdmhumm,Xdmhum,Xdmhu,Xdmhup,Xdmhupp,...
0209     XDpp_ipj,XDpp_imj,XDpm_ipj,XDpm_imj,XFp_ipj,XFp_imj,...
0210     XDmm_ijp,XDmm_ijm,XDmp_ijp,XDmp_ijm,XFm_ijp,XFm_ijm,...
0211     XDpp_ij,XDpm_ij,XDmp_ij,XFp_ij);
0212 %
0213 %Sink and source coefficients like magnetic ripple losses,described by a Krook term
0214 %
0215     XXc(:,:,ir) = XXc(:,:,ir) +  Zmripple.XXlossripple(:,:,ir).*Xpn2;
0216 %
0217 %Boundaries of the bounce domain for f0
0218 %
0219     XXc(:,1,ir) = XXc(:,1,ir) + XXa(:,1,ir);%Because of cross-derivatives
0220     XXd(:,1,ir) = XXd(:,1,ir) + XXi(:,1,ir);%Because of cross-derivatives
0221     XXb(:,1,ir) = XXb(:,1,ir) + XXf(:,1,ir);%Because of cross-derivatives
0222 %
0223     XXa(:,1,ir) = 0;%Because of cross-derivatives
0224     XXi(:,1,ir) = 0;%Because of cross-derivatives
0225     XXf(:,1,ir) = 0;%Because of cross-derivatives
0226 %
0227     XXc(:,nmhu,ir) = XXc(:,nmhu,ir) + XXe(:,nmhu,ir);%Because of cross-derivatives
0228     XXd(:,nmhu,ir) = XXd(:,nmhu,ir) + XXg(:,nmhu,ir);%Because of cross-derivatives
0229     XXb(:,nmhu,ir) = XXb(:,nmhu,ir) + XXh(:,nmhu,ir);%Because of cross-derivatives
0230 %
0231     XXe(:,nmhu,ir) = 0;%Because of cross-derivatives
0232     XXg(:,nmhu,ir) = 0;%Because of cross-derivatives
0233     XXh(:,nmhu,ir) = 0;%Because of cross-derivatives
0234 %
0235     XXc(1,:,ir) = XXc(1,:,ir) + fliplr(XXb(1,:,ir));%Because of cross-derivatives
0236     XXa(1,:,ir) = XXa(1,:,ir) + fliplr(XXf(1,:,ir));%Because of cross-derivatives
0237     XXe(1,:,ir) = XXe(1,:,ir) + fliplr(XXh(1,:,ir));%Because of cross-derivatives
0238 %
0239     XXb(1,:,ir) = 0;%Because of cross-derivatives
0240     XXh(1,:,ir) = 0;%Because of cross-derivatives
0241     XXf(1,:,ir) = 0;%Because of cross-derivatives
0242 %
0243 %Edge boundary conditions for radial transport (Maxwellian solution enforced where radial transport is non-zero)
0244 %
0245     if ~isempty(Xmask_r_edge) & ir == nr,
0246         XXa(:,:,ir) = ~Xmask_r_edge.*XXa(:,:,ir);
0247         XXb(:,:,ir) = ~Xmask_r_edge.*XXb(:,:,ir);
0248         XXc(:,:,ir) = ~Xmask_r_edge.*XXc(:,:,ir);
0249         XXd(:,:,ir) = ~Xmask_r_edge.*XXd(:,:,ir);
0250         XXe(:,:,ir) = ~Xmask_r_edge.*XXe(:,:,ir);
0251         XXf(:,:,ir) = ~Xmask_r_edge.*XXf(:,:,ir);
0252         XXg(:,:,ir) = ~Xmask_r_edge.*XXg(:,:,ir);
0253         XXh(:,:,ir) = ~Xmask_r_edge.*XXh(:,:,ir);
0254         XXi(:,:,ir) = ~Xmask_r_edge.*XXi(:,:,ir);
0255     end
0256 %
0257     if bounce_mode == 0,%No bounce averaging
0258         Xa_t = XXa(:,:,ir);
0259         Xb_t = XXb(:,:,ir);
0260         Xc_t = XXc(:,:,ir);
0261         Xd_t = XXd(:,:,ir);
0262         Xe_t = XXe(:,:,ir);
0263         Xf_t = XXf(:,:,ir);
0264         Xg_t = XXg(:,:,ir);
0265         Xh_t = XXh(:,:,ir);
0266         Xi_t = XXi(:,:,ir);              
0267         mhu_t = mhu;
0268         Xpn2_t = Xpn2;
0269     elseif bounce_mode == 1,%Full implicit bounce-averaging calculations, fifteen diagonals mode
0270          if djmhut(ir) > 3,%At least three points needed in the bounce domain
0271             Xa_t = zeros(npn,nmhu - djmhut(ir));
0272             Xb_t = zeros(npn,nmhu - djmhut(ir));
0273             Xc_t = zeros(npn,nmhu - djmhut(ir));
0274             Xd_t = zeros(npn,nmhu - djmhut(ir));
0275             Xe_t = zeros(npn,nmhu - djmhut(ir));
0276             Xf_t = zeros(npn,nmhu - djmhut(ir));
0277             Xg_t = zeros(npn,nmhu - djmhut(ir));
0278             Xh_t = zeros(npn,nmhu - djmhut(ir));
0279             Xi_t = zeros(npn,nmhu - djmhut(ir));
0280             Xj_t = zeros(npn,nmhu - djmhut(ir));
0281             Xk_t = zeros(npn,nmhu - djmhut(ir));
0282             Xl_t = zeros(npn,nmhu - djmhut(ir));
0283             Xm_t = zeros(npn,nmhu - djmhut(ir));
0284             Xn_t = zeros(npn,nmhu - djmhut(ir));
0285             Xo_t = zeros(npn,nmhu - djmhut(ir));
0286 %
0287             Xa_t(:,Zmask_f0_t_P{ir}) = XXa(:,Zmask_f0_P{ir},ir);
0288             Xb_t(:,Zmask_f0_t_P{ir}) = XXb(:,Zmask_f0_P{ir},ir);
0289             Xc_t(:,Zmask_f0_t_P{ir}) = XXc(:,Zmask_f0_P{ir},ir);
0290             Xd_t(:,Zmask_f0_t_P{ir}) = XXd(:,Zmask_f0_P{ir},ir);
0291             Xe_t(:,Zmask_f0_t_P{ir}) = XXe(:,Zmask_f0_P{ir},ir);
0292             Xf_t(:,Zmask_f0_t_P{ir}) = XXf(:,Zmask_f0_P{ir},ir);
0293             Xg_t(:,Zmask_f0_t_P{ir}) = XXg(:,Zmask_f0_P{ir},ir);
0294             Xh_t(:,Zmask_f0_t_P{ir}) = XXh(:,Zmask_f0_P{ir},ir);
0295             Xi_t(:,Zmask_f0_t_P{ir}) = XXi(:,Zmask_f0_P{ir},ir);             
0296 %
0297 %Symmetrization of the trapped region
0298 %
0299             Xa_t(:,Zmask_f0_t_Tp{ir}) = (XXa(:,Zmask_f0_Tp{ir},ir) + fliplr(XXe(:,Zmask_f0_Tm{ir},ir)))/2;
0300             Xb_t(:,Zmask_f0_t_Tp{ir}) = (XXb(:,Zmask_f0_Tp{ir},ir) + fliplr(XXb(:,Zmask_f0_Tm{ir},ir)))/2;
0301             Xc_t(:,Zmask_f0_t_Tp{ir}) = (XXc(:,Zmask_f0_Tp{ir},ir) + fliplr(XXc(:,Zmask_f0_Tm{ir},ir)))/2;
0302             Xd_t(:,Zmask_f0_t_Tp{ir}) = (XXd(:,Zmask_f0_Tp{ir},ir) + fliplr(XXd(:,Zmask_f0_Tm{ir},ir)))/2;
0303             Xe_t(:,Zmask_f0_t_Tp{ir}) = (XXe(:,Zmask_f0_Tp{ir},ir) + fliplr(XXa(:,Zmask_f0_Tm{ir},ir)))/2;
0304             Xf_t(:,Zmask_f0_t_Tp{ir}) = (XXf(:,Zmask_f0_Tp{ir},ir) + fliplr(XXh(:,Zmask_f0_Tm{ir},ir)))/2;
0305             Xg_t(:,Zmask_f0_t_Tp{ir}) = (XXg(:,Zmask_f0_Tp{ir},ir) + fliplr(XXi(:,Zmask_f0_Tm{ir},ir)))/2;
0306             Xh_t(:,Zmask_f0_t_Tp{ir}) = (XXh(:,Zmask_f0_Tp{ir},ir) + fliplr(XXf(:,Zmask_f0_Tm{ir},ir)))/2;
0307             Xi_t(:,Zmask_f0_t_Tp{ir}) = (XXi(:,Zmask_f0_Tp{ir},ir) + fliplr(XXg(:,Zmask_f0_Tm{ir},ir)))/2;
0308 %
0309             Xh_t(:,jmhut_max(ir) - djmhut(ir)) = Xh_t(:,jmhut_max(ir) - djmhut(ir))/2;
0310             Xe_t(:,jmhut_max(ir) - djmhut(ir)) = Xe_t(:,jmhut_max(ir) - djmhut(ir))/2;
0311             Xg_t(:,jmhut_max(ir) - djmhut(ir)) = Xg_t(:,jmhut_max(ir) - djmhut(ir))/2;
0312 %
0313             Xj_t(:,jmhut_max(ir) - djmhut(ir)) = Xh_t(:,jmhut_max(ir) - djmhut(ir));%f0(i-1/2,j+1/2-djmhut)
0314             Xk_t(:,jmhut_max(ir) - djmhut(ir)) = Xe_t(:,jmhut_max(ir) - djmhut(ir));%f0(i+1/2,j+1/2-djmhut)
0315             Xl_t(:,jmhut_max(ir) - djmhut(ir)) = Xg_t(:,jmhut_max(ir) - djmhut(ir));%f0(i+3/2,j+1/2-djmhut)
0316 %
0317             Xm_t(:,jmhut_min(ir) - 1) = XXh(:,jmhut_min(ir) - 1,ir);%f0(i-1/2,j+1/2+djmhut)
0318             Xn_t(:,jmhut_min(ir) - 1) = XXe(:,jmhut_min(ir) - 1,ir);%f0(i+1/2,j+1/2+djmhut)
0319             Xo_t(:,jmhut_min(ir) - 1) = XXg(:,jmhut_min(ir) - 1,ir);%f0(i+3/2,j+1/2+djmhut)
0320 %
0321             Xh_t(:,jmhut_min(ir) - 1) = 0;
0322             Xe_t(:,jmhut_min(ir) - 1) = 0;
0323             Xg_t(:,jmhut_min(ir) - 1) = 0;
0324 %
0325             Xa_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0326             Xb_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0327             Xc_t(:,nmhu/2 + 1 - djmhut(ir)) = 1;%df0(pn,mhu=pi/2) = 0
0328             Xd_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0329             Xe_t(:,nmhu/2 + 1 - djmhut(ir)) = -1;%df0(pn,mhu=pi/2) = 0
0330             Xf_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0331             Xg_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0332             Xh_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0333             Xi_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0334             Xj_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0335             Xk_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0336             Xl_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0337             Xm_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0338             Xn_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0339             Xo_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0340 %
0341             mhu_t = mhu(Zmask_f0_t{ir});
0342             Xpn2_t = Xpn2(:,Zmask_f0_t{ir});           
0343          else
0344             error(['WARNING: Not enough pitch-angle grid points in the bounce region for accurate bounce-averaged calculations (increasing nmhu may solve the problem) !'])
0345         end
0346     end
0347     %
0348     %Local boundary conditions at p = p(1/2): not necessary but help for the density conservation (Dirichlet boundary condition)
0349     %
0350     if boundary_mode_f > 0,
0351         Xa_t(1:boundary_mode_f,:) = 0;
0352         Xb_t(1:boundary_mode_f,:) = 0;
0353         Xc_t(1:boundary_mode_f,:) = 0;%Maxwellian solution enforced at p = p(1/2 -> boundary_mode_f-1/2) (More accurate but not fully necessary)
0354         Xd_t(1:boundary_mode_f,:) = 0;
0355         Xe_t(1:boundary_mode_f,:) = 0;
0356         Xf_t(1:boundary_mode_f,:) = 0;
0357         Xg_t(1:boundary_mode_f,:) = 0;
0358         Xh_t(1:boundary_mode_f,:) = 0;
0359         Xi_t(1:boundary_mode_f,:) = 0;
0360         if bounce_mode == 1,
0361             Xj_t(1:boundary_mode_f,:) = 0;
0362             Xk_t(1:boundary_mode_f,:) = 0;
0363             Xl_t(1:boundary_mode_f,:) = 0;
0364             Xm_t(1:boundary_mode_f,:) = 0;
0365             Xn_t(1:boundary_mode_f,:) = 0;
0366             Xo_t(1:boundary_mode_f,:) = 0;   
0367 %
0368             Xpn2_t(boundary_mode_f+1:npn,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0 only for p values where the Maxwellian solution is not enforced
0369         end 
0370     else
0371         if bounce_mode == 1,
0372             Xpn2_t(:,nmhu/2 + 1 - djmhut(ir)) = 0;%df0(pn,mhu=pi/2) = 0
0373         end
0374     end
0375 %
0376     Zpn2_t{ir} = Xpn2_t;    
0377 %
0378 %Momentum flux matrix calculation
0379 %
0380     if bounce_mode == 0,
0381         MXc = reshape(Xc_t',(nmhu-ns_f(ir))*npn,1);
0382         MXP_f_t = sparse([1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn],MXc,(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);
0383         MXa = reshape(Xa_t',(nmhu-ns_f(ir))*npn,1);
0384         MXP_f_t = MXP_f_t + sparse([1+1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-1],MXa(1+1:length(MXa)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);
0385         MXe = reshape(Xe_t',(nmhu-ns_f(ir))*npn,1);
0386         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-1],[1+1:(nmhu-ns_f(ir))*npn],MXe(1:length(MXe)-1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);       
0387         MXb = reshape(Xb_t',(nmhu-ns_f(ir))*npn,1);
0388         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir)):(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))],MXb(1+(nmhu-ns_f(ir)):length(MXb)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);           
0389         MXd = reshape(Xd_t',(nmhu-ns_f(ir))*npn,1);
0390         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))],[1+(nmhu-ns_f(ir)):(nmhu-ns_f(ir))*npn],MXd(1:length(MXd)-(nmhu-ns_f(ir))),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);           
0391         MXf = reshape(Xf_t',(nmhu-ns_f(ir))*npn,1);
0392         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))+1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-1],MXf(1+(nmhu-ns_f(ir))+1:length(MXf)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0393         MXh = reshape(Xh_t',(nmhu-ns_f(ir))*npn,1);
0394         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))-1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+1],MXh(1+(nmhu-ns_f(ir))-1:length(MXh)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0395         MXi = reshape(Xi_t',(nmhu-ns_f(ir))*npn,1);
0396         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+1],[1+(nmhu-ns_f(ir))-1:(nmhu-ns_f(ir))*npn],MXi(1:length(MXi)-(nmhu-ns_f(ir))+1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0397         MXg = reshape(Xg_t',(nmhu-ns_f(ir))*npn,1);
0398         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-1],[1+(nmhu-ns_f(ir))+1:(nmhu-ns_f(ir))*npn],MXg(1:length(MXg)-(nmhu-ns_f(ir))-1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0399     elseif bounce_mode == 1,
0400         MXc = reshape(Xc_t',(nmhu-ns_f(ir))*npn,1);%c(i+1/2,j+1/2)
0401         MXP_f_t = sparse([1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn],MXc,(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);
0402         MXa = reshape(Xa_t',(nmhu-ns_f(ir))*npn,1);%a(i+1/2,j-1/2)
0403         MXP_f_t = MXP_f_t + sparse([1+1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-1],MXa(1+1:length(MXa)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);
0404         MXe = reshape(Xe_t',(nmhu-ns_f(ir))*npn,1);%e(i+1/2,j+3/2)
0405         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-1],[1+1:(nmhu-ns_f(ir))*npn],MXe(1:length(MXe)-1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);           
0406         MXj = reshape(Xj_t',(nmhu-ns_f(ir))*npn,1);%j(i-1/2,j+1/2-djmhut)
0407         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))+djmhut(ir):(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-djmhut(ir)],MXj(1+(nmhu-ns_f(ir))+djmhut(ir):length(MXj)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);              
0408         MXk = reshape(Xk_t',(nmhu-ns_f(ir))*npn,1);%k(i+1/2,j+1/2-djmhut)
0409         MXP_f_t = MXP_f_t + sparse([1+djmhut(ir):(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-djmhut(ir)],MXk(1+djmhut(ir):length(MXk)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);        
0410         MXl = reshape(Xl_t',(nmhu-ns_f(ir))*npn,1);%l(i+3/2,j+1/2-djmhut)
0411         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+djmhut(ir)],[1+(nmhu-ns_f(ir))-djmhut(ir):(nmhu-ns_f(ir))*npn],MXl(1:length(MXl)-(nmhu-ns_f(ir))+djmhut(ir)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);      
0412         MXf = reshape(Xf_t',(nmhu-ns_f(ir))*npn,1);%f(i-1/2,j-1/2)
0413         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))+1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-1],MXf(1+(nmhu-ns_f(ir))+1:length(MXf)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0414         MXb = reshape(Xb_t',(nmhu-ns_f(ir))*npn,1);%b(i-1/2,j+1/2)
0415         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir)):(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))],MXb(1+(nmhu-ns_f(ir)):length(MXb)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);           
0416         MXh = reshape(Xh_t',(nmhu-ns_f(ir))*npn,1);%h(i-1/2,j+3/2)
0417         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))-1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+1],MXh(1+(nmhu-ns_f(ir))-1:length(MXh)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0418         MXi = reshape(Xi_t',(nmhu-ns_f(ir))*npn,1);%i(i+3/2,j-1/2)
0419         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+1],[1+(nmhu-ns_f(ir))-1:(nmhu-ns_f(ir))*npn],MXi(1:length(MXi)-(nmhu-ns_f(ir))+1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0420         MXd = reshape(Xd_t',(nmhu-ns_f(ir))*npn,1);%d(i+3/2,j+1/2)
0421         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))],[1+(nmhu-ns_f(ir)):(nmhu-ns_f(ir))*npn],MXd(1:length(MXd)-(nmhu-ns_f(ir))),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);                   
0422         MXg = reshape(Xg_t',(nmhu-ns_f(ir))*npn,1);%g(i+3/2,j+3/2)
0423         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-1],[1+(nmhu-ns_f(ir))+1:(nmhu-ns_f(ir))*npn],MXg(1:length(MXg)-(nmhu-ns_f(ir))-1),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);     
0424         MXm = reshape(Xm_t',(nmhu-ns_f(ir))*npn,1);%m(i-1/2,j+1/2+djmhut)
0425         MXP_f_t = MXP_f_t + sparse([1+(nmhu-ns_f(ir))-djmhut(ir):(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))+djmhut(ir)],MXm(1+(nmhu-ns_f(ir))-djmhut(ir):length(MXm)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);      
0426         MXn = reshape(Xn_t',(nmhu-ns_f(ir))*npn,1);%n(i+1/2,j+1/2+djmhut)
0427         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-djmhut(ir)],[1+djmhut(ir):(nmhu-ns_f(ir))*npn],MXn(1:length(MXn)-djmhut(ir)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);        
0428         MXo = reshape(Xo_t',(nmhu-ns_f(ir))*npn,1);%o(i+3/2,j+1/2+djmhut)
0429         MXP_f_t = MXP_f_t + sparse([1:(nmhu-ns_f(ir))*npn-(nmhu-ns_f(ir))-djmhut(ir)],[1+(nmhu-ns_f(ir))+djmhut(ir):(nmhu-ns_f(ir))*npn],MXo(1:length(MXo)-(nmhu-ns_f(ir))-djmhut(ir)),(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);      
0430     end
0431 %
0432     sm_f(ir+1) = sm_f(ir) + blocksize_f_t(ir);
0433     [i0_f,j0_f,s0_f] = find(MMXP_f_t);      
0434     [m0_f,n0_f] = size(MMXP_f_t);
0435     [i1_f,j1_f,s1_f] = find(MXP_f_t);
0436     [m1_f,n1_f] = size(MXP_f_t);
0437     MMXP_f_t = sparse(i0_f,j0_f,s0_f,m0_f+m1_f,m0_f+m1_f) + sparse(m0_f+i1_f,m0_f+j1_f,s1_f,m0_f+m1_f,m0_f+m1_f);
0438 end
0439 %
0440 % Radial supermatrix calculation MMXR_f_t
0441 %
0442 if ~isempty(transpfaste),%boundary conditions
0443     %
0444     % Interpolation coefficients
0445     %
0446     xdpsinmm = [0,xdpsin(1:nr-1)];
0447     xdpsinpp = [xdpsin(2:nr),0];
0448     %
0449     XXdpsin = reshape(repmat(xdpsin,npn*nmhu,1),[npn,nmhu,nr]);   
0450     XXdpsinm = reshape(repmat(xdpsinm,npn*nmhu,1),[npn,nmhu,nr]); 
0451     XXdpsinp = reshape(repmat(xdpsinp,npn*nmhu,1),[npn,nmhu,nr]); 
0452     XXdpsinmm = reshape(repmat(xdpsinmm,npn*nmhu,1),[npn,nmhu,nr]);   
0453     XXdpsinpp = reshape(repmat(xdpsinpp,npn*nmhu,1),[npn,nmhu,nr]);
0454     %
0455     XXdeltalp = XXdpsinpp./(XXdpsinpp + XXdpsin);%Grid interpolation coefficients
0456     XXdeltalm = XXdpsin./(XXdpsinmm + XXdpsin);%Grid interpolation coefficients
0457     %
0458     % Diagonal calculations
0459     %
0460     XXDpr_ipj = ZXXD_c.pr_ipj + ZXXD_e.pr_ipj + ZXXD_s.pr_ipj + ZXXD_rf.pr_ipj + ZXXD_a.pr_ipj;
0461     XXDpr_imj = ZXXD_c.pr_imj + ZXXD_e.pr_imj + ZXXD_s.pr_imj + ZXXD_rf.pr_imj + ZXXD_a.pr_imj;
0462     %
0463     XXDmr_ijp = ZXXD_c.mr_ijp + ZXXD_e.mr_ijp + ZXXD_s.mr_ijp + ZXXD_rf.mr_ijp + ZXXD_a.mr_ijp;
0464     XXDmr_ijm = ZXXD_c.mr_ijm + ZXXD_e.mr_ijm + ZXXD_s.mr_ijm + ZXXD_rf.mr_ijm + ZXXD_a.mr_ijm;
0465     %
0466     XXDrr_lp = ZXXD_c.rr_lp + ZXXD_e.rr_lp + ZXXD_s.rr_lp + ZXXD_rf.rr_lp + ZXXD_a.rr_lp;
0467     XXDrr_lm = ZXXD_c.rr_lm + ZXXD_e.rr_lm + ZXXD_s.rr_lm + ZXXD_rf.rr_lm + ZXXD_a.rr_lm;
0468     XXDrp_lp = ZXXD_c.rp_lp + ZXXD_e.rp_lp + ZXXD_s.rp_lp + ZXXD_rf.rp_lp + ZXXD_a.rp_lp;
0469     XXDrp_lm = ZXXD_c.rp_lm + ZXXD_e.rp_lm + ZXXD_s.rp_lm + ZXXD_rf.rp_lm + ZXXD_a.rp_lm;
0470     XXDrm_lp = ZXXD_c.rm_lp + ZXXD_e.rm_lp + ZXXD_s.rm_lp + ZXXD_rf.rm_lp + ZXXD_a.rm_lp;
0471     XXDrm_lm = ZXXD_c.rm_lm + ZXXD_e.rm_lm + ZXXD_s.rm_lm + ZXXD_rf.rm_lm + ZXXD_a.rm_lm;
0472     XXFr_lp = ZXXF_c.r_lp + ZXXF_e.r_lp + ZXXF_s.r_lp + ZXXF_rf.r_lp + ZXXF_a.r_lp;
0473     XXFr_lm = ZXXF_c.r_lm + ZXXF_e.r_lm + ZXXF_s.r_lm + ZXXF_rf.r_lm + ZXXF_a.r_lm;
0474     %
0475     XXx0 = repmat(reshape(xx0',[1,1,nr]),[npn,nmhu,1]);
0476     XXx0mm = zeros(npn,nmhu,nr);
0477     XXx0mm(:,:,2:nr) = XXx0(:,:,1:nr-1);
0478     XXx0pp = zeros(npn,nmhu,nr);
0479     XXx0pp(:,:,1:nr-1) = XXx0(:,:,2:nr);
0480     %
0481     XXx0p = (1 - XXdeltalp).*XXx0pp + XXdeltalp.*XXx0;%Evaluation on the flux grid (linear interpolation)
0482     XXx0m = (1 - XXdeltalm).*XXx0 + XXdeltalm.*XXx0mm;%Evaluation on the flux grid (linear interpolation)
0483     %
0484     XXqtilde = repmat(reshape(xqtilde',[1,1,nr]),[npn,nmhu,1]);
0485     XXqtildemm = zeros(npn,nmhu,nr);
0486     XXqtildemm(:,:,2:nr) = XXqtilde(:,:,1:nr-1);
0487     XXqtildepp = zeros(npn,nmhu,nr);
0488     XXqtildepp(:,:,1:nr-1) = XXqtilde(:,:,2:nr);
0489     %
0490     XXqtildep = (1 - XXdeltalp).*XXqtildepp + XXdeltalp.*XXqtilde;%Evaluation on the flux grid (linear interpolation)
0491     XXqtildem = (1 - XXdeltalm).*XXqtilde + XXdeltalm.*XXqtildemm;%Evaluation on the flux grid (linear interpolation)
0492     %
0493     XXBp0 = repmat(reshape(xBp0',[1,1,nr]),[npn,nmhu,1]);
0494     XXBp0mm = zeros(npn,nmhu,nr);
0495     XXBp0mm(:,:,2:nr) = XXBp0(:,:,1:nr-1);
0496     XXBp0pp = zeros(npn,nmhu,nr);
0497     XXBp0pp(:,:,1:nr-1) = XXBp0(:,:,2:nr);
0498     %
0499     XXBp0p = (1 - XXdeltalp).*XXBp0pp + XXdeltalp.*XXBp0;%Evaluation on the flux grid (linear interpolation)
0500     XXBp0m = (1 - XXdeltalm).*XXBp0 + XXdeltalm.*XXBp0mm;%Evaluation on the flux grid (linear interpolation)
0501     %
0502     XXB0 = repmat(reshape(xB0',[1,1,nr]),[npn,nmhu,1]);
0503     XXB0mm = zeros(npn,nmhu,nr);
0504     XXB0mm(:,:,2:nr) = XXB0(:,:,1:nr-1);
0505     XXB0pp = zeros(npn,nmhu,nr);
0506     XXB0pp(:,:,1:nr-1) = XXB0(:,:,2:nr);
0507     %
0508     XXB0p = (1 - XXdeltalp).*XXB0pp + XXdeltalp.*XXB0;%Evaluation on the flux grid (linear interpolation)
0509     XXB0m = (1 - XXdeltalm).*XXB0 + XXdeltalm.*XXB0mm;%Evaluation on the flux grid (linear interpolation)
0510     %
0511     XXlambdamm = zeros(npn,nmhu,nr);
0512     XXlambdamm(:,:,2:nr) = XXlambda(:,:,1:nr-1);
0513     XXlambdapp = zeros(npn,nmhu,nr);
0514     XXlambdapp(:,:,1:nr-1) = XXlambda(:,:,2:nr);
0515     %
0516     XXlambdap = (1 - XXdeltalp).*XXlambdapp + XXdeltalp.*XXlambda;%Evaluation on the flux grid (linear interpolation)
0517     XXlambdam = (1 - XXdeltalm).*XXlambda + XXdeltalm.*XXlambdamm;%Evaluation on the flux grid (linear interpolation)
0518     %
0519     XXgradpsinp = (1.0 + XXx0p/Rp).*XXBp0p*ap^2/psia_apRp;
0520     XXgradpsinm = (1.0 + XXx0m/Rp).*XXBp0m*ap^2/psia_apRp;
0521     %
0522     XXalpharp = XXqtildep.*XXgradpsinp.*XXlambdap./XXB0p;%Evaluation on the flux grid (linear interpolation)
0523     XXalpharm = XXqtildem.*XXgradpsinm.*XXlambdam./XXB0m;%Evaluation on the flux grid (linear interpolation)
0524     XXalpharm(:,:,1) = 0;%Force zero at the plasma center, which is the natural analytical limit (In cylindrical geometry, r=0)
0525     XXalpharp(:,:,nr) = 0;%Avoid NaN
0526     %
0527     MMXR_f_t = sparse(nMMXP_f_t,nMMXP_f_t);%Create supermatrix
0528     %
0529     XXMrm = -(XXalpharm.*XXgradpsinm./XXdpsinm./XXdpsin).*XXDrr_lm - (XXalpharm./XXdpsin).*XXFr_lm.*XXdeltalm;
0530     XXMr = (XXalpharp.*XXgradpsinp./XXdpsinp./XXdpsin).*XXDrr_lp + (XXalpharp./XXdpsin).*XXFr_lp.*XXdeltalp + (XXalpharm.*XXgradpsinm./XXdpsinm./XXdpsin).*XXDrr_lm - (XXalpharm./XXdpsin).*XXFr_lm.*(1 - XXdeltalm);                    
0531     XXMrp = -(XXalpharp.*XXgradpsinp./XXdpsinp./XXdpsin).*XXDrr_lp + (XXalpharp./XXdpsin).*XXFr_lp.*(1 - XXdeltalp);
0532     %
0533     XXMrm = (XXB0./XXqtilde).*XXMrm.*XXpn2./XXlambda;%Jacobian multiplication
0534     XXMr = (XXB0./XXqtilde).*XXMr.*XXpn2./XXlambda;%Jacobian multiplication
0535     XXMrp = (XXB0./XXqtilde).*XXMrp.*XXpn2./XXlambda;%Jacobian multiplication
0536     %
0537     XXMrm(:,:,nr) = 0;%Edge boundary condition (Maxwellian solution)
0538     XXMr(:,:,nr) = 0;%Edge boundary condition (Maxwellian solution)
0539 %
0540 % Radial off-diagonal terms for fast electron transport
0541 %
0542     if bounce_mode == 0,
0543         XXMrm = permute(XXMrm,[2,1,3]);
0544         MMrm_t = reshape(XXMrm,npn*nmhu*nr,1,1);
0545         XXMr = permute(XXMr,[2,1,3]);
0546         MMr_t = reshape(XXMr,npn*nmhu*nr,1,1);
0547         XXMrp = permute(XXMrp,[2,1,3]);
0548         MMrp_t = reshape(XXMrp,npn*nmhu*nr,1,1);
0549         %
0550         MMXR_f_t = MMXR_f_t + sparse(1+nmhu*npn:nMMXP_f_t,1:nMMXP_f_t-nmhu*npn,MMrm_t(1+npn*nmhu:nMMXP_f_t),nMMXP_f_t,nMMXP_f_t);
0551         MMXR_f_t = MMXR_f_t + sparse(1:nMMXP_f_t,1:nMMXP_f_t,MMr_t,nMMXP_f_t,nMMXP_f_t);    
0552         MMXR_f_t = MMXR_f_t + sparse(1:nMMXP_f_t-nmhu*npn,1+npn*nmhu:nMMXP_f_t,MMrp_t(1:nMMXP_f_t-nmhu*npn),nMMXP_f_t,nMMXP_f_t);
0553         %
0554         if display_mode >= 1,info_dke_yp(2,['Zero order configuration space matrix calculations done in ',num2str(etime(clock,time0)),' (s)']);end
0555     elseif bounce_mode == 1,
0556         %
0557         % First magnetic flux surface
0558         %
0559         ir = 1;
0560         ib0 = 0;
0561         ibp = max(cumsum(blocksize_f_t(1:ir)));
0562         %
0563         Xmask_f0 = (Xmhu2 < xmhubounce2(ir)) & (Xmhu < 0);%Only the half-bounce region is removed
0564         Xmask_f0 = nonzeros(~(Xmask_f0(2,:)==1).*cumsum(ones(1,nmhu)))';
0565         Xmask_f0p = (Xmhu2 < xmhubounce2(ir+1)) & (Xmhu < 0);%Only the half-bounce region is removed
0566         Xmask_f0p = nonzeros(~(Xmask_f0p(2,:)==1).*cumsum(ones(1,nmhu)))';
0567         %
0568         imhu0 = sum((Xmhu(1,Xmask_f0) < 0).*(Xmhu2(1,Xmask_f0) > xmhubounce2(ir))); 
0569         imhup = sum((Xmhu(1,Xmask_f0p) < 0).*(Xmhu2(1,Xmask_f0p) > xmhubounce2(ir+1)));
0570         %
0571         XMr = XXMr(:,:,ir);
0572         XMr_t = XMr(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0573         MMr_t = reshape(XMr_t',npn*(nmhu-ns_f(ir)),1,1);
0574         %
0575         MMXR_f_t_i = zeros(nMMXP_f_t*11*nr,1);%Create oversized sparse matrix to avoid reallocation
0576         MMXR_f_t_j = zeros(nMMXP_f_t*11*nr,1);
0577         MMXR_f_t_v = zeros(nMMXP_f_t*11*nr,1);
0578         %
0579         i = 1+ib0:ibp;
0580         j = 1+ib0:ibp;
0581         ia = 1;ib = ia+length(i);                 
0582         MMXR_f_t_i(ia:ib-1) = i;
0583         MMXR_f_t_j(ia:ib-1) = j;
0584         MMXR_f_t_v(ia:ib-1) = MMr_t;
0585         %
0586         % Notations: bc: backward + circulating, t: trapped (half positive part only), fc: forward + circulating
0587         %
0588         XMrp = XXMrp(:,:,ir);
0589         XMrp_t = XMrp(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0590         XMrp1_t = XMrp_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir+1));%bc -> bc
0591         XMrp2_t = XMrp_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir+1)).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));%bc -> t
0592         XMrp3_t = XMrp_t.*(Xmhu(:,Xmask_f0) > 0);%fc + t -> fc + t
0593         %
0594         for ip = 1:npn, 
0595             i = ib0+(nmhu-ns_f(ir))*(ip-1)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhup;
0596             j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+1:ibp+(nmhu-ns_f(ir+1))*(ip-1)+imhup;
0597             ia = ib;ib = ia+length(i);
0598             MMXR_f_t_i(ia:ib-1) = i;
0599             MMXR_f_t_j(ia:ib-1) = j;
0600             MMXR_f_t_v(ia:ib-1) = XMrp1_t(ip,1:imhup)';%bc -> bc
0601             %
0602             i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhup+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0;
0603             j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+((nmhu-ns_f(ir+1))-imhup):-1:ibp+(nmhu-ns_f(ir+1))*(ip-1)+((nmhu-ns_f(ir+1))-imhu0)+1;
0604             ia = ib;ib = ia+length(i);
0605             MMXR_f_t_i(ia:ib-1) = i;
0606             MMXR_f_t_j(ia:ib-1) = j;
0607             MMXR_f_t_v(ia:ib-1) = XMrp2_t(ip,imhup+1:imhu0)';%bc -> t
0608             %
0609             i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0+1:ib0+(nmhu-ns_f(ir))*ip;
0610             j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+imhup+1:ibp+(nmhu-ns_f(ir+1))*ip;
0611             ia = ib;ib = ia+length(i);
0612             MMXR_f_t_i(ia:ib-1) = i;
0613             MMXR_f_t_j(ia:ib-1) = j;
0614             MMXR_f_t_v(ia:ib-1) = XMrp3_t(ip,imhu0+1:nmhu-ns_f(ir))';%fc + t -> fc + t
0615         end 
0616         %
0617         % All magnetic flux surface except the first and the last ones
0618         %
0619         for ir = 2:nr-1,
0620             if ir == 2,
0621                 ibm = 0;
0622             else
0623                 ibm = max(cumsum(blocksize_f_t(1:ir-2)));%Block index reference
0624             end
0625             %
0626             ib0 = max(cumsum(blocksize_f_t(1:ir-1)));%Block index reference
0627             ibp = max(cumsum(blocksize_f_t(1:ir)));
0628             %
0629             Xmask_f0m = (Xmhu2 < xmhubounce2(ir-1)) & (Xmhu < 0);%Only the half-bounce region is removed
0630             Xmask_f0m = nonzeros(~(Xmask_f0m(2,:)==1).*cumsum(ones(1,nmhu)))';
0631             Xmask_f0 = (Xmhu2 < xmhubounce2(ir)) & (Xmhu < 0);%Only the half-bounce region is removed
0632             Xmask_f0 = nonzeros(~(Xmask_f0(2,:)==1).*cumsum(ones(1,nmhu)))';
0633             Xmask_f0p = (Xmhu2 < xmhubounce2(ir+1)) & (Xmhu < 0);%Only the half-bounce region is removed
0634             Xmask_f0p = nonzeros(~(Xmask_f0p(2,:)==1).*cumsum(ones(1,nmhu)))';
0635             %
0636             imhum = sum((Xmhu(1,Xmask_f0m) < 0).*(Xmhu2(1,Xmask_f0m) > xmhubounce2(ir-1)));
0637             imhu0 = sum((Xmhu(1,Xmask_f0) < 0).*(Xmhu2(1,Xmask_f0) > xmhubounce2(ir))); 
0638             imhup = sum((Xmhu(1,Xmask_f0p) < 0).*(Xmhu2(1,Xmask_f0p) > xmhubounce2(ir+1)));
0639             %
0640             XMr = XXMr(:,:,ir);
0641             XMr_t = XMr(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0642             MMr_t = reshape(XMr_t',npn*(nmhu-ns_f(ir)),1,1);
0643             %
0644             i = 1+ib0:ibp;
0645             j = 1+ib0:ibp;
0646             ia = ib;ib = ia+length(i); 
0647             MMXR_f_t_i(ia:ib-1) = i;
0648             MMXR_f_t_j(ia:ib-1) = j;
0649             MMXR_f_t_v(ia:ib-1) = MMr_t;
0650             %
0651             % Notations: bc: backward + circulating, t: trapped (half positive part only), fc: forward + circulating
0652             %
0653             XMrp = XXMrp(:,:,ir);
0654             XMrp_t = XMrp(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0655             XMrp1_t = XMrp_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir+1));%bc -> bc
0656             XMrp2_t = XMrp_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir+1)).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));%bc -> t
0657             XMrp3_t = XMrp_t.*(Xmhu(:,Xmask_f0) > 0);%t+fc -> t+fc
0658             %
0659             XMrm = XXMrm(:,:,ir);
0660             XMrm_t = XMrm(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0661             XMrm1_t = XMrm_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));%bc -> bc
0662             XMrm2_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir-1));%t -> t
0663             XMrm3_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir-1)).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir));%t -> fc
0664             XMrm4_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));%fc -> fc
0665             %
0666             for ip = 1:npn, 
0667                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhup;
0668                 j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+1:ibp+(nmhu-ns_f(ir+1))*(ip-1)+imhup;
0669                 ia = ib;ib = ia+length(i);     
0670                 MMXR_f_t_i(ia:ib-1) = i;
0671                 MMXR_f_t_j(ia:ib-1) = j;
0672                 MMXR_f_t_v(ia:ib-1) = XMrp1_t(ip,1:imhup)';%bc -> bc
0673                 %
0674                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhup+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0;
0675                 j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+((nmhu-ns_f(ir+1))-imhup):-1:ibp+(nmhu-ns_f(ir+1))*(ip-1)+((nmhu-ns_f(ir+1))-imhu0)+1;
0676                 ia = ib;ib = ia+length(i);
0677                 MMXR_f_t_i(ia:ib-1) = i;
0678                 MMXR_f_t_j(ia:ib-1) = j;
0679                 MMXR_f_t_v(ia:ib-1) = XMrp2_t(ip,imhup+1:imhu0)';  
0680                 %
0681                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0+1:ib0+(nmhu-ns_f(ir))*ip;
0682                 j = ibp+(nmhu-ns_f(ir+1))*(ip-1)+imhup+1:ibp+(nmhu-ns_f(ir+1))*ip;
0683                 ia = ib;ib = ia+length(i);
0684                 MMXR_f_t_i(ia:ib-1) = i;
0685                 MMXR_f_t_j(ia:ib-1) = j;
0686                 MMXR_f_t_v(ia:ib-1) = XMrp3_t(ip,imhu0+1:nmhu-ns_f(ir))'; 
0687                 %
0688                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0;
0689                 j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhu0;
0690                 ia = ib;ib = ia+length(i);
0691                 MMXR_f_t_i(ia:ib-1) = i;
0692                 MMXR_f_t_j(ia:ib-1) = j;
0693                 MMXR_f_t_v(ia:ib-1) = XMrm1_t(ip,1:imhu0)'; 
0694                 %
0695                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0+1:ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhum);
0696                 j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhum+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhum);
0697                 ia = ib;ib = ia+length(i);
0698                 MMXR_f_t_i(ia:ib-1) = i;
0699                 MMXR_f_t_j(ia:ib-1) = j;
0700                 MMXR_f_t_v(ia:ib-1) = XMrm2_t(ip,imhu0+1:((nmhu-ns_f(ir))-imhum))'; 
0701                 %
0702                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhum)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhu0);
0703                 j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhum)+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhu0);               
0704                 ia = ib;ib = ia+length(i);
0705                 MMXR_f_t_i(ia:ib-1) = i;
0706                 MMXR_f_t_j(ia:ib-1) = j;
0707                 MMXR_f_t_v(ia:ib-1) = XMrm3_t(ip,((nmhu-ns_f(ir))-imhum)+1:((nmhu-ns_f(ir))-imhu0))'/2; 
0708                 j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhum:-1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhu0+1;
0709                 ia = ib;ib = ia+length(i);
0710                 MMXR_f_t_i(ia:ib-1) = i;
0711                 MMXR_f_t_j(ia:ib-1) = j;
0712                 MMXR_f_t_v(ia:ib-1) = XMrm3_t(ip,((nmhu-ns_f(ir))-imhum)+1:((nmhu-ns_f(ir))-imhu0))'/2; 
0713                 %
0714                 i = ib0+(nmhu-ns_f(ir))*(ip-1)+(nmhu-ns_f(ir))-imhu0+1:ib0+(nmhu-ns_f(ir))*ip;
0715                 j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+(nmhu-ns_f(ir-1))-imhu0+1:ibm+(nmhu-ns_f(ir-1))*ip;
0716                 ia = ib;ib = ia+length(i);
0717                 MMXR_f_t_i(ia:ib-1) = i;
0718                 MMXR_f_t_j(ia:ib-1) = j;
0719                 MMXR_f_t_v(ia:ib-1) = XMrm4_t(ip,(nmhu-ns_f(ir))-imhu0+1:nmhu-ns_f(ir))';
0720             end 
0721         end
0722         %
0723         % Last magnetic flux surface
0724         %
0725         ir = nr;
0726         ibm = max(cumsum(blocksize_f_t(1:ir-2)));%Block index reference
0727         ib0 = max(cumsum(blocksize_f_t(1:ir-1)));%Block index reference
0728         ibp = max(cumsum(blocksize_f_t(1:ir)));
0729         %
0730         Xmask_f0m = (Xmhu2 < xmhubounce2(ir-1)) & (Xmhu < 0);%Only the half-bounce region is removed
0731         Xmask_f0m = nonzeros(~(Xmask_f0m(2,:)==1).*cumsum(ones(1,nmhu)))';
0732         Xmask_f0 = (Xmhu2 < xmhubounce2(ir)) & (Xmhu < 0);%Only the half-bounce region is removed
0733         Xmask_f0 = nonzeros(~(Xmask_f0(2,:)==1).*cumsum(ones(1,nmhu)))';
0734         %
0735         imhum = sum((Xmhu(1,Xmask_f0m) < 0).*(Xmhu2(1,Xmask_f0m) > xmhubounce2(ir-1)));
0736         imhu0 = sum((Xmhu(1,Xmask_f0) < 0).*(Xmhu2(1,Xmask_f0) > xmhubounce2(ir))); 
0737         %
0738         XMr = XXMr(:,:,ir);
0739         XMr_t = XMr(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0740         MMr_t = reshape(XMr_t',npn*(nmhu-ns_f(ir)),1,1);
0741         %
0742         i = 1+ib0:ibp;
0743         j = 1+ib0:ibp;
0744         ia = ib;ib = ia+length(i);
0745         MMXR_f_t_i(ia:ib-1) = i;
0746         MMXR_f_t_j(ia:ib-1) = j;
0747         MMXR_f_t_v(ia:ib-1) = MMr_t;
0748         %
0749         XMrm = XXMrm(:,:,ir);
0750         XMrm_t = XMrm(:,Xmask_f0);%Matrix reduction along the angular direction (half-bounce region is removed)
0751         XMrm1_t = XMrm_t.*(Xmhu(:,Xmask_f0) < 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));
0752         %         XMrm2_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0);
0753         XMrm2_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir-1));%t -> t
0754         XMrm3_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir-1)).*(Xmhu2(:,Xmask_f0) < xmhubounce2(ir));%t -> fc
0755         XMrm4_t = XMrm_t.*(Xmhu(:,Xmask_f0) > 0).*(Xmhu2(:,Xmask_f0) > xmhubounce2(ir));%fc -> fc
0756 %
0757         for ip = 1:npn, 
0758             i = ib0+(nmhu-ns_f(ir))*(ip-1)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0;
0759             j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhu0;
0760             ia = ib;ib = ia+length(i);
0761             MMXR_f_t_i(ia:ib-1) = i;
0762             MMXR_f_t_j(ia:ib-1) = j;
0763             MMXR_f_t_v(ia:ib-1) = XMrm1_t(ip,1:imhu0)'; 
0764             %
0765             i = ib0+(nmhu-ns_f(ir))*(ip-1)+imhu0+1:ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhum);
0766             j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhum+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhum);
0767             ia = ib;ib = ia+length(i);
0768             MMXR_f_t_i(ia:ib-1) = i;
0769             MMXR_f_t_j(ia:ib-1) = j;
0770             MMXR_f_t_v(ia:ib-1) = XMrm2_t(ip,imhu0+1:((nmhu-ns_f(ir))-imhum))'; 
0771             %
0772             i = ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhum)+1:ib0+(nmhu-ns_f(ir))*(ip-1)+((nmhu-ns_f(ir))-imhu0);
0773             j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhum)+1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+((nmhu-ns_f(ir-1))-imhu0);
0774             ia = ib;ib = ia+length(i);
0775             MMXR_f_t_i(ia:ib-1) = i;
0776             MMXR_f_t_j(ia:ib-1) = j;
0777             MMXR_f_t_v(ia:ib-1) = XMrm3_t(ip,((nmhu-ns_f(ir))-imhum)+1:((nmhu-ns_f(ir))-imhu0))'/2;
0778             j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhum:-1:ibm+(nmhu-ns_f(ir-1))*(ip-1)+imhu0+1;
0779             ia = ib;ib = ia+length(i);
0780             MMXR_f_t_i(ia:ib-1) = i;
0781             MMXR_f_t_j(ia:ib-1) = j;
0782             MMXR_f_t_v(ia:ib-1) = XMrm3_t(ip,((nmhu-ns_f(ir))-imhum)+1:((nmhu-ns_f(ir))-imhu0))'/2; 
0783             %
0784             i = ib0+(nmhu-ns_f(ir))*(ip-1)+(nmhu-ns_f(ir))-imhu0+1:ib0+(nmhu-ns_f(ir))*ip;
0785             j = ibm+(nmhu-ns_f(ir-1))*(ip-1)+(nmhu-ns_f(ir-1))-imhu0+1:ibm+(nmhu-ns_f(ir-1))*ip;
0786             ia = ib;ib = ia+length(i);
0787             MMXR_f_t_i(ia:ib-1) = i;
0788             MMXR_f_t_j(ia:ib-1) = j;
0789             MMXR_f_t_v(ia:ib-1) = XMrm4_t(ip,(nmhu-ns_f(ir))-imhu0+1:nmhu-ns_f(ir))';
0790         end
0791         %
0792         MMXR_f_t = sparse(MMXR_f_t_i(1:ib-1),MMXR_f_t_j(1:ib-1),MMXR_f_t_v(1:ib-1),nMMXP_f_t,nMMXP_f_t);%Create sparse radial transport matrix
0793         %
0794         if display_mode >= 1,info_dke_yp(2,['Zero order configuration space matrix calculations done in ',num2str(etime(clock,time0)),' (s)']);end
0795     end     
0796 else
0797     MMXR_f_t = sparse(nMMXP_f_t,nMMXP_f_t);%Create empty supermatrix
0798 end
0799 %
0800 % Time diagonal supermatrix calculation MMXT_f_t
0801 %
0802 for ir = 1:nr,
0803     MXT = reshape(Zpn2_t{ir}',(nmhu-ns_f(ir))*npn,1)/dtn;%Partial momentum Jacobian
0804     MXT_f_t = sparse([1:(nmhu-ns_f(ir))*npn],[1:(nmhu-ns_f(ir))*npn],MXT,(nmhu-ns_f(ir))*npn,(nmhu-ns_f(ir))*npn);%Resizing for 3D calculations + multiplication by the remaining Jacobian
0805 %
0806     sm_T_f(ir+1) = sm_T_f(ir) + blocksize_T_f_t(ir);
0807     [i0_T_f,j0_T_f,s0_T_f] = find(MMXT_f_t);        
0808     [m0_T_f,n0_T_f] = size(MMXT_f_t);
0809     [i1_T_f,j1_T_f,s1_T_f] = find(MXT_f_t);
0810     [m1_T_f,n1_T_f] = size(MXT_f_t);
0811     MMXT_f_t = sparse(i0_T_f,j0_T_f,s0_T_f,m0_T_f+m1_T_f,m0_T_f+m1_T_f) + sparse(m0_T_f+i1_T_f,m0_T_f+j1_T_f,s1_T_f,m0_T_f+m1_T_f,m0_T_f+m1_T_f);   
0812 %
0813     etime_fp = etime_fp + etime(clock,time0);
0814     time0 = clock;
0815 end
0816 %
0817 Zdelta.XXdeltap_imj = XXdeltap_imj; 
0818 Zdelta.XXdeltap_ipj = XXdeltap_ipj;
0819 Zdelta.XXdeltam_ijp = XXdeltam_ijp;
0820 Zdelta.XXdeltam_ijm = XXdeltam_ijm;   
0821 Zdelta.XXdeltap_imjmm = XXdeltap_imjmm;
0822 Zdelta.XXdeltap_imjpp = XXdeltap_imjpp;
0823 Zdelta.XXdeltap_ipjmm = XXdeltap_ipjmm;
0824 Zdelta.XXdeltap_ipjpp = XXdeltap_ipjpp;
0825 Zdelta.XXdeltam_immjm = XXdeltam_immjm;
0826 Zdelta.XXdeltam_ippjm = XXdeltam_ippjm;
0827 Zdelta.XXdeltam_immjp = XXdeltam_immjp;
0828 Zdelta.XXdeltam_ippjp = XXdeltam_ippjp;
0829 %
0830 if display_mode >= 1,info_dke_yp(2,['Zero order momentum space matrix calculations done in ',num2str(etime_fp),' (s)']);end
0831 %
0832 % DKE (first order) matrix calculation
0833 %
0834 etime_dke = 0;
0835 time0 = clock;
0836 %
0837 if dke_mode == 1,
0838     for ir_dke = 1:nr_dke,
0839 %
0840         jr1 = rdke(ir_dke)-1;
0841         jr2 = rdke(ir_dke);
0842         jr3 = rdke(ir_dke)+1;
0843 %
0844 % D and F tensors
0845 %
0846         XDpp_tp_ippj = ZXXD_c_tp.pp_ippj(:,:,ir_dke) + ZXXD_e_tp.pp_ippj(:,:,ir_dke) + ZXXD_s_tp.pp_ippj(:,:,ir_dke) + ZXXD_rf_tp.pp_ippj(:,:,ir_dke);
0847         XDpp_tp_ipj = ZXXD_c_tp.pp_ipj(:,:,ir_dke) + ZXXD_e_tp.pp_ipj(:,:,ir_dke) + ZXXD_s_tp.pp_ipj(:,:,ir_dke) + ZXXD_rf_tp.pp_ipj(:,:,ir_dke);
0848         XDpp_tp_ij = ZXXD_c_tp.pp_ij(:,:,ir_dke) + ZXXD_e_tp.pp_ij(:,:,ir_dke) + ZXXD_s_tp.pp_ij(:,:,ir_dke) + ZXXD_rf_tp.pp_ij(:,:,ir_dke);
0849         XDpp_tp_imj = ZXXD_c_tp.pp_imj(:,:,ir_dke) + ZXXD_e_tp.pp_imj(:,:,ir_dke) + ZXXD_s_tp.pp_imj(:,:,ir_dke) + ZXXD_rf_tp.pp_imj(:,:,ir_dke);
0850         XDpp_tp_immj = ZXXD_c_tp.pp_immj(:,:,ir_dke) + ZXXD_e_tp.pp_immj(:,:,ir_dke) + ZXXD_s_tp.pp_immj(:,:,ir_dke) + ZXXD_rf_tp.pp_immj(:,:,ir_dke);
0851 %
0852         XDpm_tp_ipj = ZXXD_c_tp.pm_ipj(:,:,ir_dke) + ZXXD_e_tp.pm_ipj(:,:,ir_dke) + ZXXD_s_tp.pm_ipj(:,:,ir_dke) + ZXXD_rf_tp.pm_ipj(:,:,ir_dke);
0853         XDpm_tp_ij = ZXXD_c_tp.pm_ij(:,:,ir_dke) + ZXXD_e_tp.pm_ij(:,:,ir_dke) + ZXXD_s_tp.pm_ij(:,:,ir_dke) + ZXXD_rf_tp.pm_ij(:,:,ir_dke);
0854         XDpm_tp_imj = ZXXD_c_tp.pm_imj(:,:,ir_dke) + ZXXD_e_tp.pm_imj(:,:,ir_dke) + ZXXD_s_tp.pm_imj(:,:,ir_dke) + ZXXD_rf_tp.pm_imj(:,:,ir_dke);
0855 %
0856         XDmp_tp_ijp = ZXXD_c_tp.mp_ijp(:,:,ir_dke) + ZXXD_e_tp.mp_ijp(:,:,ir_dke) + ZXXD_s_tp.mp_ijp(:,:,ir_dke) + ZXXD_rf_tp.mp_ijp(:,:,ir_dke);
0857         XDmp_tp_ij = ZXXD_c_tp.mp_ij(:,:,ir_dke) + ZXXD_e_tp.mp_ij(:,:,ir_dke) + ZXXD_s_tp.mp_ij(:,:,ir_dke) + ZXXD_rf_tp.mp_ij(:,:,ir_dke);
0858         XDmp_tp_ijm = ZXXD_c_tp.mp_ijm(:,:,ir_dke) + ZXXD_e_tp.mp_ijm(:,:,ir_dke) + ZXXD_s_tp.mp_ijm(:,:,ir_dke) + ZXXD_rf_tp.mp_ijm(:,:,ir_dke);
0859 %
0860         XDmm_tp_ijp = ZXXD_c_tp.mm_ijp(:,:,ir_dke) + ZXXD_e_tp.mm_ijp(:,:,ir_dke) + ZXXD_s_tp.mm_ijp(:,:,ir_dke) + ZXXD_rf_tp.mm_ijp(:,:,ir_dke);
0861         XDmm_tp_ijm = ZXXD_c_tp.mm_ijm(:,:,ir_dke) + ZXXD_e_tp.mm_ijm(:,:,ir_dke) + ZXXD_s_tp.mm_ijm(:,:,ir_dke) + ZXXD_rf_tp.mm_ijm(:,:,ir_dke);
0862 %
0863         XFp_tp_ippj = ZXXF_c_tp.p_ippj(:,:,ir_dke) + ZXXF_e_tp.p_ippj(:,:,ir_dke) + ZXXF_s_tp.p_ippj(:,:,ir_dke) + ZXXF_rf_tp.p_ippj(:,:,ir_dke);
0864         XFp_tp_ipj = ZXXF_c_tp.p_ipj(:,:,ir_dke) + ZXXF_e_tp.p_ipj(:,:,ir_dke) + ZXXF_s_tp.p_ipj(:,:,ir_dke) + ZXXF_rf_tp.p_ipj(:,:,ir_dke);
0865         XFp_tp_ij = ZXXF_c_tp.p_ij(:,:,ir_dke) + ZXXF_e_tp.p_ij(:,:,ir_dke) + ZXXF_s_tp.p_ij(:,:,ir_dke) + ZXXF_rf_tp.p_ij(:,:,ir_dke);
0866         XFp_tp_imj = ZXXF_c_tp.p_imj(:,:,ir_dke) + ZXXF_e_tp.p_imj(:,:,ir_dke) + ZXXF_s_tp.p_imj(:,:,ir_dke) + ZXXF_rf_tp.p_imj(:,:,ir_dke);
0867         XFp_tp_immj = ZXXF_c_tp.p_immj(:,:,ir_dke) + ZXXF_e_tp.p_immj(:,:,ir_dke) + ZXXF_s_tp.p_immj(:,:,ir_dke) + ZXXF_rf_tp.p_immj(:,:,ir_dke);
0868 %
0869         XFm_tp_ijp = ZXXF_c_tp.m_ijp(:,:,ir_dke) + ZXXF_e_tp.m_ijp(:,:,ir_dke) + ZXXF_s_tp.m_ijp(:,:,ir_dke) + ZXXF_rf_tp.m_ijp(:,:,ir_dke);
0870         XFm_tp_ijm = ZXXF_c_tp.m_ijm(:,:,ir_dke) + ZXXF_e_tp.m_ijm(:,:,ir_dke) + ZXXF_s_tp.m_ijm(:,:,ir_dke) + ZXXF_rf_tp.m_ijm(:,:,ir_dke);
0871 %
0872         [Xa_tp,Xb_tp,Xc_tp,Xd_tp,Xe_tp,Xf_tp,Xg_tp,Xh_tp,Xi_tp] = ...
0873             fpengine_dke_yp(dkeparam,1,...
0874             XXdeltap_imj(:,:,jr2),XXdeltap_ipj(:,:,jr2),XXdeltam_ijm(:,:,jr2),XXdeltam_ijp(:,:,jr2),...
0875             XXdeltap_imjmm(:,:,jr2),XXdeltap_imjpp(:,:,jr2),XXdeltap_ipjmm(:,:,jr2),XXdeltap_ipjpp(:,:,jr2),...
0876             XXdeltam_immjm(:,:,jr2),XXdeltam_ippjm(:,:,jr2),XXdeltam_immjp(:,:,jr2),XXdeltam_ippjp(:,:,jr2),...
0877             XXRlambda_m(:,:,jr2),XXRlambda_p(:,:,jr2),...
0878             Xpnmm,Xpnm,Xpn,Xpnp,Xpnpp,Xpn2m,Xpn2p,Xdpnmm,Xdpnm,Xdpn,Xdpnp,Xdpnpp,...
0879             Xmhumm,Xmhum,Xmhu,Xmhup,Xmhupp,X1mmhu2,X1mmhu2m,X1mmhu2p,Xdmhumm,Xdmhum,Xdmhu,Xdmhup,Xdmhupp,...
0880             XDpp_tp_ippj,XDpp_tp_ipj,XDpp_tp_ij,XDpp_tp_imj,XDpp_tp_immj,...
0881             XDpm_tp_ipj,XDpm_tp_ij,XDpm_tp_imj,...
0882             XDmp_tp_ijp,XDmp_tp_ij,XDmp_tp_ijm,...
0883             XDmm_tp_ijp,XDmm_tp_ijm,...
0884             XFp_tp_ippj,XFp_tp_ipj,XFp_tp_ij,XFp_tp_imj,XFp_tp_immj,...
0885             XFm_tp_ijp,XFm_tp_ijm); 
0886 %
0887 %Boundary conditions because of cross-derivatives (symmetry of f: fj-1/2=fj+1/2)
0888 %
0889         Xc_tp(:,1) = Xc_tp(:,1) + Xa_tp(:,1);
0890         Xd_tp(:,1) = Xd_tp(:,1) + Xi_tp(:,1);
0891         Xb_tp(:,1) = Xb_tp(:,1) + Xf_tp(:,1);
0892 
0893         Xa_tp(:,1) = 0;
0894         Xi_tp(:,1) = 0;
0895         Xf_tp(:,1) = 0;
0896 
0897         Xc_tp(:,nmhu) = Xc_tp(:,nmhu) + Xe_tp(:,nmhu);
0898         Xd_tp(:,nmhu) = Xd_tp(:,nmhu) + Xg_tp(:,nmhu);
0899         Xb_tp(:,nmhu) = Xb_tp(:,nmhu) + Xh_tp(:,nmhu);
0900 
0901         Xe_tp(:,nmhu) = 0;
0902         Xg_tp(:,nmhu) = 0;
0903         Xh_tp(:,nmhu) = 0;
0904 
0905         Xc_tp(1,:) = Xc_tp(1,:) + fliplr(Xb_tp(1,:));
0906         Xa_tp(1,:) = Xa_tp(1,:) + fliplr(Xf_tp(1,:));
0907         Xe_tp(1,:) = Xe_tp(1,:) + fliplr(Xh_tp(1,:));
0908 
0909         Xb_tp(1,:) = 0;
0910         Xh_tp(1,:) = 0;
0911         Xf_tp(1,:) = 0;
0912 %
0913 %Correction matrix resulting from the radial gradient of Chang and Cooper interpolation coefficients
0914 %
0915         dpsim = (xpsin(jr1) - xpsin(jr2));%Radial step for parabolic gradient calculations (dke)
0916         dpsi = (xpsin(jr3) - xpsin(jr1));%Radial step for parabolic gradient calculations (dke)
0917         dpsip = (xpsin(jr3) - xpsin(jr2));%Radial step for parabolic gradient calculations (dke)
0918 %
0919         XXH_tp_mm(:,:,jr1) = (xftp_norm(jr2)*dpsip/dpsi/dpsim)*(XXdeltap_imj(:,:,jr2) - XXdeltap_imj(:,:,jr1)).*XFp_tp_imj.*Xmhu.*Xpn2m.*Xpnm./Xdpn;
0920         XXH_tp_mp(:,:,jr1) = (xftp_norm(jr2)*dpsip/dpsi/dpsim)*(XXdeltap_ipj(:,:,jr2) - XXdeltap_ipj(:,:,jr1)).*XFp_tp_ipj.*Xmhu.*Xpn2p.*Xpnp./Xdpn;
0921         XXH_tp_m(:,:,jr1) = -XXH_tp_mm(:,:,jr1) - XXH_tp_mp(:,:,jr1);
0922 %
0923         XXH_tp_pm(:,:,jr3) = -(xftp_norm(jr2)*dpsim/dpsi/dpsip)*(XXdeltap_imj(:,:,jr2) - XXdeltap_imj(:,:,jr3)).*XFp_tp_imj.*Xmhu.*Xpn2m.*Xpnm./Xdpn;
0924         XXH_tp_pp(:,:,jr3) = -(xftp_norm(jr2)*dpsim/dpsi/dpsip)*(XXdeltap_ipj(:,:,jr2) - XXdeltap_ipj(:,:,jr3)).*XFp_tp_ipj.*Xmhu.*Xpn2p.*Xpnp./Xdpn;       
0925         XXH_tp_p(:,:,jr3) = -XXH_tp_pm(:,:,jr3) - XXH_tp_pp(:,:,jr3);
0926 %
0927 %Matrix coefficients calculations for g
0928 %
0929         Xa_g = XXa(:,:,jr2);%g(i+1/2,j-1/2)
0930         Xb_g = XXb(:,:,jr2);%g(i-1/2,j+1/2)
0931         Xc_g = XXc(:,:,jr2);%g(i+1/2,j+1/2) (Since matrix coefficients for f0 and g are the same)
0932         Xd_g = XXd(:,:,jr2);%g(i+3/2,j+1/2)
0933         Xe_g = XXe(:,:,jr2);%g(i+1/2,j+3/2)
0934         Xf_g = XXf(:,:,jr2);%g(i-1/2,j-1/2)
0935         Xg_g = XXg(:,:,jr2);%g(i+3/2,j+3/2)
0936         Xh_g = XXh(:,:,jr2);%g(i-1/2,j+3/2)
0937         Xi_g = XXi(:,:,jr2);%g(i+3/2,j-1/2)
0938 %
0939 %Boundary conditions for g (the Lorentz model is enforced at p = 0)
0940 %
0941         if boundary_mode_g > 0,
0942             Xa_g(1:boundary_mode_g,:) = 0;%g(i+1/2,j-1/2)
0943             Xb_g(1:boundary_mode_g,:) = 0;%g(i-1/2,j+1/2)
0944             Xc_g(1:boundary_mode_g,:) = 1;%g(i+1/2,j+1/2)
0945             Xd_g(1:boundary_mode_g,:) = 0;%g(i+3/2,j+1/2)
0946             Xe_g(1:boundary_mode_g,:) = 0;%g(i+1/2,j+3/2)
0947             Xf_g(1:boundary_mode_g,:) = 0;%g(i-1/2,j-1/2)
0948             Xg_g(1:boundary_mode_g,:) = 0;%g(i+3/2,j+3/2)
0949             Xh_g(1:boundary_mode_g,:) = 0;%g(i-1/2,j+3/2)
0950             Xi_g(1:boundary_mode_g,:) = 0;%g(i+3/2,j-1/2)
0951         end
0952 %
0953 %Boundaries of the bounce domain for g
0954 %
0955         Xa_g_t = Xa_g(:,Zmask_g{ir_dke});%Matrix reduction along the angular direction
0956         Xb_g_t = Xb_g(:,Zmask_g{ir_dke});
0957         Xc_g_t = Xc_g(:,Zmask_g{ir_dke});
0958         Xd_g_t = Xd_g(:,Zmask_g{ir_dke});
0959         Xe_g_t = Xe_g(:,Zmask_g{ir_dke});
0960         Xf_g_t = Xf_g(:,Zmask_g{ir_dke});
0961         Xg_g_t = Xg_g(:,Zmask_g{ir_dke});
0962         Xh_g_t = Xh_g(:,Zmask_g{ir_dke});
0963         Xi_g_t = Xi_g(:,Zmask_g{ir_dke});
0964 %
0965         Xe_g_t(:,(nmhu - ns_g(ir_dke))/2) = 0;%g(pn,mhu = -mhut) = 0
0966         Xg_g_t(:,(nmhu - ns_g(ir_dke))/2) = 0;%g(pn,mhu = -mhut) = 0
0967         Xh_g_t(:,(nmhu - ns_g(ir_dke))/2) = 0;%g(pn,mhu = -mhut) = 0
0968 %
0969         Xa_g_t(:,(nmhu - ns_g(ir_dke))/2 + 1) = 0;%g(pn,mhu = mhut) = 0
0970         Xf_g_t(:,(nmhu - ns_g(ir_dke))/2 + 1) = 0;%g(pn,mhu = mhut) = 0
0971         Xi_g_t(:,(nmhu - ns_g(ir_dke))/2 + 1) = 0;%g(pn,mhu = mhut) = 0
0972 %
0973 %Matrix calculation for ftp
0974 %
0975         MXc_tp = reshape(Xc_tp',nmhu*npn,1);
0976         MXP_tp = sparse([1:nmhu*npn],[1:nmhu*npn],MXc_tp,nmhu*npn,nmhu*npn);
0977         MXa_tp = reshape(Xa_tp',nmhu*npn,1);
0978         MXP_tp = MXP_tp + sparse([1+1:nmhu*npn],[1:nmhu*npn-1],MXa_tp(1+1:length(MXa_tp)),nmhu*npn,nmhu*npn);
0979         MXe_tp = reshape(Xe_tp',nmhu*npn,1);
0980         MXP_tp = MXP_tp + sparse([1:nmhu*npn-1],[1+1:nmhu*npn],MXe_tp(1:length(MXe_tp)-1),nmhu*npn,nmhu*npn);       
0981         MXb_tp = reshape(Xb_tp',nmhu*npn,1);
0982         MXP_tp = MXP_tp + sparse([1+nmhu:nmhu*npn],[1:nmhu*npn-nmhu],MXb_tp(1+nmhu:length(MXb_tp)),nmhu*npn,nmhu*npn);            
0983         MXd_tp = reshape(Xd_tp',nmhu*npn,1);
0984         MXP_tp = MXP_tp + sparse([1:nmhu*npn-nmhu],[1+nmhu:nmhu*npn],MXd_tp(1:length(MXd_tp)-nmhu),nmhu*npn,nmhu*npn);            
0985         MXf_tp = reshape(Xf_tp',nmhu*npn,1);
0986         MXP_tp = MXP_tp + sparse([1+nmhu+1:nmhu*npn],[1:nmhu*npn-nmhu-1],MXf_tp(1+nmhu+1:length(MXf_tp)),nmhu*npn,nmhu*npn);      
0987         MXh_tp = reshape(Xh_tp',nmhu*npn,1);
0988         MXP_tp = MXP_tp + sparse([1+nmhu-1:nmhu*npn],[1:nmhu*npn-nmhu+1],MXh_tp(1+nmhu-1:length(MXh_tp)),nmhu*npn,nmhu*npn);      
0989         MXi_tp = reshape(Xi_tp',nmhu*npn,1);
0990         MXP_tp = MXP_tp + sparse([1:nmhu*npn-nmhu+1],[1+nmhu-1:nmhu*npn],MXi_tp(1:length(MXi_tp)-nmhu+1),nmhu*npn,nmhu*npn);      
0991         MXg_tp = reshape(Xg_tp',nmhu*npn,1);
0992         MXP_tp = MXP_tp + sparse([1:nmhu*npn-nmhu-1],[1+nmhu+1:nmhu*npn],MXg_tp(1:length(MXg_tp)-nmhu-1),nmhu*npn,nmhu*npn);                
0993 %
0994         sm_tp(ir_dke+1) = sm_tp(ir_dke) + blocksize_tp(ir_dke);
0995         [i0_tp,j0_tp,s0_tp] = find(MMXP_tp);        
0996         [m0_tp,n0_tp] = size(MMXP_tp);
0997         [i1_tp,j1_tp,s1_tp] = find(MXP_tp);
0998         [m1_tp,n1_tp] = size(MXP_tp);
0999         MMXP_tp = sparse(i0_tp,j0_tp,s0_tp,m0_tp+m1_tp,m0_tp+m1_tp) + sparse(m0_tp+i1_tp,m0_tp+j1_tp,s1_tp,m0_tp+m1_tp,m0_tp+m1_tp);
1000 %
1001 %Matrix calculation for g
1002 %
1003         MXc_g_t = reshape(Xc_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1004         MXP_g_t = sparse([1:(nmhu - ns_g(ir_dke))*npn],[1:(nmhu - ns_g(ir_dke))*npn],MXc_g_t,(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);
1005         MXa_g_t = reshape(Xa_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1006         MXP_g_t = MXP_g_t + sparse([1+1:(nmhu - ns_g(ir_dke))*npn],[1:(nmhu - ns_g(ir_dke))*npn-1],MXa_g_t(1+1:length(MXa_g_t)),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);
1007         MXe_g_t = reshape(Xe_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1008         MXP_g_t = MXP_g_t + sparse([1:(nmhu - ns_g(ir_dke))*npn-1],[1+1:(nmhu - ns_g(ir_dke))*npn],MXe_g_t(1:length(MXe_g_t)-1),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);       
1009         MXb_g_t = reshape(Xb_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1010         MXP_g_t = MXP_g_t + sparse([1+(nmhu - ns_g(ir_dke)):(nmhu - ns_g(ir_dke))*npn],[1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))],MXb_g_t(1+(nmhu - ns_g(ir_dke)):length(MXb_g_t)),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);           
1011         MXd_g_t = reshape(Xd_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1012         MXP_g_t = MXP_g_t + sparse([1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))],[1+(nmhu - ns_g(ir_dke)):(nmhu - ns_g(ir_dke))*npn],MXd_g_t(1:length(MXd_g_t)-(nmhu - ns_g(ir_dke))),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);           
1013         MXf_g_t = reshape(Xf_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1014         MXP_g_t = MXP_g_t + sparse([1+(nmhu - ns_g(ir_dke))+1:(nmhu - ns_g(ir_dke))*npn],[1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))-1],MXf_g_t(1+(nmhu - ns_g(ir_dke))+1:length(MXf_g_t)),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);     
1015         MXh_g_t = reshape(Xh_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1016         MXP_g_t = MXP_g_t + sparse([1+(nmhu - ns_g(ir_dke))-1:(nmhu - ns_g(ir_dke))*npn],[1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))+1],MXh_g_t(1+(nmhu - ns_g(ir_dke))-1:length(MXh_g_t)),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);     
1017         MXi_g_t = reshape(Xi_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1018         MXP_g_t = MXP_g_t + sparse([1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))+1],[1+(nmhu - ns_g(ir_dke))-1:(nmhu - ns_g(ir_dke))*npn],MXi_g_t(1:length(MXi_g_t)-(nmhu - ns_g(ir_dke))+1),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);     
1019         MXg_g_t = reshape(Xg_g_t',(nmhu - ns_g(ir_dke))*npn,1);
1020         MXP_g_t = MXP_g_t + sparse([1:(nmhu - ns_g(ir_dke))*npn-(nmhu - ns_g(ir_dke))-1],[1+(nmhu - ns_g(ir_dke))+1:(nmhu - ns_g(ir_dke))*npn],MXg_g_t(1:length(MXg_g_t)-(nmhu - ns_g(ir_dke))-1),(nmhu - ns_g(ir_dke))*npn,(nmhu - ns_g(ir_dke))*npn);     
1021 %
1022         sm_g(ir_dke+1) = sm_g(ir_dke) + blocksize_g_t(ir_dke);
1023         [i0_g,j0_g,s0_g] = find(MMXP_g_t);      
1024         [m0_g,n0_g] = size(MMXP_g_t);
1025         [i1_g,j1_g,s1_g] = find(MXP_g_t);
1026         [m1_g,n1_g] = size(MXP_g_t);
1027         MMXP_g_t = sparse(i0_g,j0_g,s0_g,m0_g+m1_g,m0_g+m1_g) + sparse(m0_g+i1_g,m0_g+j1_g,s1_g,m0_g+m1_g,m0_g+m1_g);
1028 %
1029         etime_dke = etime_dke + etime(clock,time0);
1030        time0 = clock;
1031     end
1032       
1033       
1034 else
1035 %
1036 % for collisions
1037 %
1038     XXDpp_c_tp_imj = zeros(npn,nmhu,nr);
1039     XXDpm_c_tp_imj = zeros(npn,nmhu,nr);
1040     XXFp_c_tp_imj = zeros(npn,nmhu,nr);
1041     XXDpp_c_tp_ipj = zeros(npn,nmhu,nr);
1042     XXDpm_c_tp_ipj = zeros(npn,nmhu,nr);
1043     XXFp_c_tp_ipj = zeros(npn,nmhu,nr);
1044 %
1045     XXDmm_c_tp_ijm = zeros(npn,nmhu,nr);
1046     XXDmp_c_tp_ijm = zeros(npn,nmhu,nr);
1047     XXFm_c_tp_ijm = zeros(npn,nmhu,nr);
1048     XXDmm_c_tp_ijp = zeros(npn,nmhu,nr);
1049     XXDmp_c_tp_ijp = zeros(npn,nmhu,nr);
1050     XXFm_c_tp_ijp = zeros(npn,nmhu,nr);
1051 %
1052 % For Ohmic electric field
1053 %
1054     XXDpp_e_tp_imj = zeros(npn,nmhu,nr);
1055     XXDpm_e_tp_imj = zeros(npn,nmhu,nr);
1056     XXFp_e_tp_imj = zeros(npn,nmhu,nr);
1057     XXDpp_e_tp_ipj = zeros(npn,nmhu,nr);
1058     XXDpm_e_tp_ipj = zeros(npn,nmhu,nr);
1059     XXFp_e_tp_ipj = zeros(npn,nmhu,nr);
1060 %
1061     XXDmm_e_tp_ijm = zeros(npn,nmhu,nr);
1062     XXDmp_e_tp_ijm = zeros(npn,nmhu,nr);
1063     XXFm_e_tp_ijm = zeros(npn,nmhu,nr);
1064     XXDmm_e_tp_ijp = zeros(npn,nmhu,nr);
1065     XXDmp_e_tp_ijp = zeros(npn,nmhu,nr);
1066     XXFm_e_tp_ijp = zeros(npn,nmhu,nr);
1067 %
1068 % For Synchrotron reaction
1069 %
1070     XXDpp_s_tp_imj = zeros(npn,nmhu,nr);
1071     XXDpm_s_tp_imj = zeros(npn,nmhu,nr);
1072     XXFp_s_tp_imj = zeros(npn,nmhu,nr);
1073     XXDpp_s_tp_ipj = zeros(npn,nmhu,nr);
1074     XXDpm_s_tp_ipj = zeros(npn,nmhu,nr);
1075     XXFp_s_tp_ipj = zeros(npn,nmhu,nr);
1076 %
1077     XXDmm_s_tp_ijm = zeros(npn,nmhu,nr);
1078     XXDmp_s_tp_ijm = zeros(npn,nmhu,nr);
1079     XXFm_s_tp_ijm = zeros(npn,nmhu,nr);
1080     XXDmm_s_tp_ijp = zeros(npn,nmhu,nr);
1081     XXDmp_s_tp_ijp = zeros(npn,nmhu,nr);
1082     XXFm_s_tp_ijp = zeros(npn,nmhu,nr);
1083 %
1084 end
1085 %
1086 %Second order radial gradient correction for ftp (radial derivative of the Chang and Cooper momentum terms)
1087 %
1088 if dke_mode == 1,
1089 %    XXH_tp_mm = permute(XXH_tp_mm,[2,1,3]);
1090 %    XXH_tp_m = permute(XXH_tp_m,[2,1,3]);
1091 %    XXH_tp_mp = permute(XXH_tp_mp,[2,1,3]);
1092 %    XXH_tp_pm = permute(XXH_tp_pm,[2,1,3]);
1093 %    XXH_tp_p = permute(XXH_tp_p,[2,1,3]);
1094 %    XXH_tp_pp = permute(XXH_tp_pp,[2,1,3]);
1095 %
1096 %    MMH_tp_mm = reshape(XXH_tp_mm,npn*nmhu*nr,1,1);
1097 %    MMH_tp_m = reshape(XXH_tp_m,npn*nmhu*nr,1,1);
1098 %    MMH_tp_mp = reshape(XXH_tp_mp,npn*nmhu*nr,1,1);
1099 %    MMH_tp_pm = reshape(XXH_tp_pm,npn*nmhu*nr,1,1);
1100 %    MMH_tp_p = reshape(XXH_tp_p,npn*nmhu*nr,1,1);
1101 %    MMH_tp_pp = reshape(XXH_tp_pp,npn*nmhu*nr,1,1);
1102 %
1103 %    MMH_tp = MMH_tp + sparse([1+(nmhu-1)*npn:nMMH_tp],[1:nMMH_tp-(nmhu-1)*npn],MMH_tp_mm(1+npn*(nmhu-1):nMMH_tp),nMMH_tp,nMMH_tp);
1104 %    MMH_tp = MMH_tp + sparse([1+nmhu*npn:nMMH_tp],[1:nMMH_tp-nmhu*npn],MMH_tp_m(1+npn*nmhu:nMMH_tp),nMMH_tp,nMMH_tp);
1105 %    MMH_tp = MMH_tp + sparse([1+(nmhu+1)*npn:nMMH_tp],[1:nMMH_tp-(nmhu+1)*npn],MMH_tp_mp(1+npn*(nmhu+1):nMMH_tp),nMMH_tp,nMMH_tp);
1106 %    MMH_tp = MMH_tp + sparse([1:nMMH_tp-(nmhu-1)*npn],[1+(nmhu-1)*npn:nMMH_tp],MMH_tp_pm(1:nMMH_tp-npn*(nmhu-1)),nMMH_tp,nMMH_tp);
1107 %    MMH_tp = MMH_tp + sparse([1:nMMH_tp-nmhu*npn],[1+nmhu*npn:nMMH_tp],MMH_tp_p(1:nMMH_tp-npn*nmhu),nMMH_tp,nMMH_tp);
1108 %    MMH_tp = MMH_tp + sparse([1:nMMH_tp-(nmhu+1)*npn],[1+(nmhu+1)*npn:nMMH_tp],MMH_tp_pp(1:nMMH_tp-npn*(nmhu+1)),nMMH_tp,nMMH_tp);
1109 %
1110 %    if display_mode >= 1,info_dke_yp(2,['Second order radial gradient correction for ftp matrix calculations done in ',num2str(etime(clock,time0)),' (s)']);time0 = clock;end
1111 end
1112 %
1113 
1114