LegendreFit.icc

Go to the documentation of this file.

// -*- C++ -*-
// $Id: 
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_legendre.h>
#include <complex>
namespace Genfun {

FUNCTION_OBJECT_IMP(LegendreFit)

inline
LegendreFit::LegendreFit(unsigned int N):
N(N),coefA(N),coefASq(2*N)
{
  for (unsigned int i=0;i<N;i++) {
    std::ostringstream stream;
    stream << "Fraction " << i;
    fraction.push_back(new Parameter(stream.str(), 0.5, 0.0, 1.0));
  }
  for (unsigned int i=0;i<N;i++) {
    std::ostringstream stream;
    stream << "Phase " << i;
    phase.push_back(new Parameter(stream.str(), M_PI, 0.0, 2.0*M_PI));
  }
}

inline
LegendreFit::~LegendreFit() {
  for (unsigned int i=0;i<N;i++) {
    delete fraction[i];
    delete phase[i];
  }
}

inline
LegendreFit::LegendreFit(const LegendreFit & right):
  AbsFunction(),
  N(right.N),coefA(right.N), coefASq(2*right.N)
{
  for (unsigned int i=0;i<N;i++) {
    fraction.push_back(new Parameter(*right.fraction[i]));
    phase.push_back(new Parameter(*right.phase[i]));
  }
}

inline
double LegendreFit::operator() (double x) const {
  recomputeCoefficients();
  std::vector<double> Pk(N+1);
  gsl_sf_legendre_Pl_array(N, x, &Pk[0]);
  unsigned int n=N;
  std::complex<double> P=0.0;
  std::complex<double> I(0,1.0);
  while (1) {
    if (n==0) {
      double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
      
      P+=coefA(n)*Pn;
      break;
    }
    else {
      double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
      P+=coefA(n)*Pn;
      n--;
    }
  }
  return std::norm(P);
}

inline 
unsigned int LegendreFit::order() const{
  return N;
} 
inline
Parameter *LegendreFit::getFraction(unsigned int i) {
  return fraction[i];
}
inline 
const Parameter *LegendreFit::getFraction(unsigned int i) const{
  return fraction[i];
}
inline
Parameter *LegendreFit::getPhase(unsigned int i) {
  return phase[i];
}
inline 
const Parameter *LegendreFit::getPhase(unsigned int i) const{
  return phase[i];
}
inline
const LegendreCoefficientSet & LegendreFit::coefficientsA() const {
  recomputeCoefficients();
  return coefA;
}

inline
const LegendreCoefficientSet & LegendreFit::coefficientsASq() const {
  recomputeCoefficients();
  unsigned int LMAX=coefA.getLMax();
  for (unsigned int L=0;L<=2*LMAX;L++) {
    coefASq(L)=0.0;
    for (unsigned int l1=0;l1<=LMAX;l1++) {
      for (unsigned int l2=0;l2<=LMAX;l2++) {
        if (((l1+l2) >= L) && abs(l1-l2) <= int(L))  {
          coefASq(L) += (coefA(l1)*
                         conj(coefA(l2))*
                         sqrt((2*l1+1)*(2*l2+1)/4.0)*
                         pow(ClebschGordan(l1,l2,0,0,L,0),2));
        }
      }
    }
  }
  return coefASq;
}

inline 
void LegendreFit::recomputeCoefficients() const {
  unsigned int n=N;
  std::complex<double> P=0.0;
  std::complex<double> I(0,1.0);
  double f=1.0;
  while (1) {
    if (n==0) {
      double fn=1.0;
      coefA(n)=sqrt(f*fn);
      break;
    }
    else {
      double fn=getFraction(n-1)->getValue();
      double px=getPhase(n-1)->getValue();
      coefA(n)=exp(I*px)*sqrt(f*fn);
      f*=(1-fn);
      n--;
    }
  }
}
  

} // end namespace Genfun