Non-uniform discrete velocity

DV_NU.py

# coding: utf-8
# Generate non-uniform discrete velocity using a power law distribution
import numpy as np
    
def nu_dv(lmx, MX, xiMax):
    """
    Generate a non-uniform points with a power law distribution
    
    Parameter
    ---------
    lmx : int
        The order of the power law, can only be even
    MX  : int
        Number of the points
    xiMax : float
        The maximum points
    Return
    ------
    xi : array_like float
    wei : array_like float
    
    Examples
    --------
    >>> xi, wei = nu_dv(3, 20, 5.0)
    """
    sx = xiMax/((MX-1)**lmx)
    xi = sx*(np.arange(MX)*2+1-MX)**3
    wei = lmx*xiMax*pow((-MX + 1 + 2*np.arange(MX)), lmx-1)/pow(MX-1, lmx)*2
    return xi, wei

if __name__ == '__main__':
    def f(x):
        ''' helper Maxwell equlibrium function '''
        return 1.0/np.sqrt(2.0*np.pi)*np.exp(-x**2/2.0)
    
    xi, wei = nu_dv(3, 32, 5.0)
    print "1st order moment calculated: ", sum(wei*f(xi))
    print "2st order moment calculated: ", sum(wei*f(xi)*xi)