#include <stdlib.h>
#include <iostream>
#include "lobelevshealpix.h"
#include "lightptsrclevsinband.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 4
		// On lit la Sphere Healpix
		FitsInFile File(FileOflobeMap); 
		FITS_SphereHEALPix <r_4> FSphere;
		FSphere.Read(File);
		cerr<<"Lecture des donnŽes du lobe LevelSHealpix"<<endl;
		
		// On stocke les donnees
		pSphere=new SphereHEALPix<r_4> ((SphereHEALPix<r_4>) FSphere );

/*		// 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;
*/

		 
/*	   pSphere= new SphereHEALPix<float>;
	   *pSphere= ( SphereHEALPix<float> ) FSphere;
*/   
	   /*
		SphereHEALPix <float>* pUNSphere=pFSphere->getObj();
		long Nlat=pUNSphere->SizeIndex();
		pSphere=new SphereHEALPix <float>(Nlat);
		*pSphere=(SphereHEALPix <float>) (*pFSphere);
		delete pFSphere;
		*/
}

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.;
	if(cosinusZ==1.) return (*pSphere)(pSphere->PixIndexSph(0.,0.));
	else
	{ 	UnitVector VX=VP^VY;
		double theta = acos (cosinusZ);
		
		Vector3d Vproj=VInteg-cosinusZ*VP;
		double phi;
		
		double norm=Vproj.Norm();
		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;
	}
}