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

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

#include "sigcalctools.h"
//#include "integ.h"
 
 int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo) {
  // All coordinates are in Radian.
  /* kmg_euler.c
   *
   * Converts between different coordinate systems. 
   *
   * SELECT   From           To         |   SELECT    From          To
   *    1     RA-Dec (2000)  Galactic   |     4       Ecliptic      RA-Dec    
   *    2     Galactic       RA-DEC     |     5       Ecliptic
Galactic  
   *    3     RA-Dec         Ecliptic   |     6       Galactic
Ecliptic  
   *
   *    Date     Programmer Remarks
   * ----------- ---------- -------
   * 08-Aug-1999 K. Ganga   First version. Copied and modified EULER from 
   *                        the IDL Astrolib. 
   * May 2000, D. Yvon		Change coordinates units to Radians
   */

  /* Local Declarations */
  double       a, b, sb, cb, cbsa;
  long         i;
  const double TWOPI   = 2.0*M_PI;
  const double FOURPI  = 4.0*M_PI;
  const double DEG2RAD = 180.0/M_PI;

  /* J2000 coordinate conversions are based on the following constants
   * eps    =  23.4392911111  Obliquity of the ecliptic
   * alphaG = 192.85948       Right Ascension of Galactic North Pole
   * deltaG =  27.12825       Declination of Galactic North Pole
   * lomega =  32.93192       Galactic longitude of celestial equator
   * alphaE = 180.02322       Ecliptic longitude of Galactic North Pole
   * deltaE =  29.811438523   Ecliptic latitude of Galactic North Pole
   * Eomega =   6.3839743     Galactic longitude of ecliptic equator
   */

  const double psi[6]    = {0.57477043300, 4.9368292465 ,
                            0.00000000000, 0.0000000000 , 
                            0.11142137093, 4.71279419371};
  const double stheta[6] = {0.88998808748,-0.88998808748,
                            0.39777715593,-0.39777715593,
                            0.86766622025,-0.86766622025};    
  const double ctheta[6] = {0.45598377618, 0.45598377618, 
                            0.91748206207, 0.91748206207,
                            0.49714719172, 0.49714719172};
  const double phi[6]    = {4.9368292465 , 0.57477043300,
                            0.0000000000 , 0.00000000000, 
                            4.71279419371, 0.11142137093};

  i    = select - 1;
  a    = ai - phi[i];
  b    = bi;
  sb   = sin(b);
  cb   = cos(b);
  cbsa = cb*sin(a);
  b    = -stheta[i]*cbsa + ctheta[i]*sb;
  b    = ( (b >  1.0) ?  1.0 : b );
  b    = ( (b < -1.0) ? -1.0 : b );
  *bo  = asin(b);

  a    = atan2( ctheta[i] * cbsa + stheta[i] * sb, cb * cos(a) );
  *ao  = fmod(a + psi[i] + FOURPI, TWOPI);

  /* Later */
  return 0;
}


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);