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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include <SIOUX.h>

#ifdef __MWERKS__
   #include "mwerksmath.h"
   #include "unixmac.h"
   #include "macenvvariables.h"
#endif

#include "squarefilt.h"
#include "fitsioserver.h"

#include "alllobe.h"
#include "alllightsources.h"
#include "sigcalctools.h"

#define NLatTestMap (256)

static FitsIoServer FitsServer;
int CompareMapResults();
int TestSommeMapsInBand(char file1[], long nlat);

int CompareMapResults()
{
   cout<< "Test du fonctionnement par la comparaison systematique";
   cout<< " Avec les resultats de supertango"<<endl;



#ifndef __MWERKS__
  char* PATHResults=getenv("PATHResults"); 
  char* PATHDataLScr=getenv("PATHDataLScr");   
  char* PATHCarteLobe=getenv("PATHCarteLobe");
  char* PATHSTangoRes=getenv("PATHSTangoRes");
#endif

//    double Min, Max, Moy, sigma;
    char filename [150];
    	
    sprintf(filename,"CarteCieltot_res0256.fits");
    
	TestSommeMapsInBand(filename,NLatTestMap);

    return 0;
}

int TestSommeMapsInBand(char fileTotPower[], long nlat) // Pas fini Ne pas utiliser!
{
	double Min, Max, Moy, sigma;
	char filename[150];
	 // On lit les donnees brutes
    SphereGorski<float> PowerMap(nlat); 	// Pour se mettre a l'unite nW/m2/St
   
    // Definition channel C
    LobeGaussien LobeGaussChanC(28.2/60.,30.e9,50.e09); 
 //   LobeConique ConeLobe(1.,10.e9,3000.e9);			// pour des test simples
    SquareFilter FlatFilter(36.e9,44.e9);			// Hz
    
    // Test lobes
    // On somme les trois contributions intˇgrˇes sur la bande passante
/*    
    // OK pour le synchrotron
    {	LightSynchro Synchro(nlat);
    	SigCalcTool ChanCTool(&Synchro,&LobeGaussChanC,&FlatFilter);
    	addInInBandPowerMap(PowerMap,ChanCTool);		// OK pour le synchrotron
    }

    // On ajoute la poussiere diffuse
    {	LightDiffDust DiffDust(nlat);
    	SigCalcTool ChanCTool(&DiffDust,&LobeGaussChanC,&FlatFilter);
    	addInInBandPowerMap(PowerMap,ChanCTool);	
    }

    // On ajoute le CMB
    {	LightCMBPrim CMB(nlat);
    	SigCalcTool ChanCTool(&CMB,&LobeGaussChanC,&FlatFilter);
    	addInInBandPowerMap(PowerMap,ChanCTool);
    }
    
  
    // On ajoute du bruit 
    {	SphereGorski<float> NoiseChanCMap(nlat);
    	sprintf(filename, "%schannelC_noise_res%04i.fits",PATHSTangoRes,nlat);
    	FitsServer.load(NoiseChanCMap,filename);
    	scaleMap(1.e-9,NoiseChanCMap);		// Pour se mettre W/m2/St
    	addMap(PowerMap,NoiseChanCMap);
    }
   
    // Bilan carte sommee
    cout<<"Carte Totale"<<endl;
    MinMaxSigMap(PowerMap,Min,Max,Moy,sigma);
    cout <<"PowerMap: "<<"Min="<<Min<<" Max="<<Max;
    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;

	FitsServer.save(PowerMap,file1);
*/
/*
    // Test-Validation Ptotale si ca se passe bien
	{ 	
		SphereGorski<float> ObsTotSTango(256);
	    sprintf(filename, "%schannelC_obs_res0256.fits",PATHSTangoRes);
	    FitsServer.load(ObsTotSTango,filename);
	    scaleMap(1.e-9,ObsTotSTango);			// Pour se mettre W/m2/St
	    
	    MinMaxSigMap(ObsTotSTango,Min,Max,Moy,sigma);
	    cout <<"ObsTotSTango W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
	    
	    divMap1WithMap2(PowerMap,ObsTotSTango);
	    MinMaxSigMap(PowerMap,Min,Max,Moy,sigma);
	    cout <<"Test: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl; 
	    
		// On sauve la carte sur disque
		sprintf(filename,"SinusDivCieltot_res0256.fits");
	    FitsServer.sinus_picture_projection(PowerMap,filename);

		// On regarde la carte obtenue
		sprintf(filename,"SinusCieltot_res%04i.fits",nlat);
	    FitsServer.sinus_picture_projection(PowerMap,filename);
	 }
*/

/*
// On veut convoluer avec un lobe gaussien
	{	

		LightSrcMapPowerInband LightMapInCBand(fileTotPower,256,36.e9,44.e9);
			// On utilise une carte deja integree sur la reponse spectrale 
			// Dont les bornes en frequence sont specifiees dans le constructeur
		 	// Pour la forme en coherence avec  ChannelC 
			// On utilse un lobe (large bande)
			
		
			// Voila l'outil
		SigCalcTool ChannelCtool(&LightMapInCBand,&LobeGaussChanC,&FlatFilter);
		
			// La carte
		SphereGorski<float> ChannelCConvolue(nlat); 
		
			// On convolue
		addToSkyMap(ChannelCConvolue, ChannelCtool);
		
			// On sauvegarde
		sprintf(filename,"ChannelConvolue%04i.fits",nlat);
		FitsServer.save(ChannelCConvolue, filename);
		
		MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
	    cout <<"ChannelCConvolue: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;	

			// On regarde la carte obtenue
		sprintf(filename,"SinusCielCconvol%04i.fits",nlat);
	    FitsServer.sinus_picture_projection(ChannelCConvolue,filename);
	}
 */
 
 /*	   
	{		// On Compare
#ifndef __MWERKS__
  char* PATHSTangoRes=getenv("PATHSTangoRes");
#endif
		SphereGorski<float> ObsConvolSTango(nlat);
	    sprintf(filename, "%schannelC_convolved.fits",PATHSTangoRes);
	    FitsServer.load(ObsConvolSTango,filename);
	    	// Unite nW/m2/St
	    	
	    	// On met a la meme echelle
	    LightSrcMapPowerInband LightMapInCBand(fileTotPower,256,36.e9,44.e9);
		SigCalcTool ChannelCtool(&LightMapInCBand,&LobeGaussChanC,&FlatFilter);
		double lobeSize=ChannelCtool.CalcLobeSize(); //Steradian
	    double scaleFactor= 1.e-9 * lobeSize;	// Pour se mettre W/m2
	    scaleMap(scaleFactor,ObsConvolSTango);			

	    MinMaxSigMap(ObsConvolSTango,Min,Max,Moy,sigma);
	    cout <<"ObsConvolSTango W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;	    
		    
	    SphereGorski<float> ChannelCConvolue(nlat);
	    sprintf(filename,"ChannelConvolue%04i.fits",nlat);
	    FitsServer.load(ChannelCConvolue,filename);
	    	         
	    MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
	    cout <<"ChannelCConvolue W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;	    
 		
 			// On divise
	    divMap1WithMap2(ChannelCConvolue,ObsConvolSTango);
	    MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
	    cout <<"Test: "<<"Min="<<Min<<" Max="<<Max;
	    cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
	    
	    	// On sauvegarde la carte
	    sprintf(filename,"SinusConvolDiv%04i.fits",nlat);
	    FitsServer.sinus_picture_projection(ChannelCConvolue,filename);
    }
*/
	return 0;
}



 /*   
    // On lit les templates pour voir
  	sprintf(filename, "%sdust_res0256.fits",PATHDataLScr);
    FitsServer.load(DiffDustMap,filename);
    MinMaxSigMap(DiffDustMap,Min,Max,Moy,sigma);
    flog<<"DiffDustMap: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;

 */
 /*
	sprintf(filename, "%ssync_res0256.fits",PATHDataLScr);
    FitsServer.load(SynchroMap,filename);
    MinMaxSigMap(SynchroMap,Min,Max,Moy,sigma);
    flog<<"SynchroMap: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
 */
 
