#include <iostream>
#include "ptsrcinbandcalctools.h"
#include "lightptsrclevsinband.h"
#include "spherehealpix.h"

PtSrcInBandCalTools::PtSrcInBandCalTools(LightPtSrcLevSInBand* pLightSource,MeanFreqLobe* pLob, LevSPanckBand thisBand)
	:pLightSrc(pLightSource), Band(thisBand)
{	pLobe=pLob;	
	FreqMax=pLobe->maxFreq();
	FreqMin=pLobe->minFreq();
	pFilter=pLob->pRespShape;
	RAngComp =  pLobe->lobeResol();
}

double PtSrcInBandCalTools::compPixel(UnitVector& VP, UnitVector& VY)
{	double OuvAng=pLobe->AngleMax();
	double thetaPointe=VP.Theta();
	double PhiPointe=VP.Phi();

	double ThetaMin=thetaPointe-OuvAng; if(ThetaMin<0.) ThetaMin=0.;
	double ThetaMax=thetaPointe+OuvAng;	if(ThetaMax>M_PI) ThetaMax=M_PI;
	
/*
	cout<< "test Passage"<<endl;
	cout<<VP.Theta()<<'\t'<<VP.Phi()<<'\t'<<VY.Theta()<<'\t'<<VY.Phi()<<endl;
*/	
	double deltaPhi;
	{ 	// Compute limits in Phi
		double deltaPhi1;
		if(ThetaMin!=0.) deltaPhi1=OuvAng/sin(ThetaMin);
		else	deltaPhi1=M_PI;
		double deltaPhi2;
		if(ThetaMax!=0.) deltaPhi2=OuvAng/sin(ThetaMax);
		else deltaPhi2=M_PI;
		
		// Compute deltaPhi
		if(deltaPhi1>deltaPhi2) deltaPhi=deltaPhi1;
		else deltaPhi=deltaPhi2;
	}
	double PhiMin=PhiPointe-deltaPhi;
	double PhiMax=PhiPointe+deltaPhi;

/*	// Useless!
	if((PhiMin<-M_PI)||(PhiMax>M_PI))
	{	PhiMin=-M_PI;
		PhiMax=M_PI;
	}
*/
	double PowerTot=0.;
	
 	const UnitVector VecBidon(ThetaMin,0.);
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterlow=pLightSrc->MapOfPtSrc.lower_bound(VecBidon);
	const UnitVector VecBidon2(ThetaMax,0.);
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterhigh=pLightSrc->MapOfPtSrc.upper_bound(VecBidon2);

	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
	double PhiSrc;
		
//	long compteur=0;
//	long compteur2=0;

	for(iter=iterlow; iter!=iterhigh; iter++)
	{	PhiSrc=(*iter).first.Phi();
		if( (PhiSrc>PhiMin) && (PhiSrc<PhiMax) )
		{ 	double poid=pLobe->weigthAmpl((*iter).first,VP,VY);
			PowerTot+=(*iter).second.getPower(Band)*poid;
//			compteur++;
//			if (poid!=0.) compteur2++;
		}
	}
	return PowerTot;
}

double PtSrcInBandCalTools::CalcLobeSize(double frequency)
	// I KNOW! Awfully heavy. But It's the cost of generality.
	// This function should be called only once!
{	double AngMax=pLobe->AngleMax();

	// On cherche la resolution healpix correspondante
	double x= 0.5*M_PI/RAngComp;
	long nlat=1;
	while (x>1.)
	{x/=2.;	nlat*=2;}
	
	// On reserve une sphere pour en connaitres les coordonnŽes
	SphereHEALPix<uint_2> dummSphere(nlat);
	
	double SolidAngPixel = dummSphere.PixSolAngle();
	long Nbtranch = dummSphere.NbThetaSlices()*AngMax/M_PI;
	double thetaStep=M_PI/Nbtranch+1.;
	
	TVector<uint_2> VecValues;
	TVector<int_4> VecIndice;
	r_8 dumTheta;
	r_8 dumPhi;
	
	double AngSolid=0.;
	UnitVector VZ(0.,0.,1.);
	
	for(double ThetaAnn=0.5*thetaStep;ThetaAnn<AngMax; ThetaAnn+=thetaStep)
	{	// Obtient l'indice d'un pixel sur l'anneau theta=ThetaAnn
		int_4 index=dummSphere.PixIndexSph(ThetaAnn,0.);
		dummSphere.GetThetaSlice(index,dumTheta,dumPhi,VecIndice,VecValues);
		
		double theta,phi;
		// Somme la valeur du lobe sur les pointe d'un anneau
		for(long index=0; index<VecIndice.Size(); index++) 
		{	int_4 ind=VecIndice(index);
			dummSphere.PixThetaPhi(ind,theta,phi);
			UnitVector VInteg(theta, phi);
			UnitVector VY=VZ^VInteg;
			
			AngSolid+= pLobe->weigthAmpl(VInteg,VZ,VY)*SolidAngPixel;
		}
	}
	cout<<AngSolid<<endl;
	
	return AngSolid;
}

/*
double PtSrcInBandCalTools::powerInteg(double ThetaMin, double ThetaMax, 
				double PhiMin, double PhiMax, UnitVector& VP, UnitVector& VY, bool InvLog)
{	
	const UnitVector VecBidon(ThetaMin,0.);
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterlow=pLightSrc->MapOfPtSrc.lower_bound(VecBidon);
	const UnitVector VecBidon2(ThetaMax,0.);
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterhigh=pLightSrc->MapOfPtSrc.upper_bound(VecBidon2);

	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
	double Power=0.;
	double PhiSrc;
	bool test;
		
//	long compteur=0;
//	long compteur2=0;

	for(iter=iterlow; iter!=iterhigh; iter++)
	{	PhiSrc=(*iter).first.Phi();
		test=(PhiSrc>PhiMin) && (PhiSrc<PhiMax);
		if (InvLog) test= (!test);	
		if(test)
		{ 	double poid=pLobe->weigthAmpl((*iter).first,VP,VY);
			Power+=(*iter).second.getPower(Band)*poid;
//			compteur++;
//			if (poid!=0.) compteur2++;
		}
	}
	return Power;
}
*/