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

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

//#include "integ.h"

#include "abscalctool.h"

#include "lightdipole.h"
#include "sigcalctools.h"


double AbsCalcTool::compPixelQD(double theta,double phi) 
{	
//	cout<<"theta: "<<theta<<"\tphi: "<<phi<<endl;
	UnitVector VP(theta,phi);
	UnitVector VYbidon=VP.VperpPhi();	  
// Compute unit vector perpendicular to Vpoin at same theta
	return compPixel(VP,VYbidon);
}

double AbsCalcTool::NormKelvinRayleighJeans()
{	
 	double tempeCNoir=10000.;	
 					// Kelvin
 	double CutFreq=1.380662e-23*tempeCNoir/6.626176e-34/5.;
 	if(FreqMax>1.380662e-23*tempeCNoir/6.626176e-34/5.) 
 	{	cerr<< "RaleighJeans approximation is not valid for this frequency"<<endl;
 		cerr<< "Frequency: "<< FreqMax<<" in SigCalcTool::NormRayleighJeans"<<endl;
 	}
 	
 	LightBlackBody CorpsNoir(tempeCNoir, RAngComp);
 	SigCalcTool ToolRJ(&CorpsNoir,pLobe,pFilter);
 	double puissNorm = ToolRJ.compPixelQD(M_PI/2.,M_PI);	// Un pixel au hasard
 	return tempeCNoir/puissNorm;				// Kelvin RaleighJeans/(Watt/m2)
 
}
 
double AbsCalcTool::NormKelvinCMB()
{	
 	double deltatempeCNoir=1.;					// Kelvin
 	LightNormTCMB DeltaCorpsNoir(deltatempeCNoir, RAngComp);
 	SigCalcTool ToolDeltaCMB(&DeltaCorpsNoir,pLobe,pFilter);
 	double puissNorm = ToolDeltaCMB.compPixelQD(M_PI/2.,M_PI);	// Un pixel au hasard
 	return deltatempeCNoir/puissNorm;			// KelvinCMB/(Watt/m2)
}



/*	// Pour compilation sur cxx boeuf!
double AbsCalcTool::max(double a, double b) const{
  if(a>b) return a;
  else return b;
}

double AbsCalcTool::min(double a, double b) const{
  if(a<b) return a;
  else return b;
}
*/
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, AbsCalcTool& 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.compPixelQD(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, AbsCalcTool& 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.compPixelQD(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, AbsCalcTool& Tool);
template void addToSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
template void compSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
template void compSkyMap(PixelMap<float>& Map, AbsCalcTool& 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);