libclhep-dev 2.1.2.3-1 / usr / include / CLHEP / GenericFunctions / SphericalHarmonicFit.icc

// -*- C++ -*-
// $Id: 
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_legendre.h>
#include <complex>
#include <cstdlib>
#include <stdexcept>
#include "CLHEP/GenericFunctions/ClebschGordanCoefficientSet.hh"
namespace Genfun {

FUNCTION_OBJECT_IMP(SphericalHarmonicFit)

class SphericalHarmonicFit::Clockwork {
  
public:
  
  Clockwork(unsigned int LMAX):LMAX(LMAX),coefficientsA(LMAX),coefficientsASq(2*LMAX) {}

  struct MStruct {
    unsigned int M;
    Genfun::Parameter *fractionAbsMOrHigher;
    Genfun::Parameter *fractionMPositive;
    Genfun::Parameter *phaseMPlus;
    Genfun::Parameter *phaseMMinus;
  };
  
  struct LStruct {
    unsigned  int      L;
    Genfun::Parameter *fractionLOrHigher;
    Genfun::Parameter *phaseLM0;
    std::vector<MStruct> mstruct;
  };

  std::vector<LStruct>   lstruct;
  const unsigned   int   LMAX;


  SphericalHarmonicCoefficientSet coefficientsA;
  SphericalHarmonicCoefficientSet coefficientsASq;
  ClebschGordanCoefficientSet     ClebschGordan;

