#include "fitsspherehealpix.h"
#include "lobelevshealpix.h"


LobeLevSHealpix::LobeLevSHealpix(char FileOflobeMap[],double angMax_degre, double MeanFreq_Hz)
	:MeanFreqLobe(MeanFreq_Hz*0.85,MeanFreq_Hz*1.15)
	
{		angleMax=angMax_degre/180.*M_PI;			// Radian
		cosanglemax=cos(angleMax);
		
		sprintf(Name,"Lobe LevelS Healpix");

/*		// Version 1		
		pSphere=new SphereHEALPix<r_4>;
		cerr<<"Lecture des donnŽes du lobe LevelSHealpix"<<endl;	
		FITS_SphereHEALPix<float> FSphere(*pSphere);
		FSphere.Read(FileOflobeMap);
*/
/*		// Version 2
		FITS_SphereHEALPix<float> FSphere(FileOflobeMap);	
		cerr<<"Lecture des donnŽes du lobe LevelSHealpix"<<endl;
		pSphere=new SphereHEALPix<r_4>;
		*pSphere= (SphereHEALPix<r_4>) FSphere;
*/		
/*		// Version 3
		FITS_SphereHEALPix<float> FSphere(FileOflobeMap);	
		cerr<<"Lecture des donnŽes du lobe LevelSHealpix"<<endl;
		SphereHEALPix<r_4> & Sphere= (SphereHEALPix<r_4>) FSphere;
		pSphere=&Sphere;
*/
}

LobeLevSHealpix::~LobeLevSHealpix()
{	delete pSphere;
}

double LobeLevSHealpix::weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
  	  const UnitVector& VY) const
{	double cosinusZ=VP*VInteg;
	if(cosinusZ<cosanglemax) return 0.;
	else
	{ 	UnitVector VX=VP^VY;
		double theta = acos (cosinusZ);
		
		Vector3d Vproj=VInteg-cosinusZ*VP;
		double phi;
		
		double norm=Vproj.Norm();
		if(norm==0.)
		{	phi=0.;}
		else
		{	Vproj/=norm;
	  		phi = scangle(VY*Vproj,VX*Vproj); 	// Fonction de geometry.h
	  	}													// double scangle(double sinus, double cosinus)   	
	  // double PixNumb = pSphere->PixIndexSph(theta,phi);
	   double weigthVal= (*pSphere)(pSphere->PixIndexSph(theta,phi));
	   return weigthVal;
	}
}