stream_dke_jd

Stream function calculation for graphical flux analysis

by Joan Decker (joan.decker@cea.fr, CEA/DSM/IRFM) and Yves Peysson (yves.peysson@cea.fr, CEA/DSM/IRFM)

0001 function [XAstr_S,XAstr_cyl_S,ppar_cyl_S,pperp_cyl_S] = stream_dke_jd(dke_out,dkeparam,momentumDKE,gridDKE,ir_display_dke,dp_cyl,graph,XXmask,pmhu_crit,ifig)
0002 %
0003 %Stream function calculation for graphical flux analysis
0004 %
0005 %by Joan Decker (joan.decker@cea.fr, CEA/DSM/IRFM) and Yves Peysson (yves.peysson@cea.fr, CEA/DSM/IRFM)
0006 %
0007 if nargin < 9,
0008     ifig = 0;
0009 end
0010 %
0011 if nargin < 8,
0012     XXmask = [];
0013 end
0014 %
0015 if nargin < 7,
0016     graph = struct;
0017 end
0018 %
0019 if nargin < 6,
0020     dp_cyl = [];
0021 end
0022 %
0023 if nargin < 5,
0024     ir_display_dke = 1;
0025 end
0026 %
0027 if nargin < 4,
0028     info_dke_yp(2,'Wrong number of input arguments for stream_dke_jd');
0029     return;
0030 end
0031 %
0032 if ~isfield(graph,'pnmax_disp') || isnan(graph.pnmax_disp),
0033     graph.pnmax_disp = momentumDKE.pn_S(end);
0034 end
0035 %
0036 ZXXS = dke_out.ZXXS;
0037 xRR_fsav = 0;%dke_out.xRR_fsav(ir_display_dke);
0038 betath_ref = dke_out.betath_ref;
0039 %
0040 pn_S = momentumDKE.pn_S.';
0041 npn_S = find(pn_S <= graph.pnmax_disp,1,'last');
0042 pn_S = pn_S(1:npn_S);
0043 %
0044 mhu_S = momentumDKE.mhu_S;
0045 nmhu_S = length(mhu_S);
0046 %
0047 dmhu = momentumDKE.dmhu;
0048 nmhu = length(dmhu);
0049 xmhubounce2 = gridDKE.xmhubounce2(ir_display_dke);
0050 %
0051 XSp0_tot_ipj = ZXXS.p0_c_ipj(1:npn_S-1,:,ir_display_dke) + ZXXS.p0_e_ipj(1:npn_S-1,:,ir_display_dke) + ...
0052     ZXXS.p0_s_ipj(1:npn_S-1,:,ir_display_dke) + ZXXS.p0_rf_ipj(1:npn_S-1,:,ir_display_dke);
0053 %
0054 if ~isempty(XXmask),
0055     pn = momentumDKE.pn.';
0056     pnmask = pn <= graph.pnmax_disp;
0057     pn = pn(pnmask);
0058     XXmask = XXmask(pnmask,:,:);
0059     mhu = momentumDKE.mhu;
0060     np = length(pn);
0061 end
0062 %
0063 XAstr_S_Pp = zeros(npn_S,nmhu_S);%circulating particles (pn>=0)
0064 XAstr_S_Pm = zeros(npn_S,nmhu_S);%circulating particles (pn<0)
0065 XAstr_S_T = zeros(npn_S,nmhu_S);%trapped particles
0066 %
0067 if dkeparam.bounce_mode,
0068     %
0069     nT = sum(abs(mhu_S) <= sqrt(xmhubounce2));%Trapped particles
0070     nP = (nmhu_S - nT)/2;%Circulating particles
0071     %
0072     for j=1:nP,%ksi < 0
0073         XAstr_S_Pm(2:npn_S,j) = xRR_fsav - 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu(1:j),XSp0_tot_ipj(:,1:j),2);%Stream function with runaway
0074     end
0075     %
0076     for j=1:nP,%ksi > 0
0077         XAstr_S_Pp(2:npn_S,nmhu_S-j+1) = xRR_fsav + 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu(nmhu_S-j:nmhu_S-1),XSp0_tot_ipj(:,nmhu_S-j:nmhu_S-1),2);%Stream function with runaway
0078     end
0079     %
0080     for j=1:(nT-1)/2,%|ksi| <= ksi_T
0081         XAstr_S_T(2:npn_S,(nmhu_S-1)/2+j+1) = xRR_fsav - 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu((nmhu_S-1)/2+1:(nmhu_S-1)/2+j),XSp0_tot_ipj(:,(nmhu_S-1)/2+1:(nmhu_S-1)/2+j),2);%Stream function with runaway
0082         XAstr_S_T(2:npn_S,(nmhu_S-1)/2-j+1) = xRR_fsav + 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu((nmhu_S-1)/2+1-j:(nmhu_S-1)/2),XSp0_tot_ipj(:,(nmhu_S-1)/2+1-j:(nmhu_S-1)/2),2);%Stream function with runaway
0083     end
0084     %
0085     XAstr_S_T(2:npn_S,(nmhu_S+1)/2) = xRR_fsav;
0086     %
0087     XAstr_S_T = (XAstr_S_T + fliplr(XAstr_S_T))/2;
0088     XAstr_S = XAstr_S_Pm + XAstr_S_T + XAstr_S_Pp;
0089 else
0090     %
0091     pRR_ip_fsav = integral_dke_jd(dmhu,XSp0_tot_ipj,2)/2;
0092     %
0093     XSp0_tot_ipj = XSp0_tot_ipj - repmat(pRR_ip_fsav,[1,nmhu]);
0094     %
0095     for j=1:nmhu_S-1,
0096         XAstr_S_Pm(2:npn_S,j+1) = xRR_fsav - 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu(1:j),XSp0_tot_ipj(:,1:j),2);%Stream function with runaway
0097     end
0098     %
0099     for j=1:nmhu_S-1,
0100         XAstr_S_Pp(2:npn_S,nmhu_S-j) = xRR_fsav + 2*pi*pn_S(2:npn_S).*pn_S(2:npn_S).*integral_dke_jd(dmhu(nmhu_S-j:nmhu_S-1),XSp0_tot_ipj(:,nmhu_S-j:nmhu_S-1),2);%Stream function with runaway
0101     end
0102     %
0103     XAstr_S = (XAstr_S_Pm + XAstr_S_Pp)/2;
0104 end
0105 %
0106 if isempty(dp_cyl),
0107     dp_cyl = min(pn_S(2:end) - pn_S(1:end-1))/2;%cylindrical momentum grid step
0108 end
0109 %
0110 [XAstr_cyl_S,ppar_cyl_S,~,pperp_cyl_S] = s2c_dke_yp(XAstr_S,pn_S,mhu_S,dp_cyl);%Spherical to cylindrical coordinate transformation
0111 % [XAstr_cyl_S_Pm] = s2c_dke_yp(XAstr_S_Pm,pn_S,mhu_S,dp_cyl);%Spherical to cylindrical coordinate transformation
0112 % [XAstr_cyl_S_Pp] = s2c_dke_yp(XAstr_S_Pp,pn_S,mhu_S,dp_cyl);%Spherical to cylindrical coordinate transformation
0113 % [logXAstr_cyl_S,ppar_cyl_S,dppar_cyl_S,pperp_cyl_S,dpperp_S_cyl] = s2c_dke_yp(log(XAstr_S),pn_S,mhu_S,dp_cyl);%Spherical to cylindrical coordinate transformation
0114 % XAstr_cyl_S = exp(logXAstr_cyl_S);%For accurate representation in figures
0115 %
0116 if ~isempty(XXmask),
0117     [Xmask_cyl,ppar_cyl,~,pperp_cyl] = s2c_dke_yp(double(XXmask(:,:,ir_display_dke)),pn,mhu,dp_cyl);%Spherical to cylindrical coordinate transformation
0118 else
0119     Xmask_cyl = NaN;
0120 end
0121 %
0122 if isempty(pmhu_crit) || all(isnan(pmhu_crit)),
0123     pmhu_crit = zeros(1,np);
0124     for ipn = 1:np,
0125         pmhu_crit(ipn) = min([1,mhu(XXmask(ipn,:,ir_display_dke))]);
0126     end
0127 else
0128     pmhu_crit = pmhu_crit(pnmask);
0129 end
0130 %
0131 Xppar_cyl_S = ones(length(pperp_cyl_S),1)*ppar_cyl_S(:)';%Grid for the stream or flux functions only
0132 Xpperp_cyl_S = pperp_cyl_S(:)*ones(1,length(ppar_cyl_S));
0133 Xp_cyl_S = sqrt(Xppar_cyl_S.*Xppar_cyl_S + Xpperp_cyl_S.*Xpperp_cyl_S);
0134 %
0135 pnmax_S = max(ppar_cyl_S);
0136 %
0137 XAstr_cyl_S(Xp_cyl_S > pnmax_S) = 0;%Remove values above max(ppar)
0138 XAstr_cyl_S(isnan(XAstr_cyl_S)) = 0;%Remove NaN values
0139 XAstr_cyl_S(1,:) = 0;
0140 %
0141 % XAstr_cyl_S_Pm(Xp_cyl_S > pnmax_S) = 0;%Remove values above max(ppar)
0142 % XAstr_cyl_S_Pm(isnan(XAstr_cyl_S_Pm)) = 0;%Remove NaN values
0143 % XAstr_cyl_S_Pm(1,:) = 0;
0144 % %
0145 % XAstr_cyl_S_Pp(Xp_cyl_S > pnmax_S) = 0;%Remove values above max(ppar)
0146 % XAstr_cyl_S_Pp(isnan(XAstr_cyl_S_Pp)) = 0;%Remove NaN values
0147 % XAstr_cyl_S_Pp(1,:) = 0;
0148 %
0149 if ifig > 0,
0150     %
0151 %     x1 = ppar_cyl_S;
0152 %     xlab1 = 'p_{||}/p_{Te}';
0153 %     xlim1 = [-pnmax_S,pnmax_S];
0154     x1 = ppar_cyl_S*betath_ref;
0155     xlab1 = 'p_{||}';%/(mc)
0156     xlim1 = [-pnmax_S,pnmax_S]*betath_ref;
0157     %
0158 %     x2 = pperp_cyl_S;
0159 %     xlab2 = 'p_{\perp}/p_{Te}';
0160 %     xlim2 = [0,pnmax_S];
0161     x2 = pperp_cyl_S*betath_ref;
0162     xlab2 = 'p_{\perp}';%/(mc)
0163     xlim2 = [0,pnmax_S]*betath_ref;
0164     %
0165 %     y2 = sqrt(pnmax_S^2 - ppar_cyl_S.^2);
0166     y2 = sqrt(pnmax_S^2 - ppar_cyl_S.^2)*betath_ref;
0167     y3 = Xmask_cyl;
0168     %
0169     x4 = pn.'*betath_ref.*pmhu_crit;
0170     y4 = pn.'*betath_ref.*sqrt(1 - pmhu_crit.^2);
0171     %
0172 %     y1 = log10(XAstr_cyl_S);
0173 %     cont = linspace(min(min(y1(~isinf(y1)))),max(max(y1(~isinf(y1)))),30);
0174     %
0175     y1 = XAstr_cyl_S;
0176     if ~isfield(graph,'cont'),
0177         graph.cont = linspace(min(min(y1(~isinf(y1)))),max(max(y1(~isinf(y1)))),30);
0178     end
0179     %
0180     figure(ifig),hold on%,clf%
0181 %     %
0182 %     tit = '(m+p)/2';
0183 %     y1 = log(XAstr_cyl_S);
0184 %     %
0185 %     graph2D_jd(x1,x2,y1.','','','',xlim1,xlim2,1,NaN,cont,0,'-','k',0.5,20);hold on
0186 %     graph2D_jd(ppar_cyl,pperp_cyl,y3.','','','',xlim1,xlim2,0,NaN,0.5,0,'-','k',2,20);
0187 %     %
0188 %     graph1D_jd(x1,y2,0,0,xlab1,xlab2,tit,NaN,xlim1,xlim2,'--','none','k',2,20,gca,0.9,0.7,0.7);
0189 %     %
0190 %     axis('equal');axis([xlim1,xlim2]);
0191 %     %
0192 %     figure(ifig+1),clf
0193 %     %
0194 %     tit = 'm';
0195 %     y1 = log(XAstr_cyl_S_Pm);
0196 %     %
0197 %     graph2D_jd(x1,x2,y1.','','','',xlim1,xlim2,1,NaN,cont,0,'-','k',0.5,20);hold on
0198 %     graph2D_jd(ppar_cyl,pperp_cyl,y3.','','','',xlim1,xlim2,0,NaN,0.5,0,'-','k',2,20);
0199 %     %
0200 %     graph1D_jd(x1,y2,0,0,xlab1,xlab2,tit,NaN,xlim1,xlim2,'--','none','k',2,20,gca,0.9,0.7,0.7);
0201 %     %
0202 %     axis('equal');axis([xlim1,xlim2]);
0203 %     %
0204 %     figure(ifig+2),clf
0205     %
0206     tit = '';%'p';
0207     %
0208     graph2D_jd(x1,x2,y1.','','','',xlim1,xlim2,0,NaN,graph.cont,0,'-','k',0.5,20);hold on
0209     %
0210 %     graph2D_jd(x1,x2,y3.','','','',xlim1,xlim2,0,NaN,0.5,0,'-','k',2,20);
0211     graph1D_jd(x4,y4,0,0,'','','',NaN,xlim1,xlim2,'--','none','k',0.5,20);
0212         %
0213     graph1D_jd(x1,y2,0,0,xlab1,xlab2,tit,NaN,xlim1,xlim2,'--','none','k',2,20,gca,0.9,0.7,0.7);
0214     %
0215     axis('equal');axis([xlim1,xlim2]);
0216 end    
0217     
0218     
0219     
0220