  void recomputeCoefficients() {

    // Note, the calling sequence of the GSL Special Function forces us to 
    // transpose Plm from its "natural" order.. It is addressed as P[m][l].
    
    //  double ampSq=0.0;
    
    std::complex<double> I(0,1.0);
    double f=1.0;
    double fThisSum=0.0;
    for (unsigned int l=0;l<=LMAX;l++) {
      
      // lStructThis is zero if l==0;
      // lStructNext is zero if l==LMAX;
      const LStruct *lStructThis= (l==0    ? NULL: & lstruct[l-1]);
      const LStruct *lStructNext= (l==LMAX ? NULL: & lstruct[l]);
      double fHigher = lStructNext ? lStructNext->fractionLOrHigher->getValue() : NULL;
      double fThis   = f*(1-fHigher);
      fThisSum+=fThis;
      
      double g=1.0;
      double gThisSum=0.0;
      for (int m=0;m<=int(l);m++) {
	
	// mStructThis is zero if m==0;
	// mStructNext is zero if m==l;
	const MStruct *mStructThis= ((m==0 || !lStructThis) ? NULL: & lStructThis->mstruct[m-1]);
	const MStruct *mStructNext= (m==int(l) ? NULL: & lStructThis->mstruct[m]);
	double gHigher = mStructNext ? mStructNext->fractionAbsMOrHigher->getValue() : NULL;
	double gThis   = g*(1-gHigher);
	gThisSum+=gThis;
	
	if (fThis<0) {
	  std::cout << "L-fraction correction" << fThis << "-->0" << std::endl;
	  fThis=0.0;
	}
	if (gThis<0) {
	  std::cout << "M-fraction correction" << gThis << "-->0" << std::endl;
	  gThis=0.0;
	}
	double px=0.0; // phase
	if (m==0) {
	  if (lStructThis) {
	    double amplitude = sqrt(fThis*gThis);
	    px = lStructThis->phaseLM0->getValue();;
	    coefficientsA(l,m)=exp(I*px)*amplitude;
	  }
	  // L=0 occurs here:
	  else {
	    double amplitude = sqrt(fThis*gThis);
	    coefficientsA(l,m)=exp(I*px)*amplitude;
	  }
	}
	// Split it between positive and negative:
	else {
	  {
	    double amplitude = sqrt(fThis*gThis*mStructThis->fractionMPositive->getValue());
	    px = mStructThis->phaseMPlus->getValue();;
	    coefficientsA(l,m)=exp(I*px)*amplitude;
	  }
	  {
	    double amplitude = sqrt(fThis*gThis*(1-mStructThis->fractionMPositive->getValue()));
	    px = mStructThis->phaseMMinus->getValue();;
	    coefficientsA(l,-m)=exp(I*px)*amplitude;
	  }
	}
	g*=gHigher;
      }
      f*=fHigher;
    }
  }
};


inline
SphericalHarmonicFit::SphericalHarmonicFit(unsigned int LMAX):
  c(new Clockwork(LMAX))
{
  for (unsigned int l=1;l<=LMAX;l++) {
    Clockwork::LStruct lstruct;
    lstruct.L=l;
    {
      std::ostringstream stream; 
      stream << "Fraction L>=" << l;
      lstruct.fractionLOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
    }
    {
      std::ostringstream stream; 
      stream << "Phase L=" << l << "; M=0";
      lstruct.phaseLM0= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
    }
    for (unsigned int m=1;m<=l;m++) {
      Clockwork::MStruct mstruct;
      mstruct.M=m;
      {
	std::ostringstream stream; 
	stream << "Fraction L= " << l << "; |M| >=" << m;
	mstruct.fractionAbsMOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
      }
      {
	std::ostringstream stream; 
	stream << "Fraction L=" << l << "; M=+" << m ;
	mstruct.fractionMPositive= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
      }
      {
	std::ostringstream stream; 
	stream << "Phase L=" << l << "; M=+" << m ;
	mstruct.phaseMPlus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
      }
      {
	std::ostringstream stream; 
	stream << "Phase L=" << l << "; M=-" << m ;
	mstruct.phaseMMinus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
      }
      
      lstruct.mstruct.push_back(mstruct);
    }
    c->lstruct.push_back(lstruct);
  }
}


inline
SphericalHarmonicFit::~SphericalHarmonicFit() {
  for (unsigned int i=0;i<c->lstruct.size();i++) {
    delete c->lstruct[i].fractionLOrHigher;
    delete c->lstruct[i].phaseLM0;
    for (unsigned int j=0;j<c->lstruct[i].mstruct.size();j++) {
      delete c->lstruct[i].mstruct[j].fractionAbsMOrHigher;
      delete c->lstruct[i].mstruct[j].fractionMPositive;
      delete c->lstruct[i].mstruct[j].phaseMPlus;
      delete c->lstruct[i].mstruct[j].phaseMMinus;
    }
  }
  delete c;
}

