g_dke_yp

    Chang and Cooper function for Fokker-Planck calculations
    see ref. C.F.F.Karney in Computer Physics Reports 4 (1986) 183-244

    Input:
         
         - x: initial vector (uniform grid term) [n,m]
         - y: initial vector (non-uniform grid term) [n,m]
                (default y = 0, uniform grid)

    Output: 

         - z: Chang and Cooper function [n,m]

Last Update: 02/12/2002

by Y.PEYSSON CEA-DRFC 11/05/1993 <yves.peysson@cea.fr>
revised for MatLab 4.1 (07/09/1994) 
revised for MatLab 5.2 (04/07/2002)

0001 function [z] = g_dke_yp(x,y)
0002 %
0003 %    Chang and Cooper function for Fokker-Planck calculations
0004 %    see ref. C.F.F.Karney in Computer Physics Reports 4 (1986) 183-244
0005 %
0006 %    Input:
0007 %
0008 %         - x: initial vector (uniform grid term) [n,m]
0009 %         - y: initial vector (non-uniform grid term) [n,m]
0010 %                (default y = 0, uniform grid)
0011 %
0012 %    Output:
0013 %
0014 %         - z: Chang and Cooper function [n,m]
0015 %
0016 %Last Update: 02/12/2002
0017 %
0018 %by Y.PEYSSON CEA-DRFC 11/05/1993 <yves.peysson@cea.fr>
0019 %revised for MatLab 4.1 (07/09/1994)
0020 %revised for MatLab 5.2 (04/07/2002)
0021 %
0022 if nargin < 1,
0023      error('Wrong number of input arguments for g_dke_yp');
0024 elseif nargin == 1,
0025       y = zeros(size(x));%Uniform grid
0026 end
0027 %
0028 z = zeros(size(x));
0029 mask = abs(x)>1e-6;
0030 z(mask) = 1.0./x(mask) - 1.0./(exp(x(mask) + y(mask)) - 1);
0031 %
0032 if sum(sum(y)) == 0,
0033     z(~mask) = 0.5 - x(~mask)/12 + (x(~mask)).^3/720;
0034 end
0035 %