/*
    sprintf(filename, "%scmb_res0256.fits",PATHDataLScr);
    FitsServer.load(CMBPrimMap,filename);
    MinMaxSigMap(CMBPrimMap,Min,Max,Moy,sigma);
    flog<<"CMBPrimMap: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
 */  
 
 /*
  // Test Module synchrotron
    SphereGorski<float> ChanCSynchro(256);
    SphereGorski<float> MYCSynchro(256);
   
   
   // Sans lobes d'abord
   		// Lecture resultats de supertango
    sprintf(filename, "%schannelC_sync_res0256.fits",PATHSTangoRes);
    FitsServer.load(ChanCSynchro,filename);
    MinMaxSigMap(ChanCSynchro,Min,Max,Moy,sigma);
    flog<<"ChanCSynchro: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
    
    	// Calcul SigPred 
    LightSynchro Synchro(256);
    SigCalcTool ChanCTool(&Synchro,&ConeLobe,&FlatFilter);
    //ChanCTool.SetLightScr(&Synchro);
   
 	addInBandPower(MYCSynchro,ChanCTool);
	scaleMap(1.e9,MYCSynchro);		// Pour se mettre a l'unite nW/m2/St
	MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
    flog<<"MYCSynchro: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
*/    
    // Validation	
/*    //substractMap(MYCSynchro,ChanCSynchro);
    divMap1WithMap2(MYCSynchro,ChanCSynchro);
    MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
    flog<<"Division des cartes: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;    
*/
/*	// test
	divMap1WithMap2(MYCSynchro,SynchroMap);
	MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
    flog<<"MycoeffIntegrale MYCSynchro: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
    
	divMap1WithMap2(ChanCSynchro,SynchroMap);
	MinMaxSigMap(ChanCSynchro,Min,Max,Moy,sigma);
    flog<<"STangoCoeffIntegrale ChanCSynchro: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
  
*/  

