DistributionFunctions.h

Go to the documentation of this file.

// Copyright András Vukics 2006–2014. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
// -*- C++ -*-
#ifndef CPPQEDCORE_QUANTUMDATA_DISTRIBUTIONFUNCTIONS_H_INCLUDED
#define CPPQEDCORE_QUANTUMDATA_DISTRIBUTIONFUNCTIONS_H_INCLUDED

#include "BlitzArray.h"
#include "ComplexExtensions.h"
#include "MathExtensions.h"

#include "ParsFwd.h"

#include <boost/math/special_functions/factorials.hpp>


namespace quantumdata {


struct ParsFunctionScan {
  
  double
    &fLimitXUL, &fLimitYUL, &fLimitXL, &fLimitXU, &fLimitYL, &fLimitYU, &fStep;
    
  int &fCutoff;
  
  ParsFunctionScan(parameters::ParameterTable& p, const std::string& mod="");

};


namespace details {

double w(size_t n, double r, size_t k);

}


template<typename DensityOperatorFunctor>
double wignerFunction(const DensityOperatorFunctor& rho, double x, double y, size_t truncatedDimension=0)
{
  using namespace mathutils; using boost::math::factorial;

  struct Helper
  {
    static dcomp _(const DensityOperatorFunctor& rho, size_t dim, double r, size_t k)
    {
      dcomp res(0.);
      for (size_t n=0; n<dim-k; ++n) res+=details::w(n,r,k)*rho(n+k)(n);
      return res;
    }
  };
  
  const double r=sqrt(sqr(x)+sqr(y)), phi=(x || y) ? atan2(y,x) : 0.;

  const size_t dim=truncatedDimension && truncatedDimension<rho.getDimension(0) ? truncatedDimension : rho.getDimension(0);

  double res=real(Helper::_(rho,dim,r,0));
  
  for (size_t k=1; k<dim; ++k) res+=2*real(Helper::_(rho,dim,r,k)*exp(-DCOMP_I*(k*phi)));

  return res;
}

template<typename DensityOperator>
double qFunction(const DensityOperator& rho, double x, double y, size_t)
{
  using namespace mathutils;

  const dcomp alpha(x,y);

  dcomp qComplex;
  
  for (size_t m=0; m<rho.getDimension(); ++m) for (size_t n=0; n<rho.getDimension(); ++n)
    qComplex+=coherentElement(n,alpha)*coherentElement(m,conj(alpha))*rho(m)(n);
  
  return exp(-sqrAbs(alpha))*real(qComplex);
}
  

class WignerFunctionKernelOld
{
public:
  typedef DArray<1> Hermites;
  
  WignerFunctionKernelOld(double x, double y, size_t dim);

  dcomp operator()(size_t m, size_t n) const;
  
private:
  const Hermites hermite_m2x_, hermite_2y_;
  
};



template<typename DensityOperatorFunctor>
double wignerFunctionOld(const DensityOperatorFunctor& rho, double x, double y, size_t truncatedDimension=0)
// NEEDS_WORK should refer rho only via some traits class to make the code really generic.
{
  using namespace mathutils; using boost::math::factorial;
  
  const size_t dim=truncatedDimension ? truncatedDimension : rho.getDimension(0);
  
  const WignerFunctionKernelOld wignerFunctionKernel(x,y,dim);
  
  double res=0;

  for (size_t m=0; m<dim; ++m) {
    res+=real(minusOneToThePowerOf(m)/factorial<double>(m)*rho(m,m)*pow(2.*DCOMP_I,-2*m)*wignerFunctionKernel(m,m));
    
    for (size_t n=m+1; n<dim; ++n)
      res+=2.*real(minusOneToThePowerOf(m)/sqrt(factorial<double>(n)*factorial<double>(m))*rho(n,m)*pow(2.*DCOMP_I,-m-n)*wignerFunctionKernel(n,m));
    
  }
  
  return 2./PI*exp(-2*(sqr(x)+sqr(y)))*res;
  
}


template<typename DistributionFunctor,typename DensityOperator>
std::ostream& scanFunction(DistributionFunctor distributionFunctor, const DensityOperator& rho, std::ostream& os, const ParsFunctionScan& pfs)
{
  if (pfs.fLimitXUL) pfs.fLimitXL=-(pfs.fLimitXU=pfs.fLimitXUL);
  if (pfs.fLimitYUL) pfs.fLimitYL=-(pfs.fLimitYU=pfs.fLimitYUL);
  for (double x=pfs.fLimitXL; x<pfs.fLimitXU; x+=pfs.fStep) {
    for (double y=pfs.fLimitYL; y<pfs.fLimitYU; y+=pfs.fStep) os<<x<<"\t"<<y<<"\t"<<distributionFunctor(rho,x,y,pfs.fCutoff)<<std::endl;
    os<<std::endl;
  }
  return os;
}



}

#endif // CPPQEDCORE_QUANTUMDATA_DISTRIBUTIONFUNCTIONS_H_INCLUDED