#include "TMath.h"
#include "TComplex.h"
#include "Math/SpecFunc.h"

//this function is very slow if loaded in interpreter-mode//
//if you want to make it faster, just let root create a dll via ACLIC (see root user guide to know how to use it//
//see YlmFit.c//

Double_t YlmFunc(Double_t *x, Double_t *par)
{
	const double angle	= 0;
	
	const double c		= TMath::Cos(angle*TMath::DegToRad());
	const double s		= TMath::Sin(angle*TMath::DegToRad());

	const double phi	= 0;
	//check your histogram whether values in x are in deg or rad//
	//convert to rad//
	const double theta	= TMath::Cos(TMath::DegToRad()*x[0]);
	
	//otherwise:
	//const double theta	= x[0];

	//if you want to fit a 2D histo, where the y-values are the the phi values//
	//access to y-values: x[1]//
	//const double phi	= x[1];//

	//Define the Spherical harmonics//
	TComplex Y11 (ROOT::Math::assoc_legendre(1,1,theta)*TMath::Cos(phi), ROOT::Math::assoc_legendre(1,1,theta)*TMath::Sin(phi));
	TComplex Y10 (ROOT::Math::assoc_legendre(1,0,theta), 0);
	TComplex Y1m1(-1.*TComplex::Conjugate(Y11));

	TComplex Y21 (ROOT::Math::assoc_legendre(2,1,theta)*TMath::Cos(phi), ROOT::Math::assoc_legendre(2,1,theta)*TMath::Sin(phi));
	TComplex Y20 (ROOT::Math::assoc_legendre(2,0,theta), 0);
	TComplex Y2m1(-1.*TComplex::Conjugate(Y21));

	TComplex Y31 (ROOT::Math::assoc_legendre(3,1,theta)*TMath::Cos(phi), ROOT::Math::assoc_legendre(3,1,theta)*TMath::Sin(phi));
	TComplex Y30 (ROOT::Math::assoc_legendre(3,0,theta), 0);
	TComplex Y3m1(-1.*TComplex::Conjugate(Y31));

	TComplex Y41 (ROOT::Math::assoc_legendre(4,1,theta)*TMath::Cos(phi), ROOT::Math::assoc_legendre(4,1,theta)*TMath::Sin(phi));
	TComplex Y40 (ROOT::Math::assoc_legendre(4,0,theta), 0);
	TComplex Y4m1(-1.*TComplex::Conjugate(Y41));

	//get the amplitudes from the parameter vector//
	TComplex A00(par[0] ,0.);
	TComplex A10(par[1] ,par[2]);
	TComplex A11(par[3] ,par[4]);
	TComplex A20(par[5] ,par[6]);
	TComplex A21(par[7] ,par[8]);
	TComplex A30(par[9] ,par[10]);
	TComplex A31(par[11],par[12]);
	TComplex A40(par[13],par[14]);
	TComplex A41(par[15],par[16]);

	//Calculate the things//
	TComplex tx	=	A00 +
					A10*Y10 +
					A20*Y20 +
					A30*Y30 +
					A40*Y40;

	TComplex ty	=	A11*(Y11+Y1m1) +
					A21*(Y21+Y2m1) +
					A31*(Y31+Y3m1) +
					A41*(Y41+Y4m1);

	TComplex tz	=	A11*(Y11-Y1m1) +
					A21*(Y21-Y2m1) +
					A31*(Y31-Y3m1) +
					A41*(Y41-Y4m1);

	//calculate the return value//
	double retval = TComplex::Abs(tx*c + ty*s) * TComplex::Abs(tx*c + ty*s);

	return retval;
}