calc_fcyl_jd

 this function transforms the distribution into cylindrical momentum space
 and calculates the parallel distribution function

 here the momentum grid pn is normalized to the reference thermal momentum

0001 function [XXf_cyl,xXFpar_cyl,ppar_cyl,xXFperp_cyl,pperp_cyl] = calc_fcyl_jd(XXf,betath_ref,pn,mhu,dp_cyl,ir_list,iir_display,pnmax_disp,XXmask,pmhucrit,ifig)
0002 %
0003 % this function transforms the distribution into cylindrical momentum space
0004 % and calculates the parallel distribution function
0005 %
0006 % here the momentum grid pn is normalized to the reference thermal momentum
0007 %
0008 if nargin < 9,
0009     ifig = [1,2,3];
0010 end
0011 %
0012 if nargin < 8,
0013     pnmax_disp = [];
0014 end
0015 %
0016 if nargin < 7,
0017     iir_display = [];
0018 end
0019 %
0020 if nargin < 6,
0021     ir_list = 1:size(XXf,3);
0022 end
0023 %
0024 if nargin < 5,
0025     dp_cyl = 0.1;
0026 end
0027 %
0028 if nargin < 4,
0029     error('Not enough input arguments');
0030 end
0031 %
0032 nr = length(ir_list);
0033 %
0034 % Cylindrical grid (also see s2c_dke_yp)
0035 %
0036 pnmask = pn <= pnmax_disp;
0037 pn = pn(pnmask);
0038 XXf = XXf(pnmask,:,:);
0039 XXmask = XXmask(pnmask,:,:);
0040 %
0041 pmax = max(pn) + max(diff(pn));
0042 np_cyl = ceil(pmax/dp_cyl);
0043 ppar_cyl = linspace(-pmax,pmax,2*np_cyl + 1);
0044 pperp_cyl = linspace(0,pmax,np_cyl + 1);
0045 %
0046 XXf_cyl = NaN(np_cyl + 1,2*np_cyl + 1,nr);
0047 XXmask_cyl = NaN(np_cyl + 1,2*np_cyl + 1,nr);
0048 %
0049 np = length(pn);
0050 Xxmhu_crit = NaN(np,nr);
0051 %
0052 % XXf(XXf <= 0) = 0;
0053 %XXf(XXf <= eps) = 0;
0054 %
0055 for iir = 1:nr, 
0056     %
0057     % spherical to cylindrical transformation
0058     %
0059     logXf_cyl = s2c_dke_yp(log(XXf(:,:,ir_list(iir))),pn,mhu,dp_cyl);%Spherical to cylindrical coordinate transformation
0060     Xf_cyl = exp(logXf_cyl);%For accurate representation in figures
0061     %
0062     Xmask_cyl = s2c_dke_yp(double(XXmask(:,:,ir_list(iir))),pn,mhu,dp_cyl);%Spherical to cylindrical coordinate transformation
0063     %
0064     Xppar_cyl = repmat(ppar_cyl,[np_cyl + 1,1]);%Grid for the distribution functions only
0065     Xpperp_cyl = repmat(pperp_cyl.',[1,2*np_cyl + 1]);
0066     Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0067     %
0068 %     Xf_cyl(Xf_cyl <= 0) = 0;%Remove negative values and values less than computer accuracy
0069     Xf_cyl(Xf_cyl <= eps) = 0;%Remove negative values and values less than computer accuracy
0070     Xf_cyl(Xp_cyl > pmax) = 0;%Remove values above max(ppar)
0071     Xf_cyl(isnan(Xf_cyl)) = 0;%Remove NaN values
0072     %
0073     XXf_cyl(:,:,iir) = Xf_cyl;
0074     XXmask_cyl(:,:,iir) = Xmask_cyl;
0075     %
0076     for ipn = 1:np,
0077         Xxmhu_crit(ipn,ir_list(iir)) = min([1,mhu(XXmask(ipn,:,ir_list(iir)))]);
0078     end
0079     %
0080 end
0081 %
0082 % parallel distribution
0083 %
0084 XXfh_perp_cyl = (XXf_cyl(1:np_cyl,:,:) + XXf_cyl(2:np_cyl + 1,:,:))/2;
0085 XXfh_par_cyl = (XXf_cyl(:,1:2*np_cyl,:) + XXf_cyl(:,2:2*np_cyl + 1,:))/2;
0086 XXpperph_cyl = repmat((pperp_cyl(1:np_cyl) + pperp_cyl(2:np_cyl + 1)).'/2,[1,2*np_cyl + 1,nr]);
0087 dpperp_cyl = diff(pperp_cyl);
0088 dppar_cyl = diff(ppar_cyl);
0089 %
0090 xXFpar_cyl = 2*pi*integral_dke_jd(dpperp_cyl,XXpperph_cyl.*XXfh_perp_cyl,1);
0091 xXFperp_cyl = 2*pi*integral_dke_jd(dppar_cyl,XXfh_par_cyl,2);
0092 %
0093 if isempty(iir_display), 
0094     inter = 1;
0095     iir_display = input_dke_yp(['Select radial location index (1:',num2str(nr),'). Enter 0 to exit. '],1,0:nr);
0096 else
0097     inter = 0;
0098 end
0099 %
0100 while iir_display > 0,
0101     %
0102     Xf_cyl = XXf_cyl(:,:,iir_display);
0103     Xmask_cyl = XXmask_cyl(:,:,iir_display);
0104     Fpar_cyl = xXFpar_cyl(:,iir_display);
0105     Fperp_cyl = xXFperp_cyl(:,iir_display);
0106     %
0107     if isempty(pnmax_disp),
0108         pnmax = pmax;
0109     elseif isnan(pnmax_disp),
0110         pnmax = min([pmax,min(min(Xp_cyl(Xf_cyl == 0)))]);
0111     else
0112         pnmax = pnmax_disp;
0113     end
0114     %
0115     Xf_cyl(Xp_cyl > pnmax) = 0;
0116     %
0117 %     x1 = ppar_cyl;
0118 %     xlab1 = 'p_{||}/p_{Te}';
0119 %     xlim1 = [-pnmax,pnmax];
0120     x1 = ppar_cyl*betath_ref;
0121     xlab1 = 'p_{||}';%/(mc)
0122     xlim1 = [-pnmax,pnmax]*betath_ref;
0123     %
0124 %     x2 = pperp_cyl;
0125 %     xlab2 = 'p_{\perp}/p_{Te}';
0126 %     xlim2 = [0,pnmax];
0127     x2 = pperp_cyl*betath_ref;
0128     xlab2 = 'p_{\perp}';%/(mc)
0129     xlim2 = [0,pnmax]*betath_ref;
0130     %
0131     y2 = sqrt(pnmax^2 - ppar_cyl.^2)*betath_ref;
0132     %
0133     if isempty(pmhucrit) || all(isnan(pmhucrit)),
0134         pmhucrit = Xxmhu_crit(:,iir_display).';
0135     else
0136         pmhucrit = pmhucrit(pnmask);
0137     end
0138     x4 = pn*betath_ref.*pmhucrit;
0139     y4 = pn*betath_ref.*sqrt(1 - pmhucrit.^2);
0140     %
0141     % ---------- 2D distribution function -----------
0142     %
0143     if ifig(1) > 0,
0144         figure(ifig(1)),clf,set(ifig(1),'Name',['2-D contour of f, iir = ',num2str(iir_display)])
0145         %figure(ifig(1)),hold on,set(1,'Name',['2-D contour of f, iir = ',num2str(iir_display)])
0146         %
0147         y1 = Xf_cyl;
0148         y3 = Xmask_cyl;
0149         %
0150         graph2D_jd(x1,x2,y1.','','','',xlim1,xlim2,0,betath_ref,0:0.2:40,0,'-','k',0.5,20,1,max(max(y1)));hold on
0151 %         graph2D_jd(x1,x2,y3.','','','',xlim1,xlim2,0,NaN,0.5,0,'-','k',2,20);
0152         graph1D_jd(x4,y4,0,0,'','','',NaN,xlim1,xlim2,'--','none','k',0.5,20);
0153         %
0154         graph1D_jd(x1,y2,0,0,xlab1,xlab2,'',NaN,xlim1,xlim2,'--','none','k',2,20,gca,0.9,0.7,0.8);
0155         %
0156         axis('equal');axis([xlim1,xlim2]);
0157     end
0158     %
0159     % ---------- 1D par distribution function -----------
0160     %
0161     if ifig(2) > 0,
0162         figure(ifig(2)),clf,set(ifig(2),'Name',['F//, iir = ',num2str(iir_display)])
0163         %
0164         y = Fpar_cyl;
0165         ylab = 'F_{||}';
0166         ylim = [eps,1];
0167         %
0168         graph1D_jd(x1,y,0,1,xlab1,ylab,'',NaN,xlim1,ylim,'-','none','r',2,20,gca,0.9,0.7,0.7);
0169         %
0170         set(gca,'ytick',10.^(-16:2:0))
0171     end
0172     %
0173     %
0174     % ---------- 1D perp distribution function -----------
0175     %
0176     if ifig(3) > 0,
0177         figure(ifig(3)),clf,set(ifig(3),'Name',['Fperp, iir = ',num2str(iir_display)])
0178         %
0179         y = Fperp_cyl;
0180         ylab = 'F_{\perp}';
0181         ylim = [eps,1];
0182         %
0183         graph1D_jd(x2,y,0,1,xlab2,ylab,'',NaN,xlim2,ylim,'-','none','r',2,20,gca,0.9,0.7,0.7);
0184         %
0185         set(gca,'ytick',10.^(-16:2:0))
0186         %
0187         if inter == 1,
0188             iir_display = input_dke_yp(['Select radial location index (1:',num2str(nr),'). Enter 0 to exit. '],0,0:nr);
0189         end
0190     end
0191     %
0192     if inter == 0,
0193         iir_display = 0;
0194     end
0195     %
0196 end