s2c_dke_yp

    Spherical to cylindrical coordinate transformation with symmetry
    enforced at abs(mhu) = 1. 

    Input:

        - f0:matrix [np,nmhu]  
        - p: momentum vector [1,np]
        - mhu: cos(theta) vector [1,nmhu]
        - dp       

    Output:

        - f1:matrix [2*np_cyl-1,np_cyl]
        - ppar: parallel momentum [1,np_cyl]
        - dppar: parallel momentum increment [1,2*np_cyl-1]
        - pperp: perpendicular momentum [1,np_cyl]
        - dpperp: perpendicular momentum increment [1,np_cyl]
        
by Y.Peysson CEA-DRFC <yves.peysson@cea.fr>

0001 function [f1,ppar,dppar,pperp,dpperp] = s2c_dke_yp(f0,p,mhu,dp)
0002 %
0003 %    Spherical to cylindrical coordinate transformation with symmetry
0004 %    enforced at abs(mhu) = 1.
0005 %
0006 %    Input:
0007 %
0008 %        - f0:matrix [np,nmhu]
0009 %        - p: momentum vector [1,np]
0010 %        - mhu: cos(theta) vector [1,nmhu]
0011 %        - dp
0012 %
0013 %    Output:
0014 %
0015 %        - f1:matrix [2*np_cyl-1,np_cyl]
0016 %        - ppar: parallel momentum [1,np_cyl]
0017 %        - dppar: parallel momentum increment [1,2*np_cyl-1]
0018 %        - pperp: perpendicular momentum [1,np_cyl]
0019 %        - dpperp: perpendicular momentum increment [1,np_cyl]
0020 %
0021 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr>
0022 %
0023 if nargin < 4,
0024     error('Wrong number of input arguments for s2c_dke_yp');
0025     return;
0026 end    
0027 %
0028 ppar0 = p(:)*mhu(:)';
0029 pperp0 = p(:)*sqrt(1-(mhu(:).*mhu(:)))';
0030 sf0 = size(f0,2);
0031 %
0032 %Symmetry at abs(mhu) = 1 enforced. Several point used for the cubic spline
0033 %method, otherwise, the derivative is not correctly enforced (numerical
0034 %noise)
0035 %
0036 nextra = 3;%Change extra point for avoiding bad contour plot with refine mesh
0037 ppar00 = [ppar0(:,nextra:-1:1),ppar0,ppar0(:,sf0:-1:sf0-nextra+1)];
0038 pperp00 = [-pperp0(:,nextra:-1:1),pperp0,-pperp0(:,sf0:-1:sf0-nextra+1)];
0039 f00 = [f0(:,nextra:-1:1),f0,f0(:,sf0:-1:sf0-nextra+1)];
0040 %
0041 pmax = max(p) + max(diff(p));
0042 % to ensure 0 is a ppar grid point
0043 np = ceil(pmax/dp);
0044 ppar = linspace(-pmax,pmax,2*np + 1);
0045 % pperp = [0:0.05:9.95,linspace(10,pmax,np + 1)];
0046 pperp = linspace(0,pmax,np + 1);
0047 dp = pperp(2) - pperp(1);
0048 %
0049 [ppar1,pperp1] = meshgrid(ppar,pperp);
0050 f1 = griddata(ppar00,pperp00,f00,ppar1,pperp1,'cubic');
0051 %
0052 dppar = dp*ones(size(ppar));
0053 dpperp = dp*ones(size(pperp));
0054 
0055  
0056  
0057  
0058  
0059  
0060