  inline
  SphericalHarmonicFit::SphericalHarmonicFit(const SphericalHarmonicFit & right):
    AbsFunction(),
    c(new Clockwork(right.c->LMAX))
  {
    for (unsigned int i=0;i<right.c->lstruct.size();i++) {
      Clockwork::LStruct lstruct;
      lstruct.L= right.c->lstruct[i].L;
      lstruct.fractionLOrHigher = new Parameter(*right.c->lstruct[i].fractionLOrHigher);
      lstruct.phaseLM0          = new Parameter(*right.c->lstruct[i].phaseLM0);
      for (unsigned int j=0;j<right.c->lstruct[i].mstruct.size();j++) {
	Clockwork::MStruct mstruct;
	mstruct.M=right.c->lstruct[i].mstruct[j].M;
	mstruct.fractionAbsMOrHigher=new Parameter(*right.c->lstruct[i].mstruct[j].fractionAbsMOrHigher);
	mstruct.fractionMPositive   =new Parameter(*right.c->lstruct[i].mstruct[j].fractionMPositive);
	mstruct.phaseMPlus          =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMPlus);
	mstruct.phaseMMinus         =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMMinus);
	lstruct.mstruct.push_back(mstruct);
      }
      c->lstruct.push_back(lstruct);
    }
  }
  
  inline
  double SphericalHarmonicFit::operator() (double ) const {
    throw std::runtime_error("Dimensionality error in SphericalHarmonicFit");
    return 0;
  }
  
  inline
  double SphericalHarmonicFit::operator() (const Argument & a ) const {
    unsigned int LMAX=c->LMAX;
    double x = a[0];
    double phi=a[1];

    // Note, the calling sequence of the GSL Special Function forces us to 
    // transpose Plm from its "natural" order.. It is addressed as P[m][l].

    //double Plm[LMAX+1][LMAX+1];  
    std::vector< std::vector<double> > Plm(LMAX+1);
    for (int m=0;m<=int(LMAX);m++) {
      Plm[m].resize(LMAX+1);
      gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
    }

    c->recomputeCoefficients();
    std::complex<double> P=0.0;
    std::complex<double> I(0,1.0);
    for (unsigned int l=0;l<=LMAX;l++) {
      for (int m=0;m<=int(l);m++) {
	{
	  int LP=l-abs(m);
	  double Pn= Plm[abs(m)][LP];
	  
	  if (!finite(Pn)) return 0.0;

	  // Once for positive m (in all cases):
	  P+=(c->coefficientsA(l,m)*Pn*exp(I*(m*phi)));
	  // Once for negative m (skip if m==0);
	  if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficientsA(l,-m)*Pn*exp(-I*(m*phi)));
	}
      }
    }
    
    double retVal=std::norm(P);
    if (!finite(retVal)) {
      return 0.0;
    }
    return retVal;
  }
  
  inline
  unsigned int SphericalHarmonicFit::lMax() const {
    return c->LMAX;
  }
  
  inline
  unsigned int SphericalHarmonicFit::numComponents() const {
    return (c->LMAX+1)*(c->LMAX+1)-1;
  }
  
  // The fraction of Amplitude sq which is L or higher
  inline
  Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L){
    return c->lstruct[L-1].fractionLOrHigher;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L) const {
    return c->lstruct[L-1].fractionLOrHigher;
  }

  // The phase of coefficient L, M=0;
  inline
  Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L){
    return c->lstruct[L-1].phaseLM0;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L) const{
    return c->lstruct[L-1].phaseLM0;
  }
  
  // The fraction of amplitude sq which is L which is +- M OR HIGHER
  inline
  Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
  }
  

  // The fraction of amplitude sq which is +- M, which is positive
  inline
  Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
  }

  inline
  const Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
  }


  // The phase of the positive M coefficient
  inline
  Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
  }
  
  
  // The phase of the negative M coefficient
  inline
  Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
  }
  
  inline
  const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsA() const {
    c->recomputeCoefficients();
    return c->coefficientsA;
  }

  inline
  const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsASq() const{
    c->recomputeCoefficients();
    unsigned int LMAX=c->coefficientsA.getLMax();
    for (unsigned int L=0;L<=2*LMAX;L++) {
      for (int M=-L; M<=int(L); M++) {
	c->coefficientsASq(L,M)=0.0;
	for (unsigned int l1=0;l1<=LMAX;l1++) {
	  for (unsigned int l2=0;l2<=LMAX;l2++) {
	    for (int m1=-l1;m1<=int(l1);m1++) {
	      for (int m2=-l2;m2<=int(l2);m2++) {
		  if (m1-m2==M) {
		    if (((l1+l2) >= L) && abs(l1-l2) <=int(L))  {
		    c->coefficientsASq(L,M) += (c->coefficientsA(l1,m1)*
						conj(c->coefficientsA(l2,m2))*
						(m2%2 ? -1.0:1.0) * 
						sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
						c->ClebschGordan(l1,l2,0,0,L,0)*c->ClebschGordan(l1,l2,m1,-m2,L,M));
		    }
		  }
	      }
	    }
	  }
	}
      }
    }
    return c->coefficientsASq;
  }

  inline
  void SphericalHarmonicFit::recomputeCoefficients() const {
    c->recomputeCoefficients();
  }

} // end namespace Genfun