/*
      // Test Module CMB
    SphereGorski<float> ChanCCMB(256);
    SphereGorski<float> MYCCMB(256);
    		
  // Sans lobes d'abord: Integrale sur les dependances spectrales
  //degrˇs, freqmin, Freqmax
    LightCMBPrim CMB(256);
    SigCalcTool ChanCTool(&CMB,&ConeLobe,&FlatFilter);
    
  		// Lecture resultats de supertango
    sprintf(filename,"%schannelC_cmb_res0256.fits",PATHSTangoRes);
    FitsServer.load(ChanCCMB,filename);
    MinMaxSigMap(ChanCCMB,Min,Max,Moy,sigma);
    flog<<"ChanCCMB: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
*/   
/*   	// Calcul SigPred   
 	addInBandPower(MYCCMB,ChanCTool);
	scaleMap(1.e9,MYCCMB);		// Pour se mettre a l'unite nW/m2/St
	MinMaxSigMap(MYCCMB,Min,Max,Moy,sigma);
    flog<<"MYCCMB: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
*/   
    	// Validation
    	
/*	// test
	divMap1WithMap2(MYCCMB,CMBPrimMap);
	MinMaxSigMap(MYCCMB,Min,Max,Moy,sigma);
    flog<<" MycoeffIntegrale MYCCMB: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
    
	divMap1WithMap2(ChanCCMB,CMBPrimMap);
	MinMaxSigMap(ChanCCMB,Min,Max,Moy,sigma);
    flog<<"STangoCoeffIntegrale ChanCCMB: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
*/	

/*		
  // Sans lobes d'abord: Integrale sur les dependances spectrales
  //degrˇs, freqmin, Freqmax
    LightDiffDust DiffDust(256);
    SigCalcTool ChanCTool(&DiffDust,&ConeLobe,&FlatFilter);
   
      // Test Module poussi¸re diffuse
    SphereGorski<float> ChanCDiffDust(256);
    SphereGorski<float> MYCDiffDust(256);
    
    
  		// Lecture resultats de supertango
    sprintf(filename,"%schannelC_dust_res0256.fits",PATHSTangoRes);
    FitsServer.load(ChanCDiffDust,filename);
    MinMaxSigMap(ChanCDiffDust,Min,Max,Moy,sigma);
    flog<<"ChanCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
*/  
 
/*   	// Calcul SigPred   
 	addInBandPower(MYCDiffDust,ChanCTool);
	scaleMap(1.e9,MYCDiffDust);		// Pour se mettre a l'unite nW/m2/St
	MinMaxSigMap(MYCDiffDust,Min,Max,Moy,sigma);
    flog<<"MYCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
*/   
    	// Validation
    	
/*	// test
	divMap1WithMap2(MYCDiffDust,DiffDustMap);
	MinMaxSigMap(MYCDiffDust,Min,Max,Moy,sigma);
    flog<<"MycoeffIntegrale MYCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
    
	divMap1WithMap2(ChanCDiffDust,DiffDustMap);
	MinMaxSigMap(ChanCDiffDust,Min,Max,Moy,sigma);
    flog<<"STangoCoeffIntegrale ChanCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
    flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
*/	


