// Dominique YVON, CEA/DAPNIA/SPP 02/2000

#include <math.h>
#include <iostream>
#include <fstream>
#ifdef __MWERKS__
   #include "mwerksmath.h"
   #include "unixmac.h"
#endif

#include "sigcalctools.h"
//#include "integ.h"
 
template <class T> void addToSkyMap(PixelMap<T>& Map, SigCalcTool& Tool) 
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  cout<<"Nbre de pIxel a calculer dans addToSkyMap: "<<PixelNumber<<endl;
  T temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
    temp=(T) Tool.compPixel(theta,phi); 
//    if (temp<1.e-35) temp=0.;		// BUGG XXXXXX
    								// pour vivre avec les pb de fichiers fits D.Y. 
     Map(k)+=temp;
     if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
  }
   return;
}

template <class T> void compSkyMap(PixelMap<T>& Map, SigCalcTool& Tool)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
  T temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
	temp= (T) Tool.compPixel(theta,phi);
//    if (temp<1.e-35) temp=0.;
				// BUGG XXXXXX
     			// pour vivre avec les pb de fichiers fits D.Y.
	//if (out!=0.) cout<< out <<'\t'<< temp<<endl;		 
    Map(k)= temp;
    if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
  }
   return;
}

template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
   // temp=(T1) Map2.PixValSph(theta,phi);
    if(Map2.ContainsSph(theta,phi)) Map(k) += (T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
 //   temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
    if(Map2.ContainsSph(theta,phi)) Map(k) = Map(k)-(T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
     Map.PixThetaPhi(k,theta,phi);
 //   temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
     if(Map2.ContainsSph(theta,phi)) 
    	 if (Map2.PixValSph(theta,phi)!=0.) 
    		 Map(k) = Map(k)/(T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map)
{
	long PixelNumber= Map.NbPixels();
	//T temp;
	for(long k=0; k<PixelNumber; k++) 
	{	
//		temp= (T) ( scalefactor*Map2(k) );
		Map(k) = (T) ( scalefactor*Map(k) );	
	}
}


template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min, 
	double& Max, double& Moy, double& sigma)
{	
	long PixelNumber= Map.NbPixels();
	double val=0.;
	double variance=0;
	Moy=0.;
	Min=1.e36;
	Max=-1.e-36;
	sigma=0.;
	
	for(long k=0; k<PixelNumber; k++) 
	{	val=Map(k); // cout<<k<<'\t'<<val<<'\n';
	    if(val>Max) Max=val;
		if(val<Min) Min=val;
		Moy+=val;	
	}
	
	Moy/=PixelNumber;
	for(long k=0; k<PixelNumber; k++) variance+=(val-Moy)*(val-Moy);
	variance/=PixelNumber;
	sigma=sqrt(variance);
	
	if(Min==Max) 
	{	cerr<< "bizarre, une carte vide ou uniforme"<<endl;
		return 1;
	}
	return 0;
}

// On instancie les plus courants
template void addToSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
template void addToSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);
template void compSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
template void compSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);	

// Idem pour les outils de cartes

template void scaleMap(double scalefactor, PixelMap<double>& Map2);
template void scaleMap(double scalefactor, PixelMap<float>& Map2);

template int MinMaxSigMap(PixelMap<float>& Map, double& Min, 
	double& Max, double& Moy, double& sigma);
template int MinMaxSigMap(PixelMap<double>& Map, double& Min, 
	double& Max, double& Moy, double& sigma);


template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
template  void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);

template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);

template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);

template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);