5#include <gsl/gsl_sf_legendre.h>
18 Clockwork(
unsigned int LMAX):LMAX(LMAX),coefficientsA(LMAX),coefficientsASq(2*LMAX) {}
50 std::complex<double> I(0,1.0);
53 for (
unsigned int l=0;l<=LMAX;l++) {
57 const LStruct *lStructThis= (l==0 ? NULL: & lstruct[l-1]);
58 const LStruct *lStructNext= (l==LMAX ? NULL: & lstruct[l]);
60 double fThis = f*(1-fHigher);
65 for (
int m=0;m<=int(l);m++) {
69 const MStruct *mStructThis= ((m==0 || !lStructThis) ? NULL: & lStructThis->
mstruct[m-1]);
70 const MStruct *mStructNext= (m==int(l) ? NULL: & lStructThis->
mstruct[m]);
72 double gThis = g*(1-gHigher);
76 std::cout <<
"L-fraction correction" << fThis <<
"-->0" << std::endl;
80 std::cout <<
"M-fraction correction" << gThis <<
"-->0" << std::endl;
86 double amplitude = sqrt(fThis*gThis);
88 coefficientsA(l,m)=exp(I*px)*amplitude;
92 double amplitude = sqrt(fThis*gThis);
93 coefficientsA(l,m)=exp(I*px)*amplitude;
101 coefficientsA(l,m)=exp(I*px)*amplitude;
106 coefficientsA(l,-m)=exp(I*px)*amplitude;
121 for (
unsigned int l=1;l<=LMAX;l++) {
125 std::ostringstream stream;
126 stream <<
"Fraction L>=" << l;
130 std::ostringstream stream;
131 stream <<
"Phase L=" << l <<
"; M=0";
134 for (
unsigned int m=1;m<=l;m++) {
138 std::ostringstream stream;
139 stream <<
"Fraction L= " << l <<
"; |M| >=" << m;
143 std::ostringstream stream;
144 stream <<
"Fraction L=" << l <<
"; M=+" << m ;
148 std::ostringstream stream;
149 stream <<
"Phase L=" << l <<
"; M=+" << m ;
153 std::ostringstream stream;
154 stream <<
"Phase L=" << l <<
"; M=-" << m ;
158 lstruct.
mstruct.push_back(mstruct);
167 for (
unsigned int i=0;i<c->
lstruct.size();i++) {
168 delete c->
lstruct[i].fractionLOrHigher;
170 for (
unsigned int j=0;j<c->
lstruct[i].mstruct.size();j++) {
171 delete c->
lstruct[i].mstruct[j].fractionAbsMOrHigher;
172 delete c->
lstruct[i].mstruct[j].fractionMPositive;
173 delete c->
lstruct[i].mstruct[j].phaseMPlus;
174 delete c->
lstruct[i].mstruct[j].phaseMMinus;
185 for (
unsigned int i=0;i<right.c->
lstruct.size();i++) {
190 for (
unsigned int j=0;j<right.c->
lstruct[i].mstruct.size();j++) {
192 mstruct.
M=right.c->
lstruct[i].mstruct[j].M;
197 lstruct.
mstruct.push_back(mstruct);
205 throw std::runtime_error(
"Dimensionality error in SphericalHarmonicFit");
211 unsigned int LMAX=c->
LMAX;
219 std::vector< std::vector<double> > Plm(LMAX+1);
220 for (
int m=0;m<=int(LMAX);m++) {
221 Plm[m].resize(LMAX+1);
222 gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
226 std::complex<double> P=0.0;
227 std::complex<double> I(0,1.0);
228 for (
unsigned int l=0;l<=LMAX;l++) {
229 for (
int m=0;m<=int(l);m++) {
232 double Pn= Plm[abs(m)][LP];
234 if (!finite(Pn))
return 0.0;
239 if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->
coefficientsA(l,-m)*Pn*exp(-I*(m*phi)));
244 double retVal=std::norm(P);
245 if (!finite(retVal)) {
264 return c->
lstruct[L-1].fractionLOrHigher;
269 return c->
lstruct[L-1].fractionLOrHigher;
275 return c->
lstruct[L-1].phaseLM0;
280 return c->
lstruct[L-1].phaseLM0;
286 return c->
lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
291 return c->
lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
298 return c->
lstruct[L-1].mstruct[M-1].fractionMPositive;
303 return c->
lstruct[L-1].mstruct[M-1].fractionMPositive;
310 return c->
lstruct[L-1].mstruct[M-1].phaseMPlus;
315 return c->
lstruct[L-1].mstruct[M-1].phaseMPlus;
322 return c->
lstruct[L-1].mstruct[M-1].phaseMMinus;
327 return c->
lstruct[L-1].mstruct[M-1].phaseMMinus;
340 for (
unsigned int L=0;L<=2*LMAX;L++) {
341 for (
int M=-L; M<=int(L); M++) {
343 for (
unsigned int l1=0;l1<=LMAX;l1++) {
344 for (
unsigned int l2=0;l2<=LMAX;l2++) {
345 for (
int m1=-l1;m1<=int(l1);m1++) {
346 for (
int m2=-l2;m2<=int(l2);m2++) {
348 if (((l1+l2) >= L) && abs(l1-l2) <=
int(L)) {
352 sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
353 c->
ClebschGordan(l1,l2,0,0,L,0)*c->
ClebschGordan(l1,l2,m1,-m2,L,M));
#define FUNCTION_OBJECT_IMP(classname)
virtual double getValue() const
unsigned int getLMax() const
std::vector< LStruct > lstruct
ClebschGordanCoefficientSet ClebschGordan
SphericalHarmonicCoefficientSet coefficientsASq
Clockwork(unsigned int LMAX)
void recomputeCoefficients()
SphericalHarmonicCoefficientSet coefficientsA
const SphericalHarmonicCoefficientSet & coefficientsASq() const
virtual ~SphericalHarmonicFit()
Parameter * getFractionAbsMOrHigher(unsigned int L, unsigned int M)
Parameter * getPhaseMMinus(unsigned int L, unsigned int M)
virtual double operator()(double argument) const override
unsigned int numComponents() const
const SphericalHarmonicCoefficientSet & coefficientsA() const
unsigned int lMax() const
SphericalHarmonicFit(unsigned int LMAX)
void recomputeCoefficients() const
Parameter * getPhaseMPlus(unsigned int L, unsigned int M)
Parameter * getFractionMPositive(unsigned int L, unsigned int M)
Parameter * getPhaseLM0(unsigned int L)
Parameter * getFractionLOrHigher(unsigned int L)
Genfun::Parameter * fractionLOrHigher
std::vector< MStruct > mstruct
Genfun::Parameter * phaseLM0
Genfun::Parameter * fractionMPositive
Genfun::Parameter * phaseMMinus
Genfun::Parameter * phaseMPlus
Genfun::Parameter * fractionAbsMOrHigher