#include <fstream>
#include "unitvector.h"
#include "lightptsrclevsinband.h"


LightPtSrcLevSInBand::LightPtSrcLevSInBand(char FileOfPtSrc[])
{	// le fichier de pointsource ordonnŽ en Theta
	sprintf(Name, "Sources Ponctuelles du Level S");
	readLevSPtSrc(FileOfPtSrc);
}

void LightPtSrcLevSInBand::readLevSPtSrc(char PtSrcfile[])
{	char Line[200];
	
	double theta, phi; 
	float f30, f44, f70, f100, f143, f217, f352, f545, f857;
	LevSPtSrcData UneSource;
	UnitVector UnCoord;
	
	ifstream HFI(PtSrcfile);
	if (!HFI)
	{	cout<< "Erreur a l'ouverture du fichier de donnees :"<<PtSrcfile<< endl;
		exit(-1);
	}
	
	long compteur=0;
	long compteur2=0;
	while (HFI)
	{	// Data in file in radian and Jansky
	//	cout<<compteur<<'\t'<<HFI.fail()<<'\n'<<endl;;
		HFI>>theta;
		HFI>>phi;
		HFI>>f30;
		HFI>>f44;
		HFI>>f70;
		HFI>>f100;
		HFI>>f143;
		HFI>>f217;
		HFI>>f352;
		HFI>>f545;
		HFI>>f857;	
		
		compteur++;
		if(HFI.fail())
		{	HFI.clear();
			HFI.getline(Line, 199);	
			if(HFI.eof())		// Fin de lecture du fichier
			{	cerr<< "Fin de lecture du fichier: "<<PtSrcfile<<endl;
				cerr<< "Point source number : "<<MapOfPtSrc.size()<<endl;
				HFI.close();
				break;
			}
			else cerr<<"erreur lecture ligne: "<<Line<<endl;
		}
		else
		{	// cout<<theta<<'\t'<<phi<<'\t'<<f30<<'\t'<<f44<<'\t'<<f70<<'\t'<<f100;
			// cout<<'\t'<<f143<<'\t'<<f217<<'\t'<<f352<<'\t'<<f545<<'\t'<<f857<<endl<<endl;
			
			// In Memory Power is sored in Watt, angles in Radian
			UnitVector UnCoord(theta, phi);
			UneSource.flux30G=f30*1.e-26;
			UneSource.flux44G=f44*1.e-26;
			UneSource.flux70G=f70*1.e-26;
			UneSource.flux100G=f100*1.e-26;
			UneSource.flux143G=f143*1.e-26;
			UneSource.flux217G=f217*1.e-26;
			UneSource.flux352G=f352*1.e-26;
			UneSource.flux545G=f545*1.e-26;
			UneSource.flux857G=f857*1.e-26;
			
			const pair<const UnitVector,LevSPtSrcData>UneValeur=make_pair(UnCoord, UneSource);
			
			MapOfPtSrc.insert(MapOfPtSrc.end(),UneValeur);
			compteur2++;
		}
	}
	cerr<< compteur<<'\t'<<compteur2<<endl;
	return;
}

#include "fitsspherehealpix.h"

void LightPtSrcLevSInBand::DumpMap(LevSPanckBand band)
{	// Write a Mollweide projection of the Point Sources Sky
	// power in the LevSPanckBand band
	SphereHEALPix <r_4> Sphere(128);

	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterBegin=MapOfPtSrc.begin();
	multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterEnd=MapOfPtSrc.end();
	
	long indexSphere;
	double theta;
	double phi;
	double power;
	
	for(iter=iterBegin; !(iter==iterEnd); iter++)
	{	theta=(*iter).first.Theta();
		phi=(*iter).first.Phi();
		power=(*iter).second.getPower(band);
		
		indexSphere=Sphere.PixIndexSph(theta,phi);
		Sphere(indexSphere) += power;
	}
	
	FITS_SphereHEALPix <r_4> FSphere(Sphere);
	FSphere.Mollweide_picture_projection("SkyOfPtSrcInBandMollweide.fits");

	return;
}