int Basictests()
{
// 	MakeTrap99Timeline(); 
//  cerr.setf(ios::scientific);
  cerr<< "C'est parti"<<endl;


 
  SphereGorski<float>* pTestMap;
  pTestMap= new SphereGorski<float>(64);
  SphereGorski<float>& TestMap=(*pTestMap);

#ifndef __MWERKS__
  char* PATHResults=getenv("PATHResults"); 
  char* PATHDataLScr=getenv("PATHDataLScr");   
  char* PATHCarteLobe=getenv("PATHCarteLobe");
#endif

//_________________ Test implementation des lobes ___________________
  LobeGaussien GaussLobe(0.3,75.e9,125.e9);  //degrˇs, freqmin, Freqmax
  LobeConique ConeLobe(1.,75.e9,125.e9);
  Lobe4PiGaussien Lobe4Pi1(1.,75.e9,125.e9);
   
  // Test de la classe Cartelobe
  char lobeFileNameRoot[150];
  sprintf(lobeFileNameRoot, "%sdetector100_1_cf",PATHCarteLobe);
  char LobeFitsFile[150];
  sprintf(LobeFitsFile, "%sCarteLobe100GHz_1_cf.fits",PATHResults);
  
 
//  Test LobeCartoMoyen ___________________________________________ 
//  CarteLobe CarteLobe100GHzCF(lobeFileNameRoot,100.e9);  // Hz
//  LobeCartoMoyen LobeCartoTest(&CarteLobe100GHzCF,75.e9,125.e9);

  // Un filtre bete........
  SquareFilter FlatFilter(85.e9, 115.e9);		// Hz
  
  cout<<"param specifiques a BoloFilterFlat"<<endl;
  cout<<"mean: "<< FlatFilter.meanFreq()<<"/tPeakTransmission: "<<
    FlatFilter.transmission(FlatFilter.meanFreq())<<flush;
  cout<<"Width: "<<FlatFilter.maxFreq()-FlatFilter.minFreq()<<endl;
/* 
  cout<<" Lobe conique:"<<endl;
  cout<<" Nlat gorsky: "<<LobeCartoTest.gorskyNlatFromRes();
  cout<<" Resolution du lobe: "<<LobeCartoTest.lobeResol();
  cout<<" Radian" <<endl; 
 */
  char FileSinus[150];
   
  // Du plus simple au plus compliquˇ
/* 
 // Dipole 
  LightDipole Dipole;
  SigCalcTool DipoleTool(&Dipole,&ConeLobe,&FlatFilter);
//  SigCalcTool DipoleTool(&Dipole,&GaussLobe,&FlatFilter);
  compSkyMap(TestMap,DipoleTool);
  sprintf(FileSinus,"%sMapSinus_Dipole.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);
   
  double MaxDipole=0;
  for (long i=0; i<TestMap.NbPixels(); i++) 
  	if(TestMap(i)>MaxDipole) MaxDipole=TestMap(i);
  cout<<"Maximum du dipole: "<<MaxDipole<<endl;
*/
 /* 
  // Quadrupole 
  LightQuadrupole Quadrupole;
  SigCalcTool QuadTool(&Quadrupole,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,QuadTool);
  sprintf(FileSinus,"%sMapSinus_quadrupole.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);

  // Total des contributions du dipole
  addToSkyMap(TestMap,DipoleTool);
  sprintf(FileSinus, "%sMapSinus_Dip_Quad.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);
*/

/*  // Effet du soleil
  double SunAngsize = 4.833e-3;    // Radian;	
 //  double SunAngsize = 4.833e-2;  	 // Debug
  QuasiPtSources Sun(5800.,SunAngsize,M_PI/2.,M_PI);
  SigCalcTool SunTool(&Sun,&GaussLobe,&FlatFilter);
  //SigCalcTool SunTool(&Sun,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,SunTool);
  
//  sprintf(FileSinus,"%sSphereGorsky_Soleil.fits",PATHResults);
//  FitsServer.save(TestMap, FileSinus);
  
  sprintf(FileSinus,"%sMapSinus_Soleil.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);

  double MaxSun=0.;
  for (long i=0; i<TestMap.NbPixels(); i++) 
  	if(TestMap(i)>MaxSun) MaxSun=TestMap(i);
  cout<<"Maximum du soleil: "<<MaxSun<<endl;
*/
    
/*
  // Fluctuation primordiales
  LightCMBPrim CMBFlucLight(256,1.e-5);
  SigCalcTool CMBGraalTool(&CMBFlucLight,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,CMBGraalTool);
  sprintf(FileSinus,"%sMapSinus_CMBPrimordial.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus); 
*/

/*
  // Carte synchrotron
  LightSynchro Synchro(256);
  SigCalcTool SynchroTool(&Synchro,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,SynchroTool);

  sprintf(FileSinus,"%sMapSinus_Synchro.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);
*/
  
/*
 // Diffuse Dust
  LightDiffDust DiffDust(256);
  SigCalcTool DiffDustTool(&DiffDust,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,DiffDustTool);
  sprintf(FileSinus,"%sMapSinus_DiffDust.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);
  
*/  

/*
  LightGalaxResol Galaxresol(128);
  SigCalcTool GalacResolTool(&Galaxresol,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,GalacResolTool);
  sprintf(FileSinus,"%sMapSinus_GalaxResol.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);
*/
	
/*
  // Fond galactique non resolues
  LightGalaxNoResol GalaxNoresol(128);
  SigCalcTool GalacNoResolTool(&GalaxNoresol,&ConeLobe,&FlatFilter);
  compSkyMap(TestMap,GalacNoResolTool);
  sprintf(FileSinus,"%sMapSinus_GalaxNOResol.fits",PATHResults);
  FitsServer.sinus_picture_projection(TestMap,FileSinus);	
  
*/ 

return 0;
}