In [6]:

    

import numpy as np
from scipy import (zeros, array, arange, exp, real, conj, pi,
                   copy, sqrt, meshgrid, size, polyval, fliplr, conjugate,
                   cos, sin)
import scipy.sparse as sp
import scipy.fftpack as ft
import scipy.linalg as la
from scipy.special import genlaguerre
from scipy.special import binom
from scipy.special import sph_harm

from qutip.qobj import Qobj, isket, isoper
from qutip.states import ket2dm
from qutip.parallel import parfor
from qutip.utilities import clebsch
from scipy.misc import factorial

from qutip import *


    

In [7]:

    

def _wigner_clenshaw(rho, xvec, yvec, g=2**0.5):
    """
    Using Clenshaw summation - numerically stable and efficient
    iterative algorithm to evaluate polynomial series.
    
    The Wigner function is calculated as
    :math:`W = e^(-0.5*x^2)/pi * \sum_{L} c_L (2x)^L / sqrt(L!)` where 
    :math:`c_L = \sum_n \\rho_{n,L+n} LL_n^L` where
    :math:`LL_n^L = (-1)^n sqrt(L!n!/(L+n)!) LaguerreL[n,L,x]`
    
    """

    M = np.prod(rho.shape[0])
    X,Y = np.meshgrid(xvec, yvec)
    #A = 0.5 * g * (X + 1.0j * Y)
    A2 = g * (X + 1.0j * Y) #this is A2 = 2*A
    
    rho = rho.full() * (2*np.ones((M,M)) - np.diag(np.ones(M)))
    
    B = np.abs(A2)
    B *= B
    
    #calculation of \sum_{L} c_L (2x)^L / sqrt(L!)
    #using Horner's method
    w0 = rho[0,-1] + A2*0
    L = M-1
    
    while L > 0:
        L -= 1
        #here c_L = _wig_laguerre_val(L, B, np.diag(rho, L))
        w0 = _wig_laguerre_val(L, B, np.diag(rho, L)) + w0 * A2 * (L+1)**-0.5
        
    return w0.real * np.exp(-B*0.5) * (g*g*0.5 / pi)


def _wig_laguerre_val(L, x, c, tensor=True):
    """
    this is evaluation of polynomial series inspired by hermval from numpy.    
    Returns polynomial series
    \sum_n b_n LL_n^L,
    where
    LL_n^L = (-1)^n sqrt(L!n!/(L+n)!) LaguerreL[n,L,x]    
    The evaluation uses Clenshaw recursion
    """

    if len(c) == 1:
        y0 = c[0]
        y1 = 0
    elif len(c) == 2:
        y0 = c[0]
        y1 = c[1]
    else:
        k = len(c)
        y0 = c[-2]
        y1 = c[-1]
        for i in range(3, len(c) + 1):
            k -= 1
            y0,    y1 = c[-i] - y1 * (float((k - 1)*(L + k - 1))/((L+k)*k))**0.5, \
            y0 - y1 * ((L + 2*k -1) - x) * ((L+k)*k)**-0.5
            
    return y0 - y1 * ((L + 1) - x) * (L + 1)**-0.5


    

In [9]:

    

N = 20
xvec = np.linspace(-10, 10, 128)
for i in range(5):
    rho = rand_dm(N)

    Wfft = wigner(rho, xvec, xvec, method='iterative')#'laguerre')
    W = _wigner_clenshaw(rho, xvec, xvec)

    Wdiff = abs(W - Wfft)
    print(np.sum(abs(Wdiff)))


    

    




2.67938058273e-12
1.83599538991e-12
2.21552148199e-12
3.02522354444e-12
1.55810088134e-12

In [ ]: