Index: /tags/SigPred_V1/SigPredictor/AutresASiGPred.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/AutresASiGPred.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/AutresASiGPred.cc	(revision 799)
@@ -0,0 +1,82 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include "makebolotimeline.h"
+#include "bruit.h"
+//#include "fitsio.h"
+
+#ifdef __MWERKS__
+   #include "macenvvariables.h"
+   #define SetNlatLSource (256)		// version debug
+   #define NbSAMPLE (2000)				//	Pour le debug
+#else
+   # define ERROR_LABEL -32768
+   #define SetNlatLSource (1024)		// version de production CXX
+   #define NbSAMPLE (2684300)
+#endif
+
+
+#include "alllobe.h"
+#include "alllightsources.h"
+#include "sigcalctools.h"
+#include "allfilter.h"
+#include "fitsioserver.h"
+
+int MakeGalaxyMap(double Resolution); 	// Resolution en degres
+
+int MakeGalaxyMap(double Pas)
+{	
+// On passe aux cartes de fonds physiques
+// Initialisations
+   FitsIoServer FitsServer;
+   char filename[150]="";
+     
+//	cerr.setf(ios::scientific);
+  	cerr<< "C'est parti"<<endl;
+  	
+ 	// Des Lobes
+ 	LobeGaussien GLobeTrapChan1(8./60.,75.e9,600.e9);  //degrs, freqmin, Freqmax
+ 	
+#ifndef __MWERKS__
+  char* PathInstrumentData=getenv("PathInstrumentData");
+#endif
+
+ 	// Des Filtres
+ 	sprintf(filename, "%sfiltretrap99chan1.dat",PathInstrumentData);
+ 	TabulatedFilter FiltreTrapChan1(filename); 	
+
+	// Les Processus Physiques
+ 	LightDiffDust DiffDust(SetNlatLSource);
+
+	// Bolo 1_1 soit bolo11.
+ 	sprintf(filename, "ASigPred_MapGalaxie");
+ 	MakeBoloTimeline MapGalaxie(filename,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 		
+ 	// On definit les timelines desirees
+ 	double alphadeg, deltadeg;
+ 	double alpharad, deltarad;
+ 	
+ 	double NbreDePoint=360*40/(Pas*Pas);
+ 	cout<<"Nbre de point a calculer: "<<NbreDePoint<<endl;
+ 	
+ 	long compteur=0;
+ 	for(alphadeg=0; alphadeg<360.; alphadeg+=Pas)
+ 	{	for(deltadeg=-20.; deltadeg<20.; deltadeg+=Pas)
+ 		{	alpharad= alphadeg/180.*M_PI;
+ 			deltarad= deltadeg/180.*M_PI;		
+ 			MapGalaxie.addToStream(alpharad,deltarad);
+ 			
+ 			compteur++;
+ 			if(compteur%5000==0) 
+ 			{	double pourcent=100.*compteur/NbreDePoint;
+ 				cout<<"Nbre points calcules: "<<compteur<<endl;
+ 				cout<<"Travail effectue a :"<<pourcent<<" %"<<endl;
+ 			}
+ 		}	
+ 	}
+	
+
+
+	return 0;
+}
Index: /tags/SigPred_V1/SigPredictor/Lightsource_test.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/Lightsource_test.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/Lightsource_test.cc	(revision 799)
@@ -0,0 +1,48 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+//#include <fstream>
+//#include <SIOUX.h>
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+
+#include "alllobe.h"
+#include "alllightsources.h"
+#include "sigcalctools.h"
+#include "allfilter.h"
+#include "fitsioserver.h"
+
+
+#ifndef __MWERKS__
+  char* PATHResults=getenv("PATHResults"); 
+  char* PATHDataLScr=getenv("PATHDataLScr");   
+  char* PATHCarteLobe=getenv("PATHCarteLobe");
+  char* PathInstrumentData=getenv("PathInstrumentData");
+#endif
+
+
+
+int MakeArcheopsTimelines(char FilePointes[]);
+int MakeGalaxyMap(double Resolution); 		// Resolution en degres
+
+int main()
+{   
+ 	// CompareMapResults();
+ 	// Basictests();
+ 	
+ 	//char FilePointes[32]="bolopoin.ascii";
+	//MakeArcheopsTimelines(FilePointes);
+	
+  	MakeGalaxyMap(1./6); 
+  	
+ 	return 0;
+}
+
+
+
Index: /tags/SigPred_V1/SigPredictor/Tabulatedfilter.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/Tabulatedfilter.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/Tabulatedfilter.cc	(revision 799)
@@ -0,0 +1,75 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "macenvvariables.h"
+#include "tabulatedfilter.h"
+
+NumRecipes  TabulatedFilter::NR;
+static unsigned long LastfreqIndex;
+
+TabulatedFilter::TabulatedFilter(char filename[])
+{	LastfreqIndex=0;
+	FILE* pfile;
+//	char fileCur[150];
+	
+	// On lit le tableau des abscisses
+	pfile=fopen(filename,"r");
+	if (pfile==NULL)
+	{	cerr<< "Erreur a l'ouverture du fichier de donnees :"<<filename<< endl;
+    	exit(-1);
+	}
+	// On lit le nombre d'echantillons
+	double nbreEchantillon; 
+	fscanf(pfile,"%le",&nbreEchantillon);
+	nbSample=(long) nbreEchantillon;
+	Transmi= NR.vector(1,nbSample);
+	frequency=NR.vector(1,nbSample);
+	
+	double freq, filtretrans;
+	for(int i=1; i<=nbSample; i++)
+	{	fscanf(pfile,"%le%le",&freq,&filtretrans);
+		frequency[i]=freq;
+		Transmi[i]=filtretrans;
+	}
+	
+	_numin=1.e99;
+	_numax=0.;
+	for(int i=1; i<=nbSample; i++)
+	{	if(_numin>frequency[i]) _numin=frequency[i];
+		if(_numax<frequency[i]) _numax=frequency[i];
+	}
+}
+
+TabulatedFilter::~TabulatedFilter()	
+{	NR.free_vector(Transmi,1,nbSample);
+	NR.free_vector(frequency,1,nbSample);
+}
+
+double TabulatedFilter::transmission(double freq) const
+{	
+	if((freq<_numin)||(freq>_numax)) return 0.;
+	
+	NR.hunt(frequency,nbSample,freq,&LastfreqIndex);
+	if((LastfreqIndex==0)||(LastfreqIndex==nbSample)) return 0.;
+	else
+	{	double Transinf=Transmi[LastfreqIndex];
+		double Transsup=Transmi[LastfreqIndex+1];
+		double slope=(Transsup-Transinf)/(frequency[LastfreqIndex+1]-frequency[LastfreqIndex]);
+		double retVal=Transinf+slope*(freq-frequency[LastfreqIndex]);
+		return retVal;
+	}
+}	
+ 
+#include <fstream>
+
+void TabulatedFilter::TraceFiltre(ostream& strOut, long nbPoint)
+{	cout<<"On trace un TabulatedFilter"<<endl;
+	strOut<<"Frequence(Hz)"<<'\t'<<"Transmission"<<endl;
+	
+	double pasFreq= (_numax-_numin)/nbPoint;
+	double freq;
+	for(long i=0; i<nbPoint; i++)
+	{	freq=_numin+pasFreq*i;
+		strOut<<freq<<'\t'<<transmission(freq)<<endl;
+	}
+	return;
+}
Index: /tags/SigPred_V1/SigPredictor/TestSources.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/TestSources.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/TestSources.cc	(revision 799)
@@ -0,0 +1,504 @@
+			// 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 intgres 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
+  //degrs, 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
+  //degrs, freqmin, Freqmax
+    LightDiffDust DiffDust(256);
+    SigCalcTool ChanCTool(&DiffDust,&ConeLobe,&FlatFilter);
+   
+      // Test Module poussire 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);  //degrs, 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;
+}
Index: /tags/SigPred_V1/SigPredictor/abslightsource.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/abslightsource.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/abslightsource.cc	(revision 799)
@@ -0,0 +1,50 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+
+#include "abslightsource.h"
+#include "lsrcradspec.h"
+
+
+
+UnitVector AbsLightSource::GetVSrcCenter() {
+  	//Make sense only for QuasiPointSource
+  		cerr<<"Call for GetVSrcCenter() for a non QuasiPointSource: exit"<<endl;
+  		exit(-1);
+  		UnitVector vec;
+  		return vec;
+}
+
+double AbsLightSource::planckSpecDiff(double DTempe, double freq, double Tempe)
+ { // W/m2/st/Hz
+   double h_kT=4.800e-11/Tempe;
+   double prefac= 2*6.626e-34*6.626e-34/(9.e16*1.38e-23*Tempe*Tempe);
+   //
+   return prefac*pow(freq,4)*exp(h_kT*freq)/pow((exp(h_kT*freq)-1.),2)*DTempe;
+}
+
+double AbsLightSource::planckSpectrum(double freq, double Tempe){
+   double h_kT=4.800e-11/Tempe;
+   double prefac=2.*6.626e-34/9.e16;			//2*h/c^2
+   return prefac*freq*freq*freq/(exp(h_kT*freq)-1.);
+      // units (J m-3 Hz-1)
+      // Very usefull spectrum!
+	 	 // W/m2/st/Hz
+}
+
+LSrcRadSpec AbsLightSource::getPowerSpec(double theta,double phi) {
+  	// Return Object including power spectrum versus frequency in direction theta, phi
+  	// (W/m2/Hz/St)  
+  	// Not recommended for intensive use. May Be inefficient or inaccurrate
+  	return LSrcRadSpec(this, theta, phi); 
+}
+
+void AbsLightSource::print(ostream& OutStr) const
+{	OutStr<<"Source de lumiere: "<<Name<<endl;
+	OutStr<<"Resolution: "<< resolution<<endl;
+	return ;
+}
+
Index: /tags/SigPred_V1/SigPredictor/abslightsource.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/abslightsource.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/abslightsource.h	(revision 799)
@@ -0,0 +1,88 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef AbsLightSource_H
+#define AbsLightSource_H
+
+#include <iostream>
+#include "unitvector.h"
+#include "radspec.h"
+#include "lsrcradspec.h"
+#include "specresp.h"
+
+#define T_CMBR (2.726)		 // Kelvin
+
+class AbsLightSource {
+//friend class SigCalcTool;
+public:
+  AbsLightSource() :PtSourceS(false), QPtSource(false), LScrFSep(false), IsLSrcMappedPower (false)
+  {	sprintf(Name, "");
+  }
+  ~AbsLightSource() {	}
+
+  virtual double powSpecDens(double theta,double phi,double freq)=0;
+    // Return power density at coordinates
+    // (W/m2/Hz/St)
+  LSrcRadSpec getPowerSpec(double theta,double phi);
+  	// return Object including power spectrum versus frequency in direction theta, phi
+  	// (W/m2/Hz/St)  
+  	// Not recommended for intensive use. May Be inefficient or inaccurrate..
+  	
+ //virtual double powerDens(double theta,double phi)=0;
+    // Return power Dens at coordinates, integrated over frequency
+    // May Require initialisation
+    // (Watt/m2/St)
+  double LSrcResol() const {
+     return resolution;     // constructeur have to initialise res
+  }
+  
+  inline bool IsPtSourceS() { return PtSourceS;}
+  inline bool IsQPtSrc() { return QPtSource;} 
+  inline bool IsFreqSep() { return LScrFSep;}
+  inline bool IsMappedPowerSrc() {return IsLSrcMappedPower;}
+  	// Optimisation pour calculs avec Lobes separables facon R. Teyssier
+  virtual void print(ostream& OutStr) const;
+/*  
+  virtual void Reset(const SpectralResponse& Filtre) {
+  	cerr<<"you call reset function when not allowed: program exit"<<endl;
+  	exit(-1);
+  	return;
+  }
+*/
+// Make sense only for QuasiPointSource
+  virtual UnitVector GetVSrcCenter(); 
+  virtual double getAngSize() {
+  	cerr<<"call for getAngSize in a non Quasipt light source"<<endl;
+  	exit(-1);
+  	return 0;
+  }
+  
+// Optimisation for LightSource separable in frequency
+  virtual double spectre(double frequence) {
+  	frequence=0;
+        // Return the frequency spectrum independent of coordinates
+  	cerr<<"call for spectrum in an non separable light source"<<endl;
+  	exit(-1);
+  	return 0;
+  }
+  
+  virtual double powerDensAmpli(double theta,double phi) {
+  	theta=phi; //Pour eviter un warning
+  	// Return power density Amplidude at coordinates (W/m2/GHz/St) 
+  	cerr<<"call for powerDensAmpli in an non separable LightSource"<<endl;
+  	exit(-1);
+  	return 0;
+  }
+  double planckSpectrum(double freq, double Tempe=T_CMBR);
+  // W/m2/st/Hz
+  double planckSpecDiff(double DTempe, double freq, double Tempe=T_CMBR);
+  // W/m2/St/Hz
+protected:
+  bool PtSourceS;  // Is a set of PtSources (Galaxie IR ponctuelles?)
+  bool QPtSource;    // Is this a QuasiPointSource. (Sun, Jupiter, Mars et al.)
+  bool LScrFSep;	 // Are spatial coordinated and frequency separable variables?
+  bool IsLSrcMappedPower;	// Is the LSrc a in band map power. Valid only for Lobes separabe in freq.
+  double resolution; // Map Pixellisation resolution.;
+  char Name[32];
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/abslobenopolar.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/abslobenopolar.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/abslobenopolar.cc	(revision 799)
@@ -0,0 +1,23 @@
+// Dominique YVON, CEA/DAPNIA/SPP
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+//   #include "unixmac.h"
+#endif
+
+#include "abslobenopolar.h"
+
+int AbsLobeNoPolar::gorskyNlatFromRes() const {
+	double rlat= M_PI/2./lobeResol();
+    int nlat=2;
+    while(nlat<rlat) nlat*=2;
+    if (nlat>8192) cerr<<"warning: taille gorskyNlat: "<< nlat<<endl;
+    return nlat;
+}
+
+void AbsLobeNoPolar::print(ostream& OutStr) const
+{	OutStr<<"Lobe: "<<Name<<endl;
+	OutStr<<"Angle Max: "<< angleMax<<endl;
+	return;
+}
Index: /tags/SigPred_V1/SigPredictor/abslobenopolar.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/abslobenopolar.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/abslobenopolar.h	(revision 799)
@@ -0,0 +1,93 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef AbsLobeNoPolar_H
+#define AbsLobeNoPolar_H
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+//   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+
+#include <stdlib.h>
+#include "unitvector.h"
+	
+#include "lobespecresp.h"
+#include "radspec.h"
+
+class AbsLobeNoPolar{
+//friend class SigCalcTool;
+// Base class for Implementing Lobes 
+
+ public:
+  AbsLobeNoPolar(double FMin,double FMax):FreqMin(FMin),FreqMax(FMax) 
+  { cosanglemax=-1.;
+  	sprintf(Name,"");
+  }
+  virtual double weigth(const UnitVector& VInteg, const UnitVector& VP, 
+  	const UnitVector& VY, double freq) const =0;
+    // Return relative weigth of incoming power from direction VInteg
+    // at frequency freq no dimensions
+  virtual double lobeResol() const =0;
+    // Return reasonable computing resolution for this lobe
+  virtual double ResolutionCurve(double angle) const=0;
+ 	// Compute angular resolution degradation factor as a function of angle
+ 	// Usefull to optimize the speed of computation at large angle
+    // Overloaded when optimisation is required.
+  virtual void print(ostream& OutStr) const;
+  inline LobeSpecResp WeigthVsFreq(const UnitVector& VInteg, 
+  	const UnitVector& VP, const UnitVector& VY) 
+  {
+  	LobeSpecResp TheSpecResp(this,VInteg,VP,VY,this->minFreq(),this->maxFreq());
+  	return TheSpecResp;
+  }
+	// Return an object including WeigthVsFreq  with no dimensions.
+	
+/*  virtual DPixelMap& getLobeMap(double ResLSrcMap)=0;  
+ 	// Return Map of Lobe to allow multiplication of map 
+ 	// To be implemented
+*/  
+  inline double AngleMax() const{return angleMax;}
+  inline double minFreq() const {return FreqMin;}
+  inline double maxFreq() const {return FreqMax;}
+  inline bool IsFreqSep() const {return LobeFSep;}
+
+  int gorskyNlatFromRes() const;
+    // Compute Gorsky Nlat to prepare a map with full resolution of the lobe
+
+  virtual UnitVector VecShift(const UnitVector& VPointe,const UnitVector& VY) const{
+    // Usefull when peak sensitivity is not in boresight direction
+  	return VPointe;
+  }
+
+  //	Optimisations tools for frequency separable Lobe computations.
+  virtual double spectre(double freq) const{
+  	freq=0.;	//evite un warning
+  	if(!LobeFSep) cerr<< "appel a lobe.spectre() in a non separable lobe"<<endl;
+  	exit(-1);
+  	return 0;
+  }
+  virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  	const UnitVector& VY) const{
+  	if(!LobeFSep) cerr<< "appel a lobe.weigthAmpl() in a non separable lobe"<<endl;
+  	exit(-1);
+  	return 0;
+  }
+
+ protected:
+  char Name[32];
+  double FreqMin;
+  double FreqMax;
+  double angleMax;
+  double cosanglemax; 
+    // Angle beyond which constributions is neglected for this lobe
+    // Lobe 4 M_PI angleMax=180.
+  bool LobeFSep;	
+    // Are spatial coordinated and frequency separable variables?
+  inline double minToRad(double minute) { return minute/(180.*60.)*M_PI;}
+  inline double radToMin(double rad) { return rad/M_PI*180.*60.; }
+  //SpectralResponse* pSpecRep;
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/allfilter.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/allfilter.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/allfilter.h	(revision 799)
@@ -0,0 +1,4 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "squarefilt.h"
+#include "tabulatedfilter.h"
Index: /tags/SigPred_V1/SigPredictor/alllightsources.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/alllightsources.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/alllightsources.h	(revision 799)
@@ -0,0 +1,10 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "lightcmbprim.h"
+#include "lightdipole.h"
+#include "lightdiffdust.h"
+//#include "lightgalaxnoresol.h"
+//#include "lightgalaxresol.h"
+#include "lightsynchro.h"
+#include "lightsrcmappowerinband.h"
+#include "quasiptsources.h"
Index: /tags/SigPred_V1/SigPredictor/alllobe.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/alllobe.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/alllobe.h	(revision 799)
@@ -0,0 +1,4 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+#include "fastlobes.h"
+#include "lobes4pi.h"
+#include "lobecartographmoyen.h"
Index: /tags/SigPred_V1/SigPredictor/cartelobe.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/cartelobe.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/cartelobe.cc	(revision 799)
@@ -0,0 +1,200 @@
+// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stddef.h>
+#include <string.h>
+#include <math.h>
+#include <iostream>
+//#include "machine.h"
+#include "nbconst.h"
+#include "cimage.h"
+#include "cartelobe.h"
+
+
+ static NumRecipes NR;  // pour avoir acces au librairies NR
+
+//  Tout cet Objet est programme dans le philosophie des tableaux NR
+//  J'essaye d'etre coherent
+
+// Last Index used Value member
+long lastXIndex, lastYIndex;
+
+CarteLobe::CarteLobe(char fileRoot[], double frequence, int_4 NbPixelX, int_4 NbPixelY){
+  freq=frequence;
+  NbPixX=(long)NbPixelX;
+  NbPixY=(long)NbPixelY;
+  XMin=XMax=YMin=YMax=1.;
+  lastXIndex=lastYIndex=1;
+  
+  pCarte=NR.matrix(1,NbPixX,1,NbPixY);
+      // On Reserve l'espace memoire pour stocker la carte. Librairie NR
+  pAxeX= NR.vector(1,NbPixX);
+  pAxeY= NR.vector(1,NbPixY);  
+      // La suite a-t-elle un sens avec les fonction de la NR?
+  if ((pCarte==NULL)||(pAxeX==NULL)||(pAxeY==NULL)) {
+    cerr<<"erreur memoire dans constructeur pCarte"<<endl;
+    cerr<<"espace demand: "<<NbPixX*NbPixY*4<<" Octets"<<endl;
+  }
+
+  FILE* pfile;
+  char fileCur[150];
+
+  // On lit le tableau des abscisses
+  sprintf(fileCur,"%s",fileRoot);
+  strcat(fileCur,"_U.txt");
+  pfile=fopen(fileCur,"r");
+  if (pfile==NULL){
+    cerr<< "Erreur a l'ouverture du fichier de donnees :"<<fileCur<< endl;
+    exit(-1);
+  }
+
+  double PixCur=0.;
+  char dumChar;
+  for (int_4 i=1; i<(NbPixX+1); i++) {
+    fscanf(pfile,"%le%c",&PixCur,&dumChar);
+    pAxeX[i]=(r_4)PixCur;
+  }
+
+  XMin=XMax=pAxeX[1];
+  for (int_4 i=2; i<(NbPixX+1); i++){
+    if (XMin>pAxeX[i]) XMin=pAxeX[i];
+    if (XMax<pAxeX[i]) XMax=pAxeX[i];
+  }
+  fclose(pfile);
+ 
+  // On lit le tableau des ordonnees
+  sprintf(fileCur,"%s",fileRoot);
+  strcat(fileCur,"_V.txt");
+
+  pfile=fopen(fileCur,"r");
+  if (pfile==NULL){
+    cerr<< "Erreur a l'ouverture du fichier de donnees :"<< fileCur<< endl;
+    exit(-1);
+  }
+
+  for (int_4 j=1; j<(NbPixY+1); j++) {
+    fscanf(pfile,"%le%c",&PixCur,&dumChar);
+    pAxeY[j]=(r_4)PixCur;
+    for (int_4 i=0; i<(NbPixX-1); i++) {
+      fscanf(pfile,"%le%c",&PixCur,&dumChar);
+    }
+  }
+
+  YMin=YMax=pAxeY[1];
+  for (int_4 i=0; i<NbPixY; i++){
+    if (YMin>pAxeY[i]) YMin=pAxeY[i];
+    if (YMax<pAxeY[i]) YMax=pAxeY[i];
+  }
+  fclose (pfile);
+
+  // On lit le contenu de la carte 
+  sprintf(fileCur,"%s",fileRoot);
+  strcat(fileCur,"_PowdB.txt");
+  
+  pfile=fopen(fileCur,"r");
+  if (pfile==NULL){
+    cerr<< "Erreur a l'ouverture du fichier de donnees :"<< fileCur<< endl;
+   exit(-1);
+  }
+
+  for (int j=1; j<(NbPixY+1); j++){
+    for(int i=1; i<(NbPixX+1); i++){
+      fscanf(pfile,"%le%c",&PixCur,&dumChar);
+      pCarte[i][j]=PixCur;
+    }
+  }
+  resolution=CalcResol();
+  // C'est fini!
+}
+
+CarteLobe::~CarteLobe() {
+  NR.free_vector(pAxeX,1,NbPixX);
+  NR.free_vector(pAxeY,1,NbPixY);
+  NR.free_matrix(pCarte,1,NbPixX,1,NbPixY);
+}
+
+double CarteLobe::Value(double XVal, double YVal) const {
+	// Value of lobe at coordinates
+  float X=(float) XVal;
+  float Y=(float) YVal;
+  float dBValue=0.;
+  float dBError=0.;
+
+  if ((X<=XMin)||(X>=XMax)||(Y<=YMin)||(Y>=YMax)) return 0.;
+ 
+  else {
+    unsigned long XIndex=lastXIndex;
+    unsigned long YIndex=lastYIndex;
+    
+    NR.hunt(pAxeX,NbPixX,X,&XIndex);  //  0<X<NbPixX-1
+    NR.hunt(pAxeY,NbPixY,Y,&YIndex);  //  0<Y<NbPixY-1
+
+    lastXIndex=XIndex;
+    lastYIndex=YIndex;
+
+    // Extrait le soustableau encadrant XIndex et YIndex
+    if(XIndex>(NbPixX-2)) XIndex=NbPixX-2;
+    if(YIndex>(NbPixY-2)) YIndex=NbPixY-2;
+    float ** ptableau= NR.submatrix(pCarte,XIndex,XIndex+2,YIndex,YIndex+2,1,1);
+    
+    // On interpole sous le soustableau
+    NR.polin2(pAxeX+XIndex-1,pAxeY+YIndex-1,ptableau,3,3,X,Y,&dBValue,&dBError);
+
+    // On passe en echelle lineaire 
+    return pow(10.,(double)(dBValue/10.));
+  }
+}
+
+void CarteLobe::FitsVisu(char FitsFile[],long nx,long ny) {
+  Image<float>* img;
+  img = new Image<float>(nx,ny,1);
+  double xstep=(XMax-XMin)/nx;
+  double ystep=(YMax-YMin)/ny;
+  double yPos=YMin;
+  double xPos=XMin;
+
+  for(int i=0; i<nx;i++){
+    xPos=XMin;
+    for(int j=0; j<ny; j++){
+     (*img)(i,j)=Value(xPos,yPos);
+      xPos+=xstep;
+    }
+    yPos+=ystep;
+  }
+  
+  (*img).Save(FitsFile);    // BUGG XXXXXX 
+  							// Ne sauvegarde pas dans un format FITS pour le moment!!!
+  cerr<<"Ecriture fichier :"<<FitsFile<<endl;
+  delete img;
+}
+
+void CarteLobe::UVToAng(double U, double V, double& alpha, double& beta) const {
+  // U=sin alpha cos beta
+  // V=sin alpha sin beta
+  double sinalpha=sqrt(U*U+V*V);  	// Approx des petits angles pour alpha
+                               		// Pas pour beta
+  if(U>=0) beta=asin(V/sinalpha);
+  else beta=Pi-asin(V/sinalpha);
+  alpha=asin(sinalpha);
+}
+
+double CarteLobe::halfAngAperture(){
+  double alCentre=0.; double beCentre=0.; 
+  UVToAng(UCentre(),VCentre(),alCentre,beCentre);
+  double alBord=0; double beBord=0;
+  UVToAng(XMin, YMin, alBord,beBord);
+  double CosAng=0.;
+  CosAng+=cos(alCentre)*cos(alBord);
+  CosAng+=sin(alCentre)*sin(beCentre)*sin(alBord)*sin(beBord);
+  CosAng+=sin(alCentre)*cos(beCentre)*sin(alBord)*cos(beBord);
+  return acos(CosAng);
+}
+
+double CarteLobe::CalcResol(){
+  double halfAng=halfAngAperture();
+  double NbBin=sqrt(NbPixX*NbPixX+NbPixY*NbPixY);
+  return 2*halfAng/NbBin;		// C'est delirament petit
+//  on prend
+//	return 1.e-3;  // Debug XXXX
+}
Index: /tags/SigPred_V1/SigPredictor/cartelobe.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/cartelobe.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/cartelobe.h	(revision 799)
@@ -0,0 +1,62 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef CARTELOBE_SEEN
+#define CARTELOBE_SEEN
+
+#include "abslobenopolar.h"
+
+#define NRANSI
+#include "numrecipes.h"
+/* Les lobes cartographies ne presentant aucune symetrie spatialle, il est necessaire 
+de connaitre l orientation du repere Ticra pour pouvoir calculer dans les lobes. 
+Ceci est incompatible avec l'etablissement a priori de carte du ciel mesures. 
+Il faut connaitre la strategie de scanning et alors reconstruire les cartes, bref 
+faire tout le boulot.
+*/
+
+class CarteLobe {
+    /* les cartes de lobe centraux sont stockees selon les variables 
+       U=sin(alpha)cos(beta)
+       V=sin(alpha)sin(beta)  ou alpha et beta sont les angle en coordonnees 
+       polaires dans le repere du logiciel TICRA
+       Ce sont ces variables qu'il faut passer a la fonction Value pour obtenir
+       une reponse coherente
+    */
+    // On utilise les Numerical Recepies:gestions de tableau particulieres
+ public:
+    CarteLobe(char fileRoot[]="cartelobevide",double frequence=100.E9,int_4 NbPixelX=100,
+    int_4 NbPixelY=100);
+    virtual ~CarteLobe();
+    void FitsVisu(char FitsFile[],long NbX=300,long NbY=300);
+    double Value(double UValue,double VValue) const ;
+    // Retourne la valeur de la carte interpole en UValue (X) et en VValue(Y)
+    // Echelle lineaire.
+    double getFreq() const  { return freq;}
+    void UVToAng(double U, double V, double& alpha, double& beta) const;
+       //fonction ecrite dans la limite des petits angles
+    double UCentre(){return (XMin+XMax)/2.;}
+    double VCentre(){return (YMin+YMax)/2.;}
+    void UBornes(double& UMin, double& UMax) const { 
+	UMin=(double)XMin; UMax=(double)XMax; return;}
+    void VBornes(double& VMin, double& VMax) const {
+	VMin=(double)YMin; VMax=(double)YMax; return;}
+    double halfAngAperture();
+	inline double resol() {return resolution;}
+ protected:
+    //NumRecipes NR;  // pour avoir acces au librairies NR
+    long NbPixX;
+    long NbPixY;
+    float *pAxeX;
+    float *pAxeY;
+
+    float XMin, XMax;
+    float YMin, YMax;
+    float** pCarte;  // Pointeur sur la carte du lobe,facon NR
+                     // Les valeurs sont stockees en decibel
+    double freq;     // Frequence a laquelle la carte du lobe a ete etablie
+    double resolution; //Resolution angulaire de la carte de base
+
+    double CalcResol();
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/distribflatfilter.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/distribflatfilter.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/distribflatfilter.cc	(revision 799)
@@ -0,0 +1,32 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+
+#include "distribflatfilter.h"
+
+//-------------CC Code--------------------------------------------
+
+//----------------
+// Constructor --
+//----------------
+DistribFlatFilter::DistribFlatFilter(double numin, double numax)
+  : SpectralResponse(numin, numax), UnSurdeltaFreq(1./(numax-numin))
+{	multFactor=1.;
+}
+
+//--------------
+// Destructor --
+//--------------
+DistribFlatFilter::~DistribFlatFilter()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+double DistribFlatFilter::transmission(double nu) const
+{
+  if(nu>_numin && nu<_numax) return UnSurdeltaFreq*multFactor;
+  return 0.;
+}
+
+
Index: /tags/SigPred_V1/SigPredictor/distribflatfilter.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/distribflatfilter.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/distribflatfilter.h	(revision 799)
@@ -0,0 +1,50 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+
+#ifndef DistribFlatFilter_SEEN
+#define DistribFlatFilter_SEEN
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+
+//----------------------
+// Base Class Headers --
+//----------------------
+#include "specresp.h"
+
+
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+class DistribFlatFilter:public SpectralResponse
+{
+public:
+  //Constructors
+  DistribFlatFilter (double numin, double numax); 
+  
+  // destructor
+  virtual ~DistribFlatFilter();
+
+  //
+  //  Function members
+  //
+
+  double transmission(double) const;
+  void SetMultFactor(double mult) {
+  	multFactor=mult;
+  }
+protected:
+  double UnSurdeltaFreq;
+  double multFactor;
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/fastlobes.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/fastlobes.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/fastlobes.h	(revision 799)
@@ -0,0 +1,81 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef FastLOBE_SEEN
+#define FastLOBE_SEEN
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+#include "meanlobe.h"
+
+
+class LobeGaussien:public MeanFreqLobe 
+{ 
+ public:
+    LobeGaussien(double FWHM_degre, double freqMin, double freqMax)
+      :MeanFreqLobe(freqMin,freqMax){
+		sigma=minToRad(FWHM_degre/2.354*60.);
+		angleMax=3.*sigma;
+//		angleMax=1.5*sigma; // BUGG XXXX pour debug
+		cosanglemax=cos(angleMax);
+		sprintf(Name,"Lobe Gaussien");
+    }
+    virtual ~LobeGaussien(){ };
+   	virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  	  const UnitVector& VY) const {
+   	  	double cosinus=VP*VInteg;		
+	    if(cosinus>cosanglemax) 
+	    {  if(cosinus<1.) 
+	    	{ 	double ang= acos(cosinus);
+	    		//return tabFExp(-ang*ang/(2*sigma*sigma));
+	       		return exp(-ang*ang/(2*sigma*sigma));
+	       	}
+	       	else return 1.;
+	    }
+		else return 0.;
+	}
+	virtual double lobeResol() const{
+		return sigma*2.354/3.; // On veut pixeliser/calculer  au moins FWHM/3.
+	} 
+	virtual double ResolutionCurve(double angleShift) const{
+		return 1.;
+		//return (double) (1.+(angleShift/sigma/7.));
+	}
+ protected: 
+    double sigma;
+};
+
+class LobeConique:public MeanFreqLobe {
+    // Usefull to degrade map resolutions without loosing points objects
+ public:
+    LobeConique(double halfAngle_degre, double freqMin, double freqMax)    
+    	:MeanFreqLobe(freqMin,freqMax){
+		halfAng=minToRad(halfAngle_degre*60.);
+	//	cosHalfAng=cos(halfAng);
+		angleMax=halfAng;
+		cosanglemax=cos(halfAng);
+		sprintf(Name,"Lobe Conique");
+    }
+    virtual ~LobeConique() {   }
+
+	virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  		const UnitVector& VY) const {
+		double cosinus=VP*VInteg;
+		if(cosinus>cosanglemax) return 1.;
+		else return 0.;
+	}
+   	virtual double lobeResol() const{
+		return halfAng/1.5;
+   	}
+   
+   virtual double ResolutionCurve(double angle) const {
+    	return 1.;
+    }
+ protected:
+    double  halfAng;    		// demi angle d ouverture du cone
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightcmbprim.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightcmbprim.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightcmbprim.cc	(revision 799)
@@ -0,0 +1,51 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "lightcmbprim.h"
+#include "fitsioserver.h"
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "macenvvariables.h"
+#endif
+
+
+LightCMBPrim::LightCMBPrim(int_4 nside)
+	:nlat(nside) {
+
+  resolution=2*3.14159/(nlat*4.);  			// Radian
+  sprintf(Name, "CMB Primordial");
+  cout<<"Objet LightCMBPrim :"<<endl;
+  cout<<"Resolution de travail nlat= "<<nlat<<endl;
+  
+  // On reserve la mmoire pour la carte gorski
+  try {
+  	pMap =new SphereGorski<float> (nlat);
+  }
+  catch (bad_alloc) {
+  	cerr<<" Memory booking error in LightCMBPrim"<<endl;
+  	exit(-1);
+  }
+  cout<<"Vous avez rserv: "<<12*nlat*nlat*sizeof(float)<<" Octets"<<endl;
+
+  // On lit le fichier gnr par synfast
+ 			
+  char filename[150];  
+
+#ifndef __MWERKS__
+  char* PATHDataLScr=getenv("PATHDataLScr"); 
+  sprintf(filename,"%s/cmb_res%04i.fits",PATHDataLScr,nlat);
+#else
+  sprintf(filename,"%scmb_res%04i.fits",PATHDataLScr,nlat); 
+#endif
+
+  FitsIoServer FitsServer;
+  FitsServer.load(*pMap,filename);	// Units:  Microkelvin
+}
+
+LightCMBPrim::~LightCMBPrim() {
+	delete pMap;
+}
+
Index: /tags/SigPred_V1/SigPredictor/lightcmbprim.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightcmbprim.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightcmbprim.h	(revision 799)
@@ -0,0 +1,28 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LightCMBPrim_SEEN
+#define LightCMBPrim_SEEN
+#include "spheregorski.h" 
+#include "lightsrcsep.h"
+
+class LightCMBPrim:public LightSrcDiffBlackBody
+{
+ public:
+    LightCMBPrim(int_4 nside=256);
+    	// nside de la carte Gorsky
+    	// Variane required for sky fluctuationss
+    virtual ~LightCMBPrim();
+ protected:
+ 	double Variance;
+    SphereGorski<float>*  pMap;
+    int_4 nlat;
+    
+    virtual double DeltaT(double theta,double phi) {
+    	int_4 index=pMap->PixIndexSph(theta,phi); 	
+    		// La carte est en DT/T * 1.e6 
+  		double valeur=pMap->PixVal(index)*T_CMBR*1.e-6;
+  		return valeur; // En Kelvin
+    }
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightdiffdust.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightdiffdust.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightdiffdust.cc	(revision 799)
@@ -0,0 +1,66 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+//#include "integ.h"
+#include "lightdiffdust.h"
+#include "fitsioserver.h"
+
+LightDiffDust::LightDiffDust(int_4 nside=256) {
+  nlat=nside;
+  resolution=2*3.14159/4./nlat; // Resolution angulaire de la carte lue pour ce fond.
+                      			// Radians
+  sprintf(Name,"Poussire Diffuse Schlegel");                  			
+   
+  // Modle Physique                  			
+  TDust=18.0;
+  index=2.;
+  double lref=240.e-6;		// Longueur d'onde carte schlegel
+  fRef=3.e8/lref;			// Frequence de la carte Schelgel
+							// Car nous parlons en frequence
+  index2=-0.15;	
+  Hfraction=0.05;
+  double kbol=1.38e-23;
+  double freeAtPole=5.3e-6;
+  double freeAtPoleMJy_Sr=freeAtPole*2.*kbol
+  	* pow(53.e9/3.e8,2) * pow(fRef/53.e9,index2) / 1.e-20;
+  double dustAtPoleMJy_Sr=3*0.9/1.1;
+  CorrelatedfreeWeight=freeAtPoleMJy_Sr/dustAtPoleMJy_Sr;
+//  nuPlanckfRef=nuPlanck(fRef);
+  
+  cout<<"Objet LightDiffDust :"<<endl;
+  cout <<"Resolution de travail nlat= "<<nlat<<endl;
+  
+  try { pMap=new SphereGorski<float> (nlat);	}
+  catch (bad_alloc) {
+    cerr<<"memory booking error Lightdiffdust constructor nlat= "<<nlat<<endl;
+    exit(-1); 
+  }
+  cerr<<"Vous avez rserv: "<<12*nlat*nlat*sizeof(float)<<" Octets"<<endl;
+  char filename[150];
+  
+#ifndef __MWERKS__  
+  char *PATHDataLScr;
+  PATHDataLScr=getenv("PATHDataLScr");
+#endif
+
+  sprintf(filename,"%sdust_res%04i.fits",PATHDataLScr,nlat);
+	
+  cout<<"On lit la carte "<<filename<<endl;
+  if(nlat<1024) cerr<<"la carte avec nlat=1024 est dispo"<<endl;
+  
+  FitsIoServer FitsServer;
+  FitsServer.load((*pMap),filename);
+		// unite: MJy/st. Dixit Romain
+}
+
+LightDiffDust::~LightDiffDust(){
+  delete pMap;
+}
Index: /tags/SigPred_V1/SigPredictor/lightdiffdust.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightdiffdust.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightdiffdust.h	(revision 799)
@@ -0,0 +1,44 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LightDIFFDUST_SEEN
+#define LightDIFFDUST_SEEN
+#include "spheregorski.h" 
+#include "lightsrcsep.h"
+
+
+class LightDiffDust:public LightSrcSep {
+ public:
+    LightDiffDust(int_4 nside);
+    virtual ~LightDiffDust();
+ protected: 
+    int nlat;           // Nlat de la carte gorsky de reference
+    double TDust;  		// Temperature of dust
+    double fRef;   		// Frequency Reference
+    double index;  		// Index en frequence fitt sur le modle
+    double Hfraction;
+    double index2;
+   	double CorrelatedfreeWeight;
+    double nuPlanckfRef;
+    SphereGorski<float>* pMap;         // Map units MJy from IRAS Data.
+    
+    virtual double powerDensAmpli(double theta,double phi){
+    	return (*pMap)(theta,phi)*1.e-20*fRef;
+  		// Return power density Amplidude at coordinates
+    	// (W/m2/Hz/st)
+    }
+    virtual double spectre(double freq){
+    /*	double temp=planckSpectrum(freq,TDust)/planckSpectrum(fRef,TDust)
+    			*(1-Hfraction)* freq*freq/fRef/fRef;
+    	return temp;
+    */
+  	double nuPnu= nuPlanck(freq)/nuPlanck(fRef) * (1.-Hfraction) * pow(freq/fRef,index) 
+    	+ pow(freq/fRef,1+index2) * CorrelatedfreeWeight;	
+    	return nuPnu/freq;	
+    }
+    inline double nuPlanck(double freq){
+    	return freq*planckSpectrum(freq,TDust);
+    }
+};
+
+#endif
+ 
Index: /tags/SigPred_V1/SigPredictor/lightdipole.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightdipole.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightdipole.cc	(revision 799)
@@ -0,0 +1,74 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "math.h"
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+
+#include "lightdipole.h"
+
+
+LightDipole::LightDipole(){
+  resolution=0.; // process computed analytically. 
+  betaDipole=1.231e-3; 
+  phiDipole= 264.31/180.*M_PI;
+  thetaDipole= 48.05/180.*M_PI; 
+  sprintf(Name, "Dipole du CMB");
+  
+  cout<<" Creation objet LightDipole" <<endl;
+
+  vDipole[0]= cos(phiDipole)*sin(thetaDipole);
+  vDipole[1]= sin(phiDipole)*sin(thetaDipole);
+  vDipole[2]= cos(thetaDipole);
+}
+
+
+double LightDipole::DeltaT(double theta,double phi)
+{
+//Calcule /DT/TCMB apparent du corps noir dans la direction theta, phi
+//Peebles and Wilkinson  Phys Rev 174 (1968) 2168
+  double vUnit[3];   // Vecteur unitaire correspondant a theta et phi
+  double cosinus;    // cosinus de l'angle entre le pointe et VDipole
+
+  vUnit[0]= cos((double)phi)*sin((double)theta);
+  vUnit[1]= sin((double)phi)*sin((double)theta);
+  vUnit[2]= cos((double)theta);
+
+  cosinus=0.;      //cosinus de l'angle entre le pointe et VDipole
+  for (int j=0; j<3; j++) cosinus=cosinus+vUnit[j]*vDipole[j];
+
+  return betaDipole*T_CMBR*cosinus;      // DT/TCMB
+}
+
+double LightQuadrupole::DeltaT(double theta,double phi){
+  //Calcule /DT/TCMB apparent du corps noir dans la direction theta, phi
+  //Peebles and Wilkinson  Phys Rev 174 (1968) 2168
+  double vUnit[3];   // Vecteur unitaire correspondant a theta et phi
+  double cosinus=0.;  // cosinus de l'angle entre le pointe et VDipole
+
+  vUnit[0]= cos((double)phi)*sin((double)theta);
+  vUnit[1]= sin((double)phi)*sin((double)theta);
+  vUnit[2]= cos((double)theta);
+
+   //cosinus de l'angle entre le pointe et VDipole
+  for (int j=0; j<3; j++) {cosinus+=vUnit[j]*vDipole[j];}
+
+  return betaDipole*betaDipole*T_CMBR*(cosinus*cosinus-0.5);  
+                  //DT/TCMB
+}
+
+LightBlackBody::LightBlackBody(double TempeCN, double Resol)
+	:TempeCNoir(TempeCN) 
+{	resolution=Resol;
+	sprintf(Name,"CorpsNoirTempe:%5.1f K", TempeCN);
+}
+	
+LightNormTCMB::LightNormTCMB(double DTempe, double Res)
+	:DeltaTempeCNoir(DTempe) 
+{	resolution=Res;
+	sprintf(Name, "Normalisation TCMB");
+}
Index: /tags/SigPred_V1/SigPredictor/lightdipole.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightdipole.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightdipole.h	(revision 799)
@@ -0,0 +1,63 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SKYDIPOLE_SEEN
+#define SKYDIPOLE_SEEN
+
+#include "lightsrcsep.h"
+
+class LightDipole:public LightSrcDiffBlackBody {
+ public:
+    LightDipole();
+    virtual ~LightDipole() {   };
+ protected:
+    virtual double DeltaT(double theta,double phi);
+// Additional members:
+// Cosmic Microwave Background Parameters as measured by COBE and COBRA
+    double betaDipole;
+    double phiDipole;
+    double thetaDipole;
+    double vDipole[3];
+};
+
+class LightQuadrupole:public LightDipole {
+    //strictement identique a Lightdipole sauf processValue()
+    //ou on calcule D(T)/T  quadrupolaire
+ public:
+    LightQuadrupole():LightDipole(){   };
+    virtual ~LightQuadrupole() { };
+ protected:
+     virtual double DeltaT(double theta,double phi);
+};
+
+class LightBlackBody:public LightSrcSep {
+ public:
+    LightBlackBody(double Tempe=T_CMBR, double Resolution=0.);
+    virtual ~LightBlackBody() {   };
+    
+    virtual double powerDensAmpli(double theta,double phi)
+    {	return 1.; }
+    
+    virtual double spectre(double frequence)
+    { return planckSpectrum(frequence,TempeCNoir);
+    }
+ protected:
+	double TempeCNoir;
+};
+
+
+class LightNormTCMB:public LightSrcSep {
+ public:
+    LightNormTCMB(double Tempe=1., double Resol=0.);
+    virtual ~LightNormTCMB() {   };
+    
+    virtual double powerDensAmpli(double theta,double phi)
+    {	return 1.; } 	// En Kelvin
+    
+    virtual double spectre(double frequence)
+    {	 return AbsLightSource::planckSpecDiff(DeltaTempeCNoir,frequence,T_CMBR);
+    }
+    
+ protected:
+	double DeltaTempeCNoir;
+};
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.cc	(revision 799)
@@ -0,0 +1,165 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+#include "nbrandom.h"
+#include "tmatrix.h"
+#include "numrecipes.h"
+#include "integ.h"
+#include "fitsioserver.h" 
+
+#include "lightgalaxnoresol.h"
+
+
+LightGalaxNoResol::LightGalaxNoResol(int_4 nside=512)
+{
+ 
+ /* cerr<< "LightGalaxNoResol en travaux, ne pas utiliser"<<endl;
+    exit (-1);
+ */
+  sprintf(Name,"Galaxies non resolues");
+  nlat=nside;              
+  if(nside!=2048) 
+  {  cerr<<" Warning : Donnees Guiderdonni compatible avec nlat=2048"<<endl;
+  	 cerr<<" Vous avez choisi nlat= "<<nlat<<endl;
+  }
+  nbFreq=31;			   // Nbre de frequences calcules 
+  NbPixUnCote=441;
+  nbPixel=NbPixUnCote*NbPixUnCote;    // Nbre de pixel cartes Guiderdonni
+  nbPixelKept=65535;       // Ce que lon peut stocker dans un U2
+  pNR= new NumRecipes;
+
+  // Tableau des frquences auxquelles sont calcules les cartes guiderdonni
+  dataFreq=new r_8[nbFreq]; 
+  dataFreqDegueux= new float[nbFreq];
+  // Calcul des frquences des cartes guiderdonni
+  double fmin=15.e9;
+  double fmax=15.e12;
+  double dLogFreq=log10(fmax/fmin)/(nbFreq-1.);
+  
+  for(int i=0; i<nbFreq; i++) 
+  {
+    dataFreq[i]=fmin*pow(10,i*dLogFreq);
+    dataFreqDegueux[i]=(float)dataFreq[i];
+  }
+  // Chaque pixel du ciel est caracterise par un spectre chantillonn
+  //  nbfreq frequences extrait des cartes Guiderdonni.
+  // On stocke le no de pixel de la carte guiderdonni correspondant au pixel ciel
+  
+  try {pSphereIndex=new SphereGorski<uint_2>(nlat); }
+  catch(bad_alloc) {
+  	cerr<<"Memory booking error in LightGalaxNoResol: pSphereIndex"<<endl;
+  	exit(-1);
+  }
+    // Sphre faisant la correspondance ciel, no de pixel de carte guiderdonni
+  nbPixelLight=pSphereIndex->NbPixels();
+  	// On attribue a chaque pixel ciel un pixel carte Guiderdonni.
+  for( long i=0; i< nbPixelLight; i++) 
+  { 
+    // double test= MyRan()*nbPixelKept;
+    // bogue sur Mac frand01() retourne toujours la meme valeur.
+    // On utilise ran() standard.
+    
+    (*pSphereIndex)(i)= (int_2) (MyRan()*nbPixelKept);
+  }
+  
+    // Stockage des donnes Guiderdonni petites cartes 
+        // Rserve l'espace mmoire ncessaire au stockage des cartes partielles
+  ppskyData= new r_4* [nbPixelKept];	
+  		// Un tableau pour chaque point du ciel
+  try {
+  	for(int i=0; i<nbPixelKept; i++) ppskyData[i]= new r_4 [nbFreq];
+  }
+  catch(bad_alloc) {
+  	cerr<<"Memory booking error in LightGalaxNoResol: ppskyData"<<endl;
+  	exit(-1);
+  }
+
+  
+/*// ca ne marche plus. se renseigner de comment faire!.
+  if(pSphereIndex==NULL){
+    cerr<<"Memory booking error LightGalaxNoResol constructor"<<endl;
+    exit(-1);
+  }  
+*/
+
+  int_4 size=sizeof(r_4)*(nbFreq*nbPixelKept) + sizeof(uint_2)*nbPixelLight;
+  cout<<"Objet LightGalaxNoResol: "<<endl;
+  cout<<"Vous avez reserv: "<<size<<" Octets de mmoire vive"<<endl;
+
+  // Lit les cartes Guiderdonni:les charge dans dataTampon,puis dans skyData
+  #ifndef __MWERKS__
+    char* PATHDataLScr=getenv("PATHDataLScr");
+  #endif
+  
+  FitsIoServer FitsServer;
+  
+  TMatrix<float> UneCarte(NbPixUnCote, NbPixUnCote);
+  
+  long cartefreqMHz; 
+  char filename[150];
+  for(int noFreq=0; noFreq<nbFreq; noFreq++)
+  {
+  	  // Gnre le nom des cartes
+	cartefreqMHz=(long)(dataFreq[noFreq]/1000000.+0.5);
+    sprintf(filename,"%seran102_441_%08imhz.fits",PATHDataLScr,cartefreqMHz);
+
+	FitsServer.load(UneCarte,filename);
+//	cout<<"Carte Lue: "<<filename<<'\t'<<"Frquence: "<<dataFreq[noFreq]<<endl;
+	
+	long indexUneCarte1=0;
+	long indexUneCarte2=0;
+	
+    // On remplit LightData 
+    for(int NoPixel=0;NoPixel<nbPixelKept;NoPixel++)
+    { 
+      indexUneCarte1= (long) (NoPixel/NbPixUnCote); 		// Division entiere
+      indexUneCarte2= NoPixel%NbPixUnCote;
+      ppskyData[NoPixel][noFreq]= UneCarte(indexUneCarte1,indexUneCarte2);
+      cout <<"NoPixel: "<<NoPixel<<" Valeur: ";
+      cout << ppskyData[NoPixel][noFreq]<<'\t'<<UneCarte(indexUneCarte1,indexUneCarte2)<<endl;
+    }
+    
+  }  // C'est dans la boite LightData
+  
+  // C'est fini pour le constructeur
+}
+
+LightGalaxNoResol::~LightGalaxNoResol(){
+ for(int i=0; i<nbPixelKept; i++) delete[] ppskyData[i];
+ delete[] ppskyData;
+ delete pNR;
+ delete[] dataFreq;
+ delete[] dataFreqDegueux;
+ delete[] pSphereIndex; 
+
+}
+
+double LightGalaxNoResol::powSpecDens(double theta, double phi, double freq){
+  long index=pSphereIndex->PixIndexSph(theta,phi);
+  int_2 pixNb=pSphereIndex->PixVal(index);
+
+  float* pbinSpectre= ppskyData[pixNb]; //10-26 W/m2/Hz/st
+ 
+  // On interpole au troisieme ordre.
+  float InterpRes=0.;
+  float InterpResErr=0.;
+  pNR->polint(dataFreqDegueux-1,pbinSpectre-1,3,freq,&InterpRes,&InterpResErr);
+
+  return (double)InterpRes*1.e-26;  // W/m2/Hz/st
+}
+
+
+static double randMax=RAND_MAX;
+
+double LightGalaxNoResol::MyRan() 
+{
+	return rand()/randMax;
+}
+
Index: /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightgalaxnoresol.h	(revision 799)
@@ -0,0 +1,40 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SKYGALAXNORESOL_SEEN
+#define SKYGALAXNORESOL_SEEN
+
+#include "spheregorski.h"
+#include "abslightsource.h"
+#include "numrecipes.h"
+
+
+class LightGalaxNoResol:public AbsLightSource {
+ public: 
+   LightGalaxNoResol(int_4 nside); 
+   virtual ~LightGalaxNoResol(); 
+   virtual double powSpecDens(double theta,double phi,double freq);
+    // Return power density at coordinates
+    // (W/m2/GHz/St)
+    
+ protected: 
+   int_4 nlat;          // Normalement c est 2048! Pas assez de mmoire. 
+   int_4 nbPixelLight;  // Carte du ciel complte
+   int_4 NbPixUnCote;	//Carte guiderdonni carreesq de NbPixUnCote sur un cote
+   int_4 nbFreq; 		  // Nbre de frquences chantillonnes
+   int_4 nbPixel;       // Nbre de pixel carte Guiderdonni
+   uint_2 nbPixelKept;  // Ce que on peut stocker dans un U2
+
+   NumRecipes* pNR;
+   r_8* dataFreq;     
+     // Tableau des frequences auxquelles sont calcules les cartes guiderdonni
+   float* dataFreqDegueux;
+     // Idem. Ne sert que pour appeller la fonction NR Polint
+   SphereGorski<uint_2>* pSphereIndex;
+     // Tableau de correspondance no de pixel carte sphereGorski, carte guiderdonni
+
+   r_4** ppskyData; // Espace mmoire ncessaire au stockage des cartes partielles
+	       		  // Unite Jansky.
+   double MyRan();
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightgalaxresol.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightgalaxresol.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightgalaxresol.cc	(revision 799)
@@ -0,0 +1,194 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <math.h>
+#include <iostream>
+#include <fstream>
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+
+#include "nbrandom.h"
+#include "lightgalaxresol.h"
+#include "numrecipes.h"
+#include "integ.h"
+
+
+LightGalaxResol::LightGalaxResol(int_4 nside) {
+  sprintf(Name, "Galaxies Resolues");
+  nlat=nside;
+  resolution=2*M_PI/4./nlat;		// Radians
+  nbFreq=31;
+  nbPointSource=2711;
+  PtSourceS=true;
+  
+  pNR=new NumRecipes;
+  LastFreqIndex=1;
+  
+// Calcul des frquences des cartes guiderdonni
+  dataFreq=new  r_8[nbFreq];
+  dataFreqDegueux= new float[nbFreq];
+
+  double fmin=15.e9;
+  double fmax=15.e12;
+  double dLogFreq=log10(fmax/fmin)/(nbFreq-1);
+
+  for(int i=0; i<nbFreq; i++) 
+  {
+    dataFreq[i]=fmin*pow(10,i*dLogFreq);
+    dataFreqDegueux[i]=(float)dataFreq[i];				// Pour la NR
+  }
+
+// Booking of skymap of index
+  try { pLightMap= new SphereGorski<uint_2> (nlat); }
+  catch(bad_alloc){
+    cerr<<"Erreur Reservation de memoire classe LightGalacResol, nlat= "<<nlat<<endl;
+    exit(-1);
+  }
+  
+  nbPixelLight=pLightMap->NbPixels();
+  int size=sizeof(uint_2)*nbPixelLight;
+  cout<<"Objet LightGalaxResol: "<<endl;
+  cout<<"Vous avez reserv: "<<size<<"Octets de mmoire vive"<<endl;
+  
+  
+// La sphere est initialise  nbPointSource+1, cad pas de sources
+  int_2 bidon=nbPointSource+1;
+  for(int i=0; i<nbPixelLight; i++)  (*pLightMap)(i)=bidon;
+  
+// On distribue les sources uniformement sur le ciel
+  int_4 skyIndex;
+  for(int_2 i=0; i<nbPointSource; i++) 
+  {
+    skyIndex=(int)(MyRan()*nbPixelLight);
+    if ( (*pLightMap)(skyIndex)==(nbPointSource+1)) (*pLightMap)(skyIndex)=i;
+    else 
+    { 
+      while(! ((*pLightMap)(skyIndex) == (nbPointSource+1)) ) 
+      { 
+//      float test;        
+//      test= MyRan();		// frand01() ne fonctionne pas au 05/12/99 sur Mac
+        skyIndex=(int)(MyRan()*nbPixelLight);
+      }
+  	   (*pLightMap)(skyIndex)=i;
+    } 
+  }
+
+
+// _____________ On lit les parametres des sources _______________________________________
+
+	// Reservation de l'espace memoire de stockage
+  try {	
+ 	ppPointSourceData= new float* [nbPointSource];
+  	for(int i=0; i<nbPointSource; i++) ppPointSourceData[i]= new float[nbFreq];
+  }
+  catch(bad_alloc) {
+  	cerr<<" Memory booking error in LightGalaxResol: ppskyData "<<endl;
+  	exit(-1);
+  }
+
+	// Lecture des fichiers
+  char filename[150];
+  double DumStarNb;
+  char DumChaine[20];
+  r_4 thisdata;
+
+
+#ifndef __MWERKS__
+  char* PATHDataLScr=getenv("PATHDataLScr");   
+#endif
+
+  for(int NoFreq=0; NoFreq<nbFreq; NoFreq++) 
+  {	 
+     // cout<<"NoFreq :"<<NoFreq<<"Frequence: "<<dataFreq[NoFreq]<<endl;
+     // On genere le bon nom de fichier
+  	 sprintf(filename,"%seran102_441_%08imhz.s_rare",PATHDataLScr,(int)(dataFreq[NoFreq]/1.e6+.5));
+ 
+  	 ifstream InFile(filename,ios::in);		// 	Fichier Ouvert!
+  	 
+#ifdef __MWERKS__
+ 	 if (!InFile.is_open())		// cxx DEC est trop con pour comprendre is_open()
+#else
+ 	if (!InFile)
+#endif
+  	 {	cerr<<"erreur a l'ouverture du fichier de donnees: "<<filename<<endl;
+  	 	exit(-1);
+  	 }
+
+
+  	 InFile>>DumStarNb; 
+  	 if(DumStarNb!=nbPointSource) {
+  		cerr<<"ya une bogue dans le nb de galaxptsources."<<endl;
+  		cerr<<"fichier: "<< filename<<"nb sources annoncees:"<< DumStarNb<<endl;
+  	 }
+  	 InFile>>DumChaine;
+  	 
+  	 for(long noPtSource=0; noPtSource<nbPointSource; noPtSource++) 
+  	 {
+  		InFile>> thisdata;
+  		ppPointSourceData [noPtSource] [NoFreq] =thisdata;
+//  		cout<<noPtSource<<'\t'<<ppPointSourceData [noPtSource] [NoFreq]<<endl;
+  	 }
+  	 InFile.close();
+  }
+/*     for (int NoFreq=0; NoFreq<nbFreq; NoFreq ++) 
+  	   cout<< NoFreq <<'\t'<< ppPointSourceData [] [NoFreq]<< endl;
+*/
+}
+
+
+LightGalaxResol::~LightGalaxResol() 
+{
+  // Librons la mmoire!! Bientt mai!
+  for(int i=0; i<nbPointSource; i++) delete[] ppPointSourceData[i];  
+  delete[] ppPointSourceData;
+  delete pLightMap;
+  delete[] dataFreqDegueux;
+  delete[] dataFreq;
+  delete pNR;
+}
+
+double LightGalaxResol::powSpecDens(double theta, double phi, double freq)
+	 //(Watt/m2/Hz)
+{     
+  long pixelNo= pLightMap->PixIndexSph(theta, phi);
+  int_2 NoSource=(*pLightMap)(pixelNo);
+ 
+  if ( NoSource==(nbPointSource+1)) return 0.; 	// Pas de source
+  else 
+  {  
+    float* binSpectre=(ppPointSourceData[NoSource]); 
+    float InterpRes;
+    float InterpResErr;
+    
+    float thisFreq=freq;
+
+	//Ou est la frequence dans le tableau?
+	pNR->hunt(dataFreqDegueux-1,(unsigned long) nbFreq, thisFreq, &LastFreqIndex);
+	if((LastFreqIndex==0)||(LastFreqIndex==nbFreq)) 
+	{	cerr<<"Appel a powSpecDens dans LightGalaxResol hors limites de frequence"<<endl;
+		LastFreqIndex=1;
+		return 0;
+	}
+			
+	
+      // On interpole 3 eme ordre. Num Reciepies fonction.
+    pNR->polint(dataFreqDegueux-1+LastFreqIndex,binSpectre-1+LastFreqIndex,
+    	2,thisFreq,&InterpRes,&InterpResErr);
+      // InterpRes en Jansky:10-26 W/m2/Hz           
+    return (double) (InterpRes*10.e-26); 	// W/m2/Hz
+      // RMQ: Le Jansky est une unite de source ponctuelle !
+      // Pas de steradian dans les homogeneites.
+ }
+}
+
+static double randMax=RAND_MAX;
+double LightGalaxResol::MyRan() {
+	return rand()/randMax;
+}
+
Index: /tags/SigPred_V1/SigPredictor/lightgalaxresol.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightgalaxresol.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightgalaxresol.h	(revision 799)
@@ -0,0 +1,41 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SKYGALAXRESOL_SEEN
+#define SKYGALAXRESOL_SEEN
+
+#include "spheregorski.h"
+#include "numrecipes.h"
+#include "abslightsource.h"
+
+class LightGalaxResol:public AbsLightSource {
+ public:
+    LightGalaxResol(int_4 nside);
+    virtual ~LightGalaxResol();
+    virtual double powSpecDens(double theta,double phi,double freq);
+    // Return power density at coordinates
+    // (W/m2/Hz/St)
+    
+  	
+ // virtual void Reset(const SpectralResponse& Filtre);
+    // stocke le flux de puissance (watt/m2/st) pour toutes les sources dans SourcesInBandPower
+ protected:
+    int_4 nlat;         // Normalement c est 2048!  
+    int_4 nbPixelLight; // nbre de pixel de la carte du ciel 
+    int_4 nbFreq;       // Nbre de frq ou les sources sont talonnes
+    int_4 nbPointSource;// Nbre de sources ponctuelles
+    NumRecipes* pNR;
+    unsigned long LastFreqIndex;
+
+    r_8* dataFreq;     
+    // Tableau des frequences auxquelles sont calcules les cartes guiderdonni
+    float* dataFreqDegueux;
+    // Idem. Ne sert que pour appeller la fonction NR Polint
+    r_4** ppPointSourceData;    // (Jansky, cad 10-26 W/m2/st/Hz)
+    //donnes sur chaque source a chaque frequence
+    SphereGorski <uint_2>* pLightMap;  
+    //grand (!) tableau des pixels du ciel. Contient eventuellement l index de 
+    // du tableau SourcesInBandPower, cad le no de la source.
+    double MyRan();
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.cc	(revision 799)
@@ -0,0 +1,38 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "lightsrcmappowerinband.h"
+#include "fitsioserver.h" 
+
+double LightSrcMapPowerInband::powerDensAmpli(double theta,double phi)
+{	return (*pMap)(theta,phi);
+  		// Return power density Amplidude at coordinates
+    	// (W/m2/Hz/st)
+}
+
+double LightSrcMapPowerInband::spectre(double freq)
+{	if((freq>freqMax)||(freq<freqMin)) return 0.;
+    	else return SpectreAmpli;
+}
+
+LightSrcMapPowerInband::LightSrcMapPowerInband
+	(char filename[],long nside, double MinFreq, double MaxFreq)
+  :freqMax(MaxFreq), freqMin(MinFreq), nlat(nside)
+{  
+   resolution=2*3.14159/4./nlat; // Resolution angulaire de la carte lue pour ce fond.
+                      			 // Radians
+   SpectreAmpli=1./(freqMax-freqMin);
+   if(MinFreq>MaxFreq) cerr<<"Erreur freq Min>freq Max in LightSrcMapPowerInband"<<endl;
+   
+   try {	pMap=new SphereGorski<float> (nlat);	}
+   catch(bad_alloc) {
+      cerr<<"Memory booking error LightSrcMapPowerInband, nlat= "<<nlat<<endl;
+      exit(-1); 
+   }
+   cout<<"On lit la carte "<<filename<<endl;  
+   FitsIoServer FitsServer;
+   FitsServer.load((*pMap),filename);
+}
+
+LightSrcMapPowerInband::~LightSrcMapPowerInband()
+{	delete pMap;
+}
Index: /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightsrcmappowerinband.h	(revision 799)
@@ -0,0 +1,24 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LightSsrcMap_H
+#define LightSsrcMap_H
+#include "spheregorski.h"
+#include "lightsrcsep.h"
+
+
+class LightSrcMapPowerInband:public LightSrcSep {
+  public: 
+	LightSrcMapPowerInband(char filename[],long nside, double freqMin, double freqMax);
+	virtual ~LightSrcMapPowerInband();
+  protected:
+	int nlat;
+	double freqMax;
+	double freqMin;
+ 	SphereGorski<float>* pMap;
+ 	double SpectreAmpli;
+ 	
+ 	virtual double powerDensAmpli(double theta,double phi);
+ 	virtual double spectre(double freq);
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightsrcsep.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightsrcsep.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightsrcsep.h	(revision 799)
@@ -0,0 +1,57 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LightSrcSep_H
+#define LightSrcSep_H
+
+#include "abslightsource.h"
+
+class LightSrcSep:public AbsLightSource {
+ friend class LSrcRadSpecSep;
+public: 
+  LightSrcSep(){LScrFSep=true; QPtSource=false;}
+  ~LightSrcSep() {	}
+  virtual double powSpecDens(double theta,double phi,double freq){
+    // Return power density at coordinates (W/m2/Hz/St)
+  		return spectre(freq)*powerDensAmpli(theta,phi);
+  }
+  
+  virtual double spectre(double frequence)=0; 
+        // Return the frequency spectrum independent of coordinates
+        // should have no dimension
+  virtual double powerDensAmpli(double theta,double phi)=0;
+  		// Return power density Amplidude at coordinates
+    	// (W/m2/Hz/St)
+protected:
+};
+
+class LightSrcDiffBlackBody:public LightSrcSep {
+	// Classe mettant en oeuvre les sources parametrees 
+	// Comme des ecarts de temperature a la tempe du CMB
+	// valable dans la limite des petits ecarts de tempe de TCMB
+    // DeltaT retourne la DTCMB apparente au point theta, phi.
+public:
+   LightSrcDiffBlackBody() {
+   		h_kT=6.626e-34/T_CMBR/1.38e-23;  // Hz-1
+    	Prefac=2*6.626e-34*6.626e-34/(9.e16*1.38e-23*T_CMBR*T_CMBR); 
+    	     // 2*h2/(c2*k*TCMB2) ( W/m2/st/Hz^3/K) ou (J/m2/st/Hz^4/K)
+   }
+   virtual ~LightSrcDiffBlackBody() { };
+protected:
+	double h_kT;   		// nanosecondes
+	double Prefac;
+	virtual double DeltaT(double theta,double phi)=0;  	// en Kelvin
+    virtual double spectre(double freq){		   		//	Frequence en Hz
+		return pow(freq,4)*exp(h_kT*freq)/pow((exp(h_kT*freq)-1),2);
+		// Hz^4  
+	}	
+ 	virtual double powerDensAmpli(double theta,double phi){
+  		// Return power density Amplidude at coordinates
+    	// J/m2/st/Hz^4 ou W/m2/St/Hz^5
+       	double PDens= Prefac*DeltaT(theta,phi);
+  		return PDens;
+    }
+
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lightsynchro.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightsynchro.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightsynchro.cc	(revision 799)
@@ -0,0 +1,44 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "spheregorski.h"
+#include "fitsioserver.h"
+#include "lightsynchro.h"
+
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "macenvvariables.h"
+#endif
+
+LightSynchro::LightSynchro(int_4 nside) 
+{ sprintf(Name, "Lumire Synchrotron");
+  freqRef=0.408e9;				// frequence de la carte 0.408 GHz
+  index=-0.9;
+  nlat=nside;
+  if (!(nside==1024)) {
+    cerr<<" Fond LightSynchro calcul avec Nside= "<<nside<<endl;
+    cerr<<" carte de resolution nside= 1024 disponible"<<endl;
+  } 
+  
+  resolution=2*3.14159/4./nlat;		// Radians
+  
+  try {	pMap=new SphereGorski<float>(nlat);  }	// MJy/m2/Hz/st cad 10-20 W/m2/Hz/st
+  catch(bad_alloc) {
+    cerr<<"memory booking error LightSynchro constructor nlat= "<<nlat<<endl;
+    exit(-1);
+  }
+  char filename[150];
+#ifndef __MWERKS__
+  char *PATHDataLScr=getenv("PATHDataLScr"); 
+#endif
+  sprintf(filename,"%ssync_res%04i.fits",PATHDataLScr,nside);
+  FitsIoServer FitsServer;
+  FitsServer.load((*pMap),filename);
+    // unite: MJy/m2/Hz/st. 
+};
+
+LightSynchro::~LightSynchro() {
+  delete pMap;  
+}
Index: /tags/SigPred_V1/SigPredictor/lightsynchro.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lightsynchro.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lightsynchro.h	(revision 799)
@@ -0,0 +1,35 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SKYSYNCHRO_SEEN
+#define SKYSYNCHRO_SEEN
+
+
+#include "lightsrcsep.h"
+#include "spheregorski.h"
+
+
+class LightSynchro:public LightSrcSep{
+ public:
+ LightSynchro(int_4 nside=256);
+ virtual ~LightSynchro();
+ protected:
+ int_4 nlat;
+ double index;
+ double freqRef; 
+ SphereGorski<float> * pMap; // Map units unknown to be looked for
+ 
+ virtual double spectre(double freq) 
+ {
+  	return pow(freq/freqRef,index); 		// No dimension
+ }
+ virtual double powerDensAmpli(double theta,double phi)
+ {
+   return 2*1.38e-23/(3.e8*3.e8)*freqRef*freqRef*(*pMap)(theta,phi);
+  // return 1.e-20*(*pMap)(theta,phi); 	// Ainsi parlait Haslam par Maitre R. Teyssier
+      	// (W/m2/Hz/St)  soit du MJ/St	
+ }
+
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lobecartographmoyen.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobecartographmoyen.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobecartographmoyen.cc	(revision 799)
@@ -0,0 +1,64 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+   #include "macenvvariables.h"
+#endif
+#include "unitvector.h"
+#include "lobecartographmoyen.h"
+
+LobeCartoMoyen::LobeCartoMoyen(CarteLobe* pcart,double freqMin, double freqMax)
+  :MeanFreqLobe(freqMin,freqMax), pCarte(pcart) {
+  
+  	angleMax=pCarte->halfAngAperture();
+  	cosanglemax=cos(angleMax);
+  	sprintf(Name,"Lobe Cartographi TICRA");
+}
+
+LobeCartoMoyen::~LobeCartoMoyen(){	}
+double LobeCartoMoyen::lobeResol() const{
+  // Du fait des extrapolations paraboliques, on peut calculer le lobe en
+  // tout point sans reconter de probleme de precision
+  return pCarte->resol();
+}
+
+UnitVector LobeCartoMoyen::VecShift(const UnitVector& Vpointe, const UnitVector& VY) const {
+  double alpha,beta;
+  pCarte->UVToAng(pCarte->UCentre(),pCarte->VCentre(),alpha,beta);
+  UnitVector Vec;
+  Vec=Vpointe.Rotate(VY,alpha);
+  Vec=Vec.Rotate(Vpointe,beta);
+  return Vec;
+}
+
+
+double LobeCartoMoyen::weigthAmpl(const UnitVector& Vcourant,const UnitVector& VP, 
+  		const UnitVector& VY) const{
+  UnitVector VCentreLobe;
+  VCentreLobe=VecShift(VP, VY);
+  double cosbidon=(Vcourant*VCentreLobe);
+  if (cosanglemax>cosbidon) return 0.;
+  else 
+  //return 1.;	//Debug
+  {   double cosalpha=Vcourant*VP;
+  	//  cout<<1-cosalpha*cosalpha<<endl;
+ 	  double sinalpha=sqrt(1-cosalpha*cosalpha); 	// Approx. des petits angles
+  
+	  Vector3d VProjVcourant;
+	  VProjVcourant=Vcourant-cosalpha*VP;
+	  VProjVcourant.Normalize();
+	  
+	  UnitVector VX=VY^VP;		
+	  // Produit vectoriel
+	  // On peut optimiser en sortant le calcul de VX de la fonction.
+	  double cosbeta=VProjVcourant*VX;
+	  double sinbeta=VProjVcourant*VY;
+
+	  return pCarte->Value(sinalpha*cosbeta,sinalpha*sinbeta);
+	  // Il faut comprendre les homogeneites de la carte BUGGGGGG  XXXXXXX
+	 
+  }
+  
+}
Index: /tags/SigPred_V1/SigPredictor/lobecartographmoyen.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobecartographmoyen.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobecartographmoyen.h	(revision 799)
@@ -0,0 +1,30 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LOBE_CARTO_SEEN
+#define LOBE_CARTO_SEEN
+
+
+#include "unitvector.h"
+#include "meanlobe.h"
+#include "cartelobe.h"
+
+
+
+class LobeCartoMoyen:public MeanFreqLobe{
+ public:
+    LobeCartoMoyen(CarteLobe* pcart,double freqMin, double freqMax);
+    virtual ~LobeCartoMoyen();
+    virtual double lobeResol() const;
+ private:
+    CarteLobe* pCarte;
+    virtual UnitVector VecShift(const UnitVector& Vpointe,const UnitVector& VY)const;
+    virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  		const UnitVector& VY) const;
+    // return relative weigth of incoming power from direction Vcourant
+
+    // Computes adequate map resolution in radian for this lobe to prevent loss info 
+     virtual double ResolutionCurve(double angle) const { return 1.;}
+     
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lobes4pi.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobes4pi.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobes4pi.h	(revision 799)
@@ -0,0 +1,37 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+#include "meanlobe.h"
+
+
+
+class Lobe4PiGaussien:public MeanFreqLobe 
+{ 
+ public:
+    Lobe4PiGaussien(double FWHM_degre, double freqMin, double freqMax)
+      :MeanFreqLobe(freqMin,freqMax){
+		sigma=minToRad(FWHM_degre/2.354*60.);
+		angleMax=M_PI;
+		cosanglemax=cos(angleMax);
+		sprintf(Name, "Lobe 4Pi Gaussien");
+    }
+    virtual ~Lobe4PiGaussien(){ };
+   	virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  	  const UnitVector& VY) const {
+		double ang= acos(VP*VInteg);
+		return exp(-ang*ang/(2*sigma*sigma));
+	}
+	virtual double lobeResol() const{
+		return sigma*2.354/3.; // On veut pixeliser/calculer  au moins FWHM/3.
+	} 
+	virtual double ResolutionCurve(double angleShift) const{
+		return (double) (1.+(angleShift/sigma/7.));
+	}
+	
+ protected: 
+    double sigma;
+};
Index: /tags/SigPred_V1/SigPredictor/lobesep.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobesep.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobesep.h	(revision 799)
@@ -0,0 +1,27 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LOBESEP_H
+#define LOBESEP_H
+
+#include "abslobenopolar.h"
+
+class LobeSep:public AbsLobeNoPolar {
+public:
+  LobeSep(double freqMin, double freqMax):AbsLobeNoPolar(freqMin,freqMax) {  LobeFSep=true;}
+  ~LobeSep() {	}
+  virtual double weigth(const UnitVector& VInteg, const UnitVector& VP, 
+  	const UnitVector& VY, double freq) const 
+   	// Return relative weigth of incoming power from direction VInteg
+    // at frequency freq.
+  {  return spectre(freq)*weigthAmpl(VInteg,VP,VY);}
+  virtual double spectre(double freq) const=0;
+  virtual double weigthAmpl(const UnitVector& VInteg, const UnitVector& VP, 
+  	const UnitVector& VY) const=0;
+protected:
+
+    // Return relative weigth of incoming power from direction VInteg
+    // Averaged on frequency. Mean Lobe
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/lobespecresp.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobespecresp.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobespecresp.cc	(revision 799)
@@ -0,0 +1,9 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "abslobenopolar.h"
+#include "lobespecresp.h"
+
+
+double LobeSpecResp::transmission(double freq) const {
+  	return pLobe->weigth(VIntegration,Vpointe,VaxeY,freq);
+}
Index: /tags/SigPred_V1/SigPredictor/lobespecresp.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lobespecresp.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lobespecresp.h	(revision 799)
@@ -0,0 +1,27 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LOBESPECRESP_H
+#define LOBESPECRESP_H
+
+#include "specresp.h"
+#include "unitvector.h"
+
+class AbsLobeNoPolar;
+
+
+class LobeSpecResp:public SpectralResponse {
+  public:	
+	LobeSpecResp(AbsLobeNoPolar* pLob, const UnitVector& VInteg, 
+  	  const UnitVector& VP,const UnitVector& VY, double numin=1.,
+  	  double numax=1000.)
+		:SpectralResponse(numin,numax),pLobe(pLob),Vpointe(VP),
+		 VaxeY(VY), VIntegration(VInteg) {	}
+	~LobeSpecResp() { }	 
+    virtual double transmission(double freq) const ;
+  protected:	
+	AbsLobeNoPolar* pLobe;
+	const UnitVector& Vpointe;
+	const UnitVector& VaxeY;
+	const UnitVector& VIntegration;
+};
+#endif
Index: /tags/SigPred_V1/SigPredictor/lsrcradspec.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/lsrcradspec.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lsrcradspec.cc	(revision 799)
@@ -0,0 +1,16 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include "lsrcradspec.h"
+#include "abslightsource.h"
+
+LSrcRadSpec::LSrcRadSpec(AbsLightSource* pLSrc, double theta, double phi, 
+ double FreqMin, double FreqMax) :RadSpectra(FreqMin,FreqMax)
+ {	pLightSrc=pLSrc;
+ 	Theta=theta;
+ 	Phi=phi;
+ }
+ 
+ double LSrcRadSpec::flux(double freq) const {
+	return (pLightSrc->powSpecDens(Theta,Phi,freq));
+}
+
Index: /tags/SigPred_V1/SigPredictor/lsrcradspec.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/lsrcradspec.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/lsrcradspec.h	(revision 799)
@@ -0,0 +1,30 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef LSrcRadSpec_H
+#define LSrcRadSpec_H
+
+#define FreqIntegMin (1.)
+#define FreqIntegMax (1000.)
+
+#include "radspec.h"
+
+class AbsLightSource;
+
+class  LSrcRadSpec:public RadSpectra{
+public: 
+//---------------	Constructor ---------------------------------
+	LSrcRadSpec(AbsLightSource* pLSrc, double theta, double phi, 
+		 double FreqMin=FreqIntegMin, double FreqMax=FreqIntegMax);
+	   
+//---------------   Desctructor ---------------------------------
+	~LSrcRadSpec() {}
+	
+    double flux(double freq) const;
+  
+protected:
+	AbsLightSource* pLightSrc;
+	double Theta;
+	double Phi;
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/macenvvariables.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/macenvvariables.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/macenvvariables.h	(revision 799)
@@ -0,0 +1,33 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef MacEnvVariables_H
+#define MacEnvVariables_H
+
+/*
+//	Version portable
+#define PATHDataLScr "Stock:DataDock:ArchSigPredData:templates:"
+#define PATHCarteLobe "Stock:DataDock:ArchSigPredData:ApparatusData:CartesDeLobe:"
+#define PATHSTangoRes "Stock:DataDock:ArchSigPredData:SuperTangoProg:observations:"
+#define PATHResults "Stock:DataDock:ArchSigPredData:Results:"
+#define HOME "Stock:DataDock:ArchSigPredData:"
+#define PathTimeline "Stock:DataDock:ArchSigPredData:ApparatusData:Trapani99:"
+*/
+
+
+//*
+// Version Mac Labo
+#define PATHDataLScr "Amarylis:ArchSigPredData:templates:"
+#define PATHCarteLobe "Amarylis:ArchSigPredData:CartesDeLobe:"
+#define PATHSTangoRes "Amarylis:ArchSigPredData:SuperTangoProg:observations:"
+#define PATHResults "Amarylis:ArchSigPredData:Results:"
+#define HOME "Amarylis:ArchSigPredData:"
+#define PathInstrumentData "Amarylis:ArchSigPredData:ApparatusData:ArcheopsTrapani99:"
+#define PathTimeline "Amarylis:Trapani99:"
+// */
+
+# define ERROR_LABEL -32768
+
+//#define SetNlatLSource (1024)		// version de production
+#define SetNlatLSource (256)		// version debug
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/makebolotimeline.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/makebolotimeline.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/makebolotimeline.cc	(revision 799)
@@ -0,0 +1,99 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <iomanip.h>
+#ifdef __MWERKS__
+   #include "macenvvariables.h"
+#else
+	# define ERROR_LABEL -32768
+#endif
+
+#include "makebolotimeline.h"
+
+
+MakeBoloTimeline::MakeBoloTimeline(char file[], AbsLightSource& LSrc, 
+		AbsLobeNoPolar& Lobe, SpectralResponse &Filter, Normalisation Norm)
+	:tool(&LSrc,&Lobe,&Filter), NormMethode(Norm)
+{
+#ifndef __MWERKS__
+  char* PATHResults=getenv("PATHResults"); 
+#endif
+	char filecur[150];
+// On remplit le ReadMe
+	sprintf(filecur,"%s%s.ReadMe", PATHResults,file);
+	cout<<filecur<<endl;
+	ofstream ReadMe(filecur,ios::out|ios::trunc);
+	
+	ReadMe<<" Fichier de timeline pour:"<<endl;
+	ReadMe<<" La source physique: "<<endl;
+	LSrc.print(ReadMe);
+	ReadMe<<" Le Lobe : "<<endl;
+	Lobe.print(ReadMe);
+	ReadMe<<" Le Filtre: "<<endl;
+	Filter.Print(ReadMe);	
+	ReadMe<<endl;
+	
+	NormRJ=tool.NormKelvinRayleighJeans();
+	NormTCMB=tool.NormKelvinCMB();
+	
+//	ReadMe.setf(ios::scientific);
+	ReadMe<< setprecision(3);
+	ReadMe<<" Constantes de normalisation"<<endl;
+	ReadMe<<" RayleighJeans: "<<NormRJ<< " KelvinRJ/(W/m2)"<<endl;
+	ReadMe<<" TempeCMB: "<<NormTCMB<< " Kelvin CMB/(W/m2)"<<endl;
+	
+	ReadMe<<"Cette carte est en unit: ";
+	switch(NormMethode)
+	{	
+		case TempeCMB: 
+		{	ReadMe<< "Kelvin CMB"<<endl;
+			break;
+		}
+		case RayleighJeans:
+		{ 	ReadMe<< "Kelvin RayleighJeans"<<endl;
+			break;
+		}
+		default: ReadMe<< "Watt/m2 effectif de mirroir"<<endl;
+	}
+	ReadMe<<"Fichier ecrit selon le format"<<endl;
+	ReadMe<<"Radian\tRadian\tAbove Unit"<<endl;
+	
+	// J'espere que c'est assez 
+	ReadMe.close();
+	
+	sprintf(filecur,"%s%s.dat", PATHResults,file);
+	
+	pMystr= new ofstream(filecur,ios::out|ios::trunc);
+//	pMystr->setf(ios::scientific);
+	(*pMystr)<<setprecision(5);
+}
+
+void MakeBoloTimeline::addToStreamArchTOI(double theta,double phi)
+{	//	 Theta en heures, Phi en degres
+	if((theta==ERROR_LABEL)||(phi==ERROR_LABEL)) return;
+	double thetaRad=theta*15./180*M_PI;
+	double phiRad=phi/180.*M_PI;
+	
+	addToStream(theta,phi);
+}
+
+void MakeBoloTimeline::addToStream(double theta,double phi)
+{	// theta et Phi en Radian.	
+	switch (NormMethode)
+	{
+		case RayleighJeans:
+		{	double TempeRJ= NormRJ*tool.compPixel(theta,phi);
+			(*pMystr)<<theta<<'\t'<<phi<<'\t'<<TempeRJ<<'\n';
+			return;
+		}
+		case TempeCMB:
+		{	double tempeCMB= NormTCMB*tool.compPixel(theta,phi);
+			(*pMystr)<<theta<<'\t'<<phi<<'\t'<<tempeCMB<<'\n';
+			return;
+		}
+		default:
+		{	double power=tool.compPixel(theta,phi);
+			(*pMystr)<<theta<<'\t'<<phi<<'\t'<<power<<'\n';
+			return;
+		}
+	}
+}
Index: /tags/SigPred_V1/SigPredictor/makebolotimeline.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/makebolotimeline.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/makebolotimeline.h	(revision 799)
@@ -0,0 +1,37 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef MakeBoloTimeline_H
+#define MakeBoloTimeline_H
+
+#include <fstream>
+
+#include "abslightsource.h"
+#include "abslobenopolar.h"
+#include "sigcalctools.h"
+#include "specresp.h"
+
+enum Normalisation { RayleighJeans, TempeCMB, Power};
+class MakeBoloTimeline
+{	
+public:
+	MakeBoloTimeline(char file[], AbsLightSource& LSrc, 
+	AbsLobeNoPolar& Lobe, SpectralResponse &Filter, Normalisation Norm);
+	~MakeBoloTimeline() 
+	{	pMystr->close();
+		delete pMystr;
+	}	
+	void addToStreamArchTOI(double theta, double phi);
+		//	 Theta en heures, Phi en degres
+	void addToStream(double theta,double phi);
+		// theta et phi en radian
+protected:
+	SigCalcTool tool;
+	ofstream* pMystr;
+	double NormRJ;
+	double NormTCMB;
+	Normalisation NormMethode;
+//	long flushControl;
+
+};
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/maketimeline.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/maketimeline.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/maketimeline.cc	(revision 799)
@@ -0,0 +1,294 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include "makebolotimeline.h"
+#include "bruit.h"
+//#include "fitsio.h"
+
+#ifdef __MWERKS__
+   #include "macenvvariables.h"
+   #define SetNlatLSource (256)		// version debug
+   #define NbSAMPLE (2000)			//	Pour le debug
+#else
+   # define ERROR_LABEL -32768
+   #define SetNlatLSource (1024)		// version de production CXX
+   #define NbSAMPLE (2684300)
+#endif
+
+
+#include "alllobe.h"
+#include "alllightsources.h"
+#include "sigcalctools.h"
+#include "allfilter.h"
+#include "fitsioserver.h"
+
+
+
+int MakeArcheopsTimelines(char FilePointes[]);
+int GenTimLin1_F(long Nbsample, double TStep, float fknee);
+
+int MakeArcheopsTimelines(char FilePointes[])
+{
+
+/*
+
+// On genre des timelines de bruit en 1/f
+ 	GenTimLin1_F(NbSAMPLE,5.83e-3,1.);
+ 	GenTimLin1_F(NbSAMPLE,5.83e-3,0.5);
+ 	GenTimLin1_F(NbSAMPLE,5.83e-3,0.1);
+*/
+ 	
+// On passe aux cartes de fonds physiques
+// Initialisations
+   FitsIoServer FitsServer;
+   char filename[150]="";
+     
+//	cerr.setf(ios::scientific);
+  	cerr<< "C'est parti"<<endl;
+  	
+ 	// Des Lobes
+ 	LobeGaussien GLobeTrapChan1(8./60.,75.e9,600.e9);  //degrs, freqmin, Freqmax
+ 	LobeGaussien GLobeTrapChan2(8./60.,75.e9,600.e9); 
+ 	LobeGaussien GLobeTrapChan3(8./60.,75.e9,600.e9); 
+ 	
+#ifndef __MWERKS__
+/* char* PATHResults=getenv("PATHResults"); 
+   char* PATHDataLScr=getenv("PATHDataLScr");   
+   char* PATHCarteLobe=getenv("PATHCarteLobe");
+*/
+  char* PathInstrumentData=getenv("PathInstrumentData");
+#endif
+
+ 	// Des Filtres
+ 	sprintf(filename, "%sfiltretrap99chan1.dat",PathInstrumentData);
+ 	TabulatedFilter FiltreTrapChan1(filename); 	
+ 	sprintf(filename, "%sfiltretrap99chan2.dat",PathInstrumentData);
+ 	TabulatedFilter FiltreTrapChan2(filename);
+ 	sprintf(filename, "%sfiltretrap99chan3.dat",PathInstrumentData);
+ 	TabulatedFilter FiltreTrapChan3(filename);
+
+/*	
+ 	// On imprime les filtres
+ 	ofstream Cfilter1("TraceFiltre1", ios::out);
+	FiltreTrapChan1.TraceFiltre(Cfilter1);
+	ofstream Cfilter2("TraceFiltre2", ios::out);
+	FiltreTrapChan2.TraceFiltre(Cfilter2);
+ 	ofstream Cfilter3("TraceFiltre3", ios::out);
+	FiltreTrapChan3.TraceFiltre(Cfilter3);
+*/
+
+// /*
+ 	// Les Processus Physiques
+ 	// LightSynchro Synchro(SetNlatLSource);		// Pas pour Archeops
+ 	LightDiffDust DiffDust(SetNlatLSource);
+ 	LightCMBPrim CMBPrim(SetNlatLSource);		// Pour debuter
+ 	// Le bruit reste  chacun de le gnerer
+ 	LightDipole CMBDipole;
+ 	
+ 	// On definit les timeline desirees
+ 
+ 	char fileOut[150];	
+ 	
+ 	// Bolo 1_1 soit bolo11.
+ 	sprintf(fileOut, "ASigPred_bolo1_1Dipole");
+ 	MakeBoloTimeline TL1_1Dipole(fileOut,CMBDipole,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+
+ 	sprintf(fileOut, "ASigPred_bolo1_1CMBPrim");
+ 	MakeBoloTimeline TL1_1CMBPrim(fileOut, CMBPrim,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);	
+
+ 	sprintf(fileOut, "ASigPred_bolo1_1DiffDust");
+ 	MakeBoloTimeline TL1_1DiffDust(fileOut,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+
+// Bolo 1_3 soit bolo13.
+ 	sprintf(fileOut, "ASigPred_bolo1_3Dipole");
+ 	MakeBoloTimeline TL1_3Dipole(fileOut, CMBDipole,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+
+ 	sprintf(fileOut, "ASigPred_bolo1_3CMBPrim");
+ 	MakeBoloTimeline TL1_3CMBPrim(fileOut, CMBPrim,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);	
+
+ 	sprintf(fileOut, "ASigPred_bolo1_3DiffDust");
+ 	MakeBoloTimeline TL1_3DiffDust(fileOut,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 		
+ 	// Bolo 2_4 soit bolo9.
+ 	sprintf(fileOut, "ASigPred_bolo2_4Dipole");
+ 	MakeBoloTimeline TL2_4Dipole(fileOut, CMBDipole,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+
+ 	sprintf(fileOut, "ASigPred_bolo2_4CMBPrim");
+ 	MakeBoloTimeline TL2_4CMBPrim(fileOut, CMBPrim,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);	
+ 		
+ 	sprintf(fileOut, "ASigPred_bolo2_4DiffDust");
+ 	MakeBoloTimeline TL2_4DiffDust(fileOut,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 		
+ 	// Bolo 2_5 soit bolo4.
+ 	sprintf(fileOut, "ASigPred_bolo2_5Dipole");
+ 	MakeBoloTimeline TL2_5Dipole(fileOut, CMBDipole,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+
+ 	sprintf(fileOut, "ASigPred_bolo2_5CMBPrim");
+ 	MakeBoloTimeline TL2_5CMBPrim(fileOut, CMBPrim,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);	
+
+ 	sprintf(fileOut, "ASigPred_bolo2_5DiffDust");
+ 	MakeBoloTimeline TL2_5DiffDust(fileOut,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 				
+ 	// Bolo 3_6 soit bolo15.
+ 	sprintf(fileOut, "ASigPred_bolo3_6Dipole");
+ 	MakeBoloTimeline TL3_6Dipole(fileOut, CMBDipole,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 	
+	sprintf(fileOut, "ASigPred_bolo3_6CMBPrim");
+ 	MakeBoloTimeline TL3_6CMBPrim(fileOut, CMBPrim,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);	
+
+ 	sprintf(fileOut, "ASigPred_bolo3_6DiffDust");
+ 	MakeBoloTimeline TL3_6DiffDust(fileOut,DiffDust,GLobeTrapChan1,
+ 		FiltreTrapChan1,RayleighJeans);
+ 		
+ 		
+ 	// Il faut maintenant lire le fichier des pointages
+ 	
+#ifndef __MWERKS__
+  char* PathTimeline=getenv("PathTimeline");
+#endif
+
+	char FileName[150];
+	sprintf(FileName, "%s%s",PathTimeline,FilePointes);
+	ifstream InFile(FileName);
+	
+#ifdef __MWERKS__	
+	// Cxx est trop con!
+	if (!InFile.is_open()) 
+#else
+	if (!InFile)
+#endif 
+	{	cerr<< "Erreur a l'ouverture du fichier de donnees :"<<FileName<< endl;
+		exit(-1);
+	}
+
+ 	long SamplNumber;
+ 	float alphaBolo2_5,deltaBolo2_5, alphaBolo2_4,deltaBolo2_4,
+ 		alphaBolo1_1,deltaBolo1_1, alphaBolo1_3,deltaBolo1_3, 
+ 		alphaBolo3_6, deltaBolo3_6; 
+ 	char Line[200];
+ 	
+ 	int nbParamLu;
+ 	long compteur=0;
+//	On genre les timelines
+	while(!InFile.eof())
+	{	
+	// On lit les directions pointes
+		InFile.getline(Line,200);
+		nbParamLu=sscanf(Line," %d %f %f %f %f %f %f %f %f %f %f",
+			&SamplNumber, &alphaBolo2_5,&deltaBolo2_5, &alphaBolo2_4,&deltaBolo2_4,
+ 		&alphaBolo1_1,&deltaBolo1_1, &alphaBolo1_3,&deltaBolo1_3, 
+ 		&alphaBolo3_6, &deltaBolo3_6);
+ 		if (nbParamLu!=11) {
+ 			cerr<<" Problme a la lecture du fichier: "<<FileName<<endl;
+ 			cerr<<" Ligne :"<<compteur<<endl;
+ 		}
+ 		
+ 	// On remplit les fichiers correspondants 
+ 		TL1_1Dipole.addToStreamArchTOI(alphaBolo1_1,deltaBolo1_1);
+ 		TL1_1CMBPrim.addToStreamArchTOI(alphaBolo1_1,deltaBolo1_1);
+ 		TL1_1DiffDust.addToStreamArchTOI(alphaBolo1_1,deltaBolo1_1);
+ 		
+ 		TL1_3Dipole.addToStreamArchTOI(alphaBolo1_3,deltaBolo1_3);
+		TL1_3CMBPrim.addToStreamArchTOI(alphaBolo1_3,deltaBolo1_3);
+ 		TL1_3DiffDust.addToStreamArchTOI(alphaBolo1_3,deltaBolo1_3);
+ 		
+ 		TL2_4Dipole.addToStreamArchTOI(alphaBolo2_4,deltaBolo2_4);
+ 		TL2_4CMBPrim.addToStreamArchTOI(alphaBolo2_4,deltaBolo2_4);
+ 		TL2_4DiffDust.addToStreamArchTOI(alphaBolo2_4,deltaBolo2_4);
+ 		
+ 		TL2_5Dipole.addToStreamArchTOI(alphaBolo2_5,deltaBolo2_5);
+		TL2_5CMBPrim.addToStreamArchTOI(alphaBolo2_5,deltaBolo2_5);
+ 		TL2_5DiffDust.addToStreamArchTOI(alphaBolo2_5,deltaBolo2_5);
+ 		
+ 		TL3_6Dipole.addToStreamArchTOI(alphaBolo3_6,deltaBolo3_6);
+	 	TL3_6CMBPrim.addToStreamArchTOI(alphaBolo3_6,deltaBolo3_6);
+ 		TL3_6DiffDust.addToStreamArchTOI(alphaBolo3_6,deltaBolo3_6);
+ 		
+ 		compteur++;
+ 		if((compteur%1000)==0)
+ 		{	double pourcent=100.*compteur/NbSAMPLE;
+ 			cout<<"Compteur= "<<compteur<<endl;
+ 			cout<<"Travail effectue a :"<<pourcent<<" %"<<endl;
+ 		}
+	}
+	
+	#ifndef __MWERKS__
+	char* PATHResults=getenv("PATHResults"); 
+#endif
+
+	char filecur[150];
+	
+// On remplit le ReadMe
+	sprintf(filecur,"%sArchSigPredTimeLineReadMe",PATHResults);
+	ofstream ReadMe(filecur, ios::out|ios::trunc);
+	
+	ReadMe<<"Files containning timeline stream"<<endl;
+	ReadMe<<"Generated using ArchSignalPredictor, D. Yvon /03/2000"<<endl;
+	ReadMe<<"Templates used: as provided by R. Teyssier"<<endl;
+	ReadMe<<"Resolution of computation: "<< SetNlatLSource <<" Nlat Healpix"<< endl;
+	ReadMe<<"Number of consecutive samples generated:"<<compteur<<endl;
+	ReadMe<<"Samples with undefined pointing skipped."<<endl;
+	ReadMe<<"Dominique Yvon, DAPNIA/SPP, 03/2000"<<endl;
+	ReadMe.close();	// J'espere que c'est assez 
+	
+ 	return 0;
+}
+
+
+int GenTimLin1_F(long Nbsample, double TStep, float fknee)
+{
+#ifndef __MWERKS__
+	char* PATHResults=getenv("PATHResults"); 
+#endif
+
+	char filecur[150];
+	
+// On remplit le ReadMe
+	sprintf(filecur,"%sTimLineNoisetot%2.2fReadMe", PATHResults,fknee);
+	ofstream ReadMe(filecur, ios::out|ios::trunc);
+	
+	ReadMe<<"File containning noise timeline"<<endl;
+	ReadMe<<"Generated using SumNoise et Memnoise of F. Couchot and R. Ansari"<<endl;
+	ReadMe<<"Generation of a white noise of sigma=1."<<" Any unit"<<endl;
+	ReadMe<<"1/f noise added with knee frequency :"<< fknee<<" Hz"<< endl;
+	ReadMe<<"Timestep for generation: "<<TStep<<" Second"<<endl;
+	ReadMe<<"Number of consecutive samples generated:"<<Nbsample<<endl;
+	ReadMe<<"Dominique Yvon, DAPNIA/SPP, 03/2000"<<endl;
+	
+	ReadMe.close();	// J'espere que c'est assez 
+	
+// On bosse
+	SumNoise NoiseGen(fknee);
+
+	TVector<float> Array(Nbsample);
+	cout<< " vous avez reserv un espace memoire de :"<<Nbsample*sizeof(float)<<endl;
+	for(int i=0; i<Nbsample; i++) 
+	{	Array(i)= NoiseGen.Noise();
+		if((i%1000)==0)	{
+			cout<<"GenTimLin1_F, i= "<< i <<endl;
+		}
+	}
+	
+// Write my fitFile, Version avec FitsIOServer
+	FitsIoServer FitsServer;
+	sprintf(filecur,"%sTimLineNoisetot%2.2f.fits",PATHResults,fknee);/* name for new FITS file */
+	FitsServer.save(Array,filecur);
+     
+	return 0;
+}
Index: /tags/SigPred_V1/SigPredictor/meanlobe.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/meanlobe.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/meanlobe.h	(revision 799)
@@ -0,0 +1,33 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef MeanFreqLobe_H
+#define MeanFreqLobe_H
+
+#include "lobesep.h"
+#include "distribflatfilter.h"
+
+/* Class design to implement mean lobes: 
+   Lobes you spatial response do no depend on frequency.
+   Thus behavior is completely described as behavior at mean frequency
+*/
+ 
+class MeanFreqLobe:public LobeSep {
+public:
+  MeanFreqLobe(double freqMin, double freqMax):LobeSep(freqMin,freqMax) {
+  	pRespShape=new DistribFlatFilter(freqMin,freqMax);
+  	One_ExcurFreq=1.;
+  }
+  virtual ~MeanFreqLobe() {delete pRespShape;}
+protected:
+	double One_ExcurFreq;		
+ 	virtual double spectre(double freq) const {
+ 		if(freq<FreqMin||freq>FreqMax) return 0.;
+ 		return One_ExcurFreq;
+ 	}
+ 	
+    DistribFlatFilter* pRespShape;
+    		// Units conversion tools
+};
+
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/numrecipes.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/numrecipes.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/numrecipes.cc	(revision 799)
@@ -0,0 +1,439 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+
+#include "numrecipes.h"
+
+// En tete du fichier nrutil.h
+
+static float sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+
+static double dsqrarg;
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
+
+static double dmaxarg1,dmaxarg2;
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
+        (dmaxarg1) : (dmaxarg2))
+
+static double dminarg1,dminarg2;
+#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
+        (dminarg1) : (dminarg2))
+
+static float maxarg1,maxarg2;
+#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
+        (maxarg1) : (maxarg2))
+
+static float minarg1,minarg2;
+#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
+        (minarg1) : (minarg2))
+
+static long lmaxarg1,lmaxarg2;
+#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
+        (lmaxarg1) : (lmaxarg2))
+
+static long lminarg1,lminarg2;
+#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
+        (lminarg1) : (lminarg2))
+
+static int imaxarg1,imaxarg2;
+#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
+        (imaxarg1) : (imaxarg2))
+
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
+        (iminarg1) : (iminarg2))
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+
+
+#define NR_END 1
+#define FREE_ARG char*
+
+
+void NumRecipes::nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+	fprintf(stderr,"Numerical Recipes run-time error...\n");
+	fprintf(stderr,"%s\n",error_text);
+	fprintf(stderr,"...now exiting to system...\n");
+	exit(1);
+}
+
+float* NumRecipes::vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+	float *v;
+
+	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+	if (!v) nrerror("allocation failure in vector()");
+	return v-nl+NR_END;
+}
+
+int* NumRecipes::ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+	int *v;
+
+	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+	if (!v) nrerror("allocation failure in ivector()");
+	return v-nl+NR_END;
+}
+
+unsigned char* NumRecipes::cvector(long nl, long nh)
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+	unsigned char *v;
+
+	v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+	if (!v) nrerror("allocation failure in cvector()");
+	return v-nl+NR_END;
+}
+
+unsigned long * NumRecipes::lvector(long nl, long nh)
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+	unsigned long *v;
+
+	v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
+	if (!v) nrerror("allocation failure in lvector()");
+	return v-nl+NR_END;
+}
+
+double* NumRecipes::dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
+	if (!v) nrerror("allocation failure in dvector()");
+	return v-nl+NR_END;
+}
+
+float ** NumRecipes::matrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double ** NumRecipes::dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+int ** NumRecipes::imatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	int **m;
+
+	/* allocate pointers to rows */
+	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float** NumRecipes::submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl)
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+	float **m;
+
+	/* allocate array of pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in submatrix()");
+	m += NR_END;
+	m -= newrl;
+
+	/* set pointers to rows */
+	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float** NumRecipes::convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in convert_matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* set pointers to rows */
+	m[nrl]=a-ncl;
+	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float *** NumRecipes::f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	float ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+void NumRecipes::free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void NumRecipes::free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void NumRecipes::free_cvector(unsigned char *v, long nl, long nh)
+/* free an unsigned char vector allocated with cvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void NumRecipes::free_lvector(unsigned long *v, long nl, long nh)
+/* free an unsigned long vector allocated with lvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void NumRecipes::free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void NumRecipes::free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
+/* free a float matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void NumRecipes::free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void NumRecipes::free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void NumRecipes::free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a submatrix allocated by submatrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void NumRecipes::free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a matrix allocated by convert_matrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void NumRecipes::free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh)
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+
+//_____________________________ Fichier hunt.c _________________________________________
+
+void  NumRecipes::hunt(float xx[], unsigned long n, float x, unsigned long *jlo)
+{
+	unsigned long jm,jhi,inc;
+	int ascnd;
+
+	ascnd=(xx[n] > xx[1]);
+	if (*jlo <= 0 || *jlo > n) {
+		*jlo=0;
+		jhi=n+1;
+	} else {
+		inc=1;
+		if (x >= xx[*jlo] == ascnd) {
+			if (*jlo == n) return;
+			jhi=(*jlo)+1;
+			while (x >= xx[jhi] == ascnd) {
+				*jlo=jhi;
+				inc += inc;
+				jhi=(*jlo)+inc;
+				if (jhi > n) {
+					jhi=n+1;
+					break;
+				}
+			}
+		} else {
+			if (*jlo == 1) {
+				*jlo=0;
+				return;
+			}
+			jhi=(*jlo)--;
+			while (x < xx[*jlo] == ascnd) {
+				jhi=(*jlo);
+				inc <<= 1;
+				if (inc >= jhi) {
+					*jlo=0;
+					break;
+				}
+				else *jlo=jhi-inc;
+			}
+		}
+	}
+	while (jhi-(*jlo) != 1) {
+		jm=(jhi+(*jlo)) >> 1;
+		if (x > xx[jm] == ascnd)
+			*jlo=jm;
+		else
+			jhi=jm;
+	}
+}
+// __________________________Fichier polint.c___________________________________
+	
+void NumRecipes::polint(float xa[], float ya[], int n, float x, float *y, float *dy)
+{
+	int i,m,ns=1;
+	float den,dif,dift,ho,hp,w;
+	float *c,*d;
+
+	dif=fabs(x-xa[1]);
+	c=vector(1,n);
+	d=vector(1,n);
+	for (i=1;i<=n;i++) {
+		if ( (dift=fabs(x-xa[i])) < dif) {
+			ns=i;
+			dif=dift;
+		}
+		c[i]=ya[i];
+		d[i]=ya[i];
+	}
+	*y=ya[ns--];
+	for (m=1;m<n;m++) {
+		for (i=1;i<=n-m;i++) {
+			ho=xa[i]-x;
+			hp=xa[i+m]-x;
+			w=c[i+1]-d[i];
+			if ( (den=ho-hp) == 0.0) nrerror("Error in routine polint");
+			den=w/den;
+			d[i]=hp*den;
+			c[i]=ho*den;
+		}
+		*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
+	}
+	free_vector(d,1,n);
+	free_vector(c,1,n);
+	
+}
+
+// ________________________Fichier polin2.c_______________________________________
+
+void NumRecipes::polin2(float x1a[], float x2a[], float **ya, int m, int n, float x1,
+	float x2, float *y, float *dy)
+{
+	
+	int j;
+	float *ymtmp;
+
+	ymtmp=vector(1,m);
+	for (j=1;j<=m;j++) {
+		polint(x2a,ya[j],n,x2,&ymtmp[j],dy);
+	}
+	polint(x1a,ymtmp,m,x1,y,dy);
+	free_vector(ymtmp,1,m);
+}
Index: /tags/SigPred_V1/SigPredictor/numrecipes.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/numrecipes.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/numrecipes.h	(revision 799)
@@ -0,0 +1,52 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef Numrecipes_H_
+#define Numrecipes_H_
+
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+
+
+/* 	Cree pour encapsuler les fonctions utiles de la Numerical Recipes
+	Actuellement encapsule: nrtil.h, hunt.c, polin2.c polint.c
+*/
+
+class NumRecipes{
+public:	
+	NumRecipes() {}
+	~NumRecipes() {}
+	void nrerror(char error_text[]);
+	float *vector(long nl, long nh);
+	int *ivector(long nl, long nh);
+	unsigned char *cvector(long nl, long nh);
+	unsigned long *lvector(long nl, long nh);
+	double *dvector(long nl, long nh);
+	float **matrix(long nrl, long nrh, long ncl, long nch);
+	double **dmatrix(long nrl, long nrh, long ncl, long nch);
+	int **imatrix(long nrl, long nrh, long ncl, long nch);
+	float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+		long newrl, long newcl);
+	float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
+	float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+	void free_vector(float *v, long nl, long nh);
+	void free_ivector(int *v, long nl, long nh);
+	void free_cvector(unsigned char *v, long nl, long nh);
+	void free_lvector(unsigned long *v, long nl, long nh);
+	void free_dvector(double *v, long nl, long nh);
+	void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
+	void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
+	void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
+	void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
+	void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
+	void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+		long ndl, long ndh);
+
+	void hunt(float xx[], unsigned long n, float x, unsigned long *jlo);
+  	void polin2(float *xx,float *yy, float** yaa,int,int,float,float,float*,float*);
+  	void polint(float xa[], float ya[], int n, float x, float *y, float *dy);
+ };
+ 
+ #endif
Index: /tags/SigPred_V1/SigPredictor/quasiptsources.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/quasiptsources.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/quasiptsources.cc	(revision 799)
@@ -0,0 +1,58 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+//   #include "unixmac.h"
+#endif
+#include "strutil.h"
+#include "quasiptsources.h" 
+
+
+
+QuasiPtSources::QuasiPtSources(double tempe,double angsiz,double theta,double phi)
+	:AngSize(angsiz),MyTheta(theta),MyPhi(phi),SurfaceTemperature(tempe)
+{
+  QPtSource="true";
+  resolution=angsiz/1.5;
+  h_kT=6.62e-34/1.38e-23/SurfaceTemperature;  				// Dimension 1/Hz
+  Prefac=8*M_PI*6.626e-34/9.e16;		// 8*Pi*h/c2 == J m-2 st-1 hz-3
+  												// Ou W m-2 st-1 Hz-4
+  CosAngSize=cos(AngSize);
+  
+  pVSource= new UnitVector(MyTheta,MyPhi);
+}
+
+QuasiPtSources::~QuasiPtSources(){
+  delete pVSource;
+}
+
+void QuasiPtSources::SetSourceCoor(double theta,double phi){
+  delete pVSource;
+  MyTheta=theta; 
+  MyPhi=phi;
+  pVSource= new UnitVector(MyTheta,MyPhi);
+  return;
+}
+
+UnitVector QuasiPtSources::GetVSrcCenter()
+{  return *pVSource;
+} 
+
+double QuasiPtSources::DeltaT(double theta,double phi) {
+	// sans dimensions ici
+  UnitVector Vec(theta,phi);
+  double test1= Vec*GetVSrcCenter();
+  if(test1>CosAngSize) return 1.;	
+  else return 0;
+}
+
+double QuasiPtSources::spectre(double freq){
+	// sans dimentions? Oui si on imagine que freq est en Hz!!
+	// en Fait Hz^3
+	double returnval=freq*freq*freq/(exp(h_kT*freq)-1.);
+ 	return returnval;
+}
+
Index: /tags/SigPred_V1/SigPredictor/quasiptsources.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/quasiptsources.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/quasiptsources.h	(revision 799)
@@ -0,0 +1,51 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SKYQUASIPTSOURCES_SEEN
+#define SKYQUASIPTSOURCES_SEEN
+
+#include "lightsrcsep.h"
+#include "unitvector.h"
+
+class QuasiPtSources:public LightSrcDiffBlackBody{
+ public:
+    QuasiPtSources(double tempe=T_CMBR,double angsiz=0.,double theta=0.,double phi=0.);
+    // Tempe en Kelvin, angsiz, theta et phi en radian
+    // angsiz, demi angle au sommet.
+    virtual ~QuasiPtSources();
+    void GetSourceCoor(double& theta,double& phi){
+	theta=MyTheta; phi=MyPhi; return;}
+	
+    void SetSourceCoor(double theta,double phi);
+    virtual UnitVector GetVSrcCenter(); 
+    // Retourne le vecteur unitaire pointant sur la source quasiponctuelle
+    double getAngSize() {return AngSize;}
+ protected:
+    double AngSize;             // Demie taille angulaire en Radian
+    double CosAngSize;          // Cosinus de la demi taille angulaire
+    double SurfaceTemperature;  // guess! Kelvin
+    double MyTheta;             // coordinates on the Sky Radian
+    double MyPhi;
+    UnitVector* pVSource;          // Direction de la source
+    
+ // Utilitaires de calcul
+ //	double h_kT;
+ //	double Prefac;		// 8*pi*h/c2 J m-2 s3
+ 	
+    virtual double spectre(double freq);
+    // Nous ne sommes plus dans la limite des petits carts de temperature
+    // Surcharge ncessaire
+    // GHz 3
+
+    virtual double powerDensAmpli(double theta,double phi) {
+  		// Return power density Amplidude at coordinates
+    	// W m-2 st-1 GHz-4, Vaut Prefac*DT(theta,phi)
+  		double PDens= Prefac*DeltaT(theta,phi);
+  		return PDens;
+  	}
+    double DeltaT(double theta,double phi);
+    // Retourne 1. si l on pointe sur la source, 0. sinon
+    
+    //    virtual double powerInterne(float theta, float phi, DetFilter const& Fil); 
+    // return power (Watt/m2/steradian) in direction theta,phi
+};
+#endif
Index: /tags/SigPred_V1/SigPredictor/sigcalchackingtools.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/sigcalchackingtools.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/sigcalchackingtools.cc	(revision 799)
@@ -0,0 +1,170 @@
+			// 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 calculs "<<"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 calculs "<<"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);
+
Index: /tags/SigPred_V1/SigPredictor/sigcalctools.cc
===================================================================
--- /tags/SigPred_V1/SigPredictor/sigcalctools.cc	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/sigcalctools.cc	(revision 799)
@@ -0,0 +1,516 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <math.h>
+#include <iostream>
+#include <iostream>
+#include <fstream>
+#ifdef __MWERKS__
+   #include "mwerksmath.h"
+   #include "unixmac.h"
+#endif
+#include "sigcalctools.h"
+#include "lightdipole.h"
+
+//_______________ ici toutes les frequences sont en Hz ___________________________
+
+static SigCalcTool* pSigToolcur;
+
+double SigCalGLFreqFunc1(double freq) {
+	double temp1=(pSigToolcur->pLSrc)->spectre(freq);
+	double temp2=(pSigToolcur->pLobe)->spectre(freq);
+	double temp3=(pSigToolcur->pFilter)->transmission(freq);
+	
+	return	temp1*temp2*temp3;
+}
+
+double SigCalGLFreqFunc2(double freq) 
+{ 
+	// Integration function for GLInteg
+  double temp1=
+	(pSigToolcur->pLSrc)->powSpecDens((pSigToolcur->VPointe).Theta(),(pSigToolcur->VPointe).Phi(),freq);
+  double temp2=(pSigToolcur->pLobe)->weigth(pSigToolcur->VCur,pSigToolcur->VPointe,pSigToolcur->VY,freq);
+  double temp3=(pSigToolcur->pFilter)->transmission(freq);
+  return temp1*temp2*temp3;
+}
+
+SigCalcTool::SigCalcTool(AbsLightSource* pLightSrc, AbsLobeNoPolar* pLobeNoPolar, 
+	SpectralResponse* pFilt):pLSrc(pLightSrc),pLobe(pLobeNoPolar),pFilter(pFilt)
+{	
+	SigCalcToolInit(); 
+}
+
+void SigCalcTool::SigCalcToolInit()
+{	emptySignal=false;
+// Compute frequency integration boundaries
+	cout<< "Initialisation Calctool"<<endl;
+	FreqMin=max(pLobe->minFreq(), pFilter->minFreq());
+	FreqMax=min(pLobe->maxFreq(), pFilter->maxFreq());
+	if(FreqMax<FreqMin) { 
+		emptySignal=true;
+		cerr<< "Frequency max is lower than Frequency Min in SigCalcTool"<<endl;
+		cerr<< "check consistency of lobes and Filters"<<endl;
+	}
+// Computation Options
+	if(pLSrc->IsMappedPowerSrc())
+	{  if(!pLobe->IsFreqSep())
+		{ cerr<<" Sigcalctool error: using a LightMapPowerInband with a lobe non freq separable"<<endl;
+		  cerr<<" Did you change lobe between constructing the map and running sigcalctool?"<<endl;
+		  cerr<<" Program exited"<<endl;
+		  exit(-1.);
+		}
+		Option=IsLightMapPowerInband;
+		pIntegrale= new GLInteg();
+			// Pour eviter un plantage dans ~SigCalcTool
+	}
+	
+	else if(pLSrc->IsFreqSep()&&pLobe->IsFreqSep()) {
+		Option=AllSeparable;
+		pIntegrale= new GLInteg(SigCalGLFreqFunc1,FreqMin,FreqMax); //en Hz.
+		pSigToolcur=this;	
+		pIntegrale->NStep(200);			// Integration tres srieuse
+		IntegSpectOverFreq=pIntegrale->Value();
+	}
+
+	else 
+	{	Option=NonSeparable;
+		pIntegrale= new GLInteg(SigCalGLFreqFunc2,FreqMin,FreqMax);
+	    pIntegrale->NStep(10); 		  // Pour aller plus vite. Serieux si le filtre est "compact"
+	}
+// Computation Resolution
+	RAngComp=pLSrc->LSrcResol();	  // On integre sur la resolution de la carte
+	if(RAngComp==0.) 
+	{	RAngComp=pLobe->lobeResol();
+		if(RAngComp==0.) 
+		{	cerr<<" Bizarre un lobe de resolution nulle?"<<endl;
+			RAngComp= 5.e-4; 	// Radians
+								// On prend la resolution nominale de Planck
+		}
+	}
+	if(RAngComp<pLobe->lobeResol())
+	{	cerr<<" SigCalcTool: LightSource resolution lower than expected lobe resolution"<<endl;
+		cerr<<" Not healthy: Ckeck consistency"<<endl;
+	} 
+	cout<<"Resolution de calcul: "<<RAngComp<<" Radian"<<endl<<endl;
+}
+
+double SigCalcTool::compPixel(double theta,double phi) {
+	UnitVector VP(theta,phi);
+	UnitVector VYbidon=VP.VperpPhi();	  
+// Compute unit vector perpendicular to Vpoin at same theta
+	return compPixel(VP,VYbidon);
+}
+
+double SigCalcTool::compPixel(UnitVector& VP, UnitVector& VdirectY){
+  double returnRes=0.;
+  VPointe=VP;
+  VY=VdirectY;
+  VX=VY^VP;
+  if(!emptySignal) returnRes=powerInteg(); // On integre sur la sphere 
+  return returnRes;
+}
+
+
+double SigCalcTool::calcPowerDens() const{
+// Compute the power integrated on frequency dependance, (Lobe and LightSource)
+  pSigToolcur=(SigCalcTool*) this;
+  double returnRes;
+  double poidlobe;
+  double Puiss;
+  switch (Option) 
+  {
+  	case AllSeparable:
+  	{ 
+  	  poidlobe=(pSigToolcur->pLobe)->weigthAmpl(VCur,VPointe,VY);  // ss dimensions
+/* 
+  	  if (poidlobe>.1) 
+ 	  { cout<<poidlobe<<endl;
+  	  }
+*/ 	   
+	 
+  	  Puiss=(pSigToolcur->pLSrc)->powerDensAmpli(VCur.Theta(),VCur.Phi());
+  	  		// W m-2 st-1 Hz-1 
+  	  returnRes=Puiss * poidlobe * IntegSpectOverFreq;	// W / m2 / st
+  	  return returnRes;
+  	}
+  	case IsLightMapPowerInband:
+  	{   
+  	 //	cout<<"VCur.Theta: "<<VCur.Theta()<<"VCur.Phi(): "<<VCur.Phi()<<endl;
+  		poidlobe= (pSigToolcur->pLobe)->weigthAmpl(VCur,VPointe,VY);
+  		Puiss= (pSigToolcur->pLSrc)->powerDensAmpli(VCur.Theta(),VCur.Phi());
+  		returnRes=Puiss * poidlobe;
+  		return returnRes;
+  	}
+  	
+	default: 
+	{  // Cas NonSeparable
+	   // Integration over at coordinates
+	  returnRes=pIntegrale->Value();
+	  return returnRes;
+  	}
+ 
+  }
+}
+
+
+#define NBStepCircleMin (12)
+
+double SigCalcTool::powerInteg() {
+  	// compute power on detector
+  
+  double powerInteg=0.; 
+  	// Sum of the incominig power on detector.
+  UnitVector VPoin;
+    // VPointe Boresigth du telescope microonde
+    // VPoin direction priviliegiee du lobe, autour de laquelle on calcule
+    // VCur, vecteur courant du calcul.
+//  double thetaCur, phiCur;    // Coordinates of VCur
+  							  // Units is radian
+  							  
+ 
+  
+  //------Initialisation of Lobe integration------------------------------------------
+  double angShift=0.;       // Angular distance from VPoin
+  double angShiftLimit;		// On calcule jusqu'a angShiftLimit de VPoin
+  
+  if(pLSrc->IsQPtSrc())
+  {	 double ang1=pLSrc->getAngSize()+pLobe->AngleMax();
+  	 VPoin=pLobe->VecShift(VPointe, VY); 
+     if (ang1>=M_PI) { } //rien 
+  	 else 
+  	 {  double cosinus=VPoin*pLSrc->GetVSrcCenter();
+  	    if (cosinus<cos(ang1)) return 0.;
+  		//C'est le cas ou la source est trop loin de la direction pointe
+  	 }
+  	 	// Maintenant on intgre
+  	 angShiftLimit=ang1;
+  }
+  else 
+  {
+     VPoin=pLobe->VecShift(VPointe, VY);      
+     angShiftLimit=pLobe->AngleMax();  
+  }
+  
+  // On va tourner autour de VPoin
+  // Compute unit vector perpendicular to Vpoin at same theta
+  UnitVector VPerp;
+  VPerp=VPoin.VperpPhi();
+
+  double dAngShift=AngResComp(0.)*1.1;
+    // AngleSteps are not necessarily constant.
+    // factor 1.1 to raise ambiguities in nearby pixel integration
+  double lastAngShiftMax;                            
+    // Needed to compute accurately solid angle values
+
+  VCur=VPoin;
+  
+  powerInteg+=calcPowerDens()*diffSolidAng(0.,dAngShift/2.);
+  lastAngShiftMax= dAngShift/2.;
+
+  long NbPasOneCircle;
+  long CircleNumber=0;   // no du cercle en cour: 
+  						 // Gestion des dcalages pour un echantillonnage en quinconce
+  double solidAngStepCircle; 
+  float stepAngCircle;
+
+  ///---------- Lobe integration-----------------------------------------
+  // generate vectors around VPoin at angular distance angShift.
+  // Compute power flux from foreground in this direction
+  // Weigth  with weigth function and solid angle
+  dAngShift=AngResComp(lastAngShiftMax);
+
+  while((lastAngShiftMax+dAngShift)<angShiftLimit){
+  	CircleNumber++;
+    angShift=lastAngShiftMax+dAngShift/2.;
+    
+    VCur=VPoin.Rotate(VPerp,angShift);
+    
+    // Compute number of step and associates on a circle  
+   	NbPasOneCircle=(long) (2*M_PI*sin(angShift)/sin(dAngShift));
+  	if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
+    stepAngCircle=2*M_PI/NbPasOneCircle;
+    solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShift+dAngShift/2.)/NbPasOneCircle;
+ 	// MRotAround=RotVec(VPoin,stepAngCircle);
+ 	
+    //----------- integrate on a circle -------------------
+    if((CircleNumber%2)==0) VCur=VCur.Rotate(VPoin,stepAngCircle/2.);
+    	// Pour un echantillonnage en quinconce
+    	
+    for(long i=0;i<NbPasOneCircle;i++)
+    {  
+      //cout<< "rotation numb: "<< i<<endl;	
+      powerInteg+=calcPowerDens()*solidAngStepCircle;
+      VCur=VCur.Rotate(VPoin,stepAngCircle);
+    }	// end of circle
+    
+    lastAngShiftMax+=dAngShift;
+    dAngShift=AngResComp(lastAngShiftMax);
+  }
+  
+  // On s'occupe des effets de bord: un dernier tour!
+  // On change le code pour eviter les instabilites dues a dAngShift tres petit
+  CircleNumber++;
+  angShift=(angShiftLimit+lastAngShiftMax)/2.;
+  
+  VCur=VPoin.Rotate(VPerp,angShift);
+  // Compute number of step and associates on a circle  
+  NbPasOneCircle=(long) 2*M_PI*sin(angShift)/sin(AngResComp(angShift));
+  if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
+  stepAngCircle=2*M_PI/NbPasOneCircle;
+  solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShiftLimit)/NbPasOneCircle;
+
+  //----------- integrate on last circle -------------------
+  for(long i=0;i<NbPasOneCircle;i++)
+  {
+     powerInteg+=calcPowerDens()*solidAngStepCircle;
+     VCur=VCur.Rotate(VPoin,stepAngCircle);
+  }
+  //end of last circle
+
+  //end of integration
+  
+// cout<<"On a termine un point, OUFF"<< endl;
+  return powerInteg;
+}
+
+double SigCalcTool::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.compPixel(M_PI/2.,M_PI);	// Un pixel au hasard
+ 	return tempeCNoir/puissNorm;				// Kelvin RaleighJeans/(Watt/m2)
+ 
+}
+ 
+double SigCalcTool::NormKelvinCMB()
+{	
+ 	double deltatempeCNoir=1.;					// Kelvin
+ 	LightNormTCMB DeltaCorpsNoir(deltatempeCNoir, RAngComp);
+ 	SigCalcTool ToolDeltaCMB(&DeltaCorpsNoir,pLobe,pFilter);
+ 	double puissNorm = ToolDeltaCMB.compPixel(M_PI/2.,M_PI);	// Un pixel au hasard
+ 	return deltatempeCNoir/puissNorm;			// KelvinCMB/(Watt/m2)
+ 
+}
+
+/*
+double SigCalcTool::CalcInBandPower(double theta, double phi)
+{
+  double returnRes=0.;
+  UnitVector VP(theta,phi);
+  UnitVector VYbidon=VP.VperpPhi();	  
+// Compute unit vector perpendicular to Vpoin at same theta
+  VCur=VP;
+  VPointe=VP;
+  VY=VYbidon;
+  VX=VY^VP;
+  if(!emptySignal) returnRes=calcPowerDens(); // On integre sur la frequence 
+  return returnRes;
+}
+*/
+
+double SigCalcTool::AngResComp(double angle) const
+{
+  double AngRes;
+  if(pLSrc->IsQPtSrc()) AngRes=RAngComp;
+  else AngRes=RAngComp*pLobe->ResolutionCurve(angle);
+  return AngRes;
+}
+
+
+
+double SigCalcTool::max(double a, double b) const{
+  if(a>b) return a;
+  else return b;
+}
+
+ double SigCalcTool::min(double a, double b) const{
+  if(a<b) return a;
+  else return b;
+}
+
+double SigCalcTool::CalcLobeSize(double frequency) 
+{
+  	// Compute lobe extension in steradians
+  	
+  if(frequency== -10.) frequency=(FreqMin+FreqMax)/2.;
+  
+  double SizeInteg=0.; 
+  	// Sum of the incominig power on detector.
+  UnitVector VPoin;
+    // VPointe Boresigth du telescope microonde
+    // VPoin direction priviliegiee du lobe, autour de laquelle on calcule
+    // VCur, vecteur courant du calcul.							 
+  
+  //------Initialisation of Lobe integration------------------------------------------
+  double angShift=0.;       			// Angular distance from VPoin
+  double angShiftLimit=pLobe->AngleMax();		// On calcule jusqu'a angShiftLimit de VPoin
+  
+  
+  // On va tourner autour de VPoin
+  // Compute unit vector perpendicular to Vpoin at same theta
+  UnitVector VPerp;
+  VPerp=VPoin.VperpPhi();
+
+  double dAngShift=AngResComp(0.)*1.1;
+    // AngleSteps are not necessarily constant.
+    // factor 1.1 to raise ambiguities in nearby pixel integration.
+  double lastAngShiftMax;                            
+    // Needed to compute accurately solid angle values
+  UnitVector VCur;
+  VCur=VPoin;
+  
+  SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
+  lastAngShiftMax= dAngShift/2.;
+
+  long NbPasOneCircle;
+  long CircleNumber=0;   // no du cercle en cour: 
+  						 // Gestion des dcalages pour un echantillonnage en quinconce
+  double solidAngStepCircle; 
+  float stepAngCircle;
+
+  ///---------- Lobe integration-----------------------------------------
+  // generate vectors around VPoin at angular distance angShift.
+  // Compute power flux from foreground in this direction
+  // Weigth  with weigth function and solid angle
+  dAngShift=AngResComp(lastAngShiftMax);
+
+  while((lastAngShiftMax+dAngShift)<angShiftLimit)
+  {
+    CircleNumber++;
+    angShift=lastAngShiftMax+dAngShift/2.;
+    
+    VCur=VPoin.Rotate(VPerp,angShift);
+    
+    // Compute number of step and associates on a circle  
+   	NbPasOneCircle=(long) (2*M_PI*sin(angShift)/sin(dAngShift));
+  	if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
+    stepAngCircle=2*M_PI/NbPasOneCircle;
+    solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShift+dAngShift/2.)/NbPasOneCircle;
+ 	
+    //----------- integrate on a circle -------------------
+    if((CircleNumber%2)==0) VCur=VCur.Rotate(VPoin,stepAngCircle/2.);
+    	// Pour un echantillonnage en quinconce
+    	
+    for(long i=0;i<NbPasOneCircle;i++)
+    {   	
+      SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
+      VCur=VCur.Rotate(VPoin,stepAngCircle);
+    }	// end of circle
+    
+    lastAngShiftMax+=dAngShift;
+    dAngShift=AngResComp(lastAngShiftMax);
+  }
+  
+  // On s'occupe des effets de bord: un dernier tour!
+  // On change le code pour eviter les instabilites dues a dAngShift tres petit
+  CircleNumber++;
+  angShift=(angShiftLimit+lastAngShiftMax)/2.;
+  
+  VCur=VPoin.Rotate(VPerp,angShift);
+  // Compute number of step and associates on a circle  
+  NbPasOneCircle=(long) 2*M_PI*sin(angShift)/sin(AngResComp(angShift));
+  if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
+  stepAngCircle=2*M_PI/NbPasOneCircle;
+  solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShiftLimit)/NbPasOneCircle;
+
+  //----------- integrate on last circle -------------------
+  for(long i=0;i<NbPasOneCircle;i++)
+  {
+     SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
+     VCur=VCur.Rotate(VPoin,stepAngCircle);
+  }
+  //end of last circle
+  
+  //end of integration  
+  
+  return SizeInteg;
+}
+ 
+double SigCalcTool::diffSolidAng(double ang1,double ang2) const
+{	double returnVal;		// Steradians
+ 	if(!pLSrc->IsPtSourceS()) returnVal= fabs(2*M_PI*(cos(ang1)-cos(ang2)));
+ 		// Cas d'une source tendue.
+ 	else  returnVal= 1.;
+ 		// Cas d'une source ponctuelle. Sont flux est exprime en W/m2/Hz
+ 		// Pas d'angle solide. 
+ 	return returnVal;
+ 		
+}
+
+// should be included as a class member, would template member function 
+// work on all compilers
+
+static AbsLobeNoPolar* AddInBandPowerpLobe;
+static AbsLightSource* AddInBandPowerpLSrc;
+static SpectralResponse* AddInBandPowerpFilter;
+static double AIBtheta;
+static double AIBphi;
+
+static double AddInBandPowerFreqFunc1(double freq)
+{  // Integration function for GLInteg
+  double temp1= AddInBandPowerpLSrc->powSpecDens(AIBtheta,AIBphi,freq);
+  double temp2= AddInBandPowerpLobe->spectre(freq);
+  double temp3= AddInBandPowerpFilter->transmission(freq);
+  return temp1*temp2*temp3;
+}
+
+template <class T> void addInInBandPowerMap(PixelMap<T>& Map, SigCalcTool& Tool)
+{	// No spatial integration on the lobe
+	// Valid if lobe is separable in frequency
+	// Test  
+  AddInBandPowerpLobe=Tool.getpLobe();
+  AddInBandPowerpLSrc=Tool.getpLSrc();
+  AddInBandPowerpFilter=Tool.getpFilter();
+  if(!AddInBandPowerpLobe->IsFreqSep())
+  {	cerr<<" Adding power to a map using a lobe non separable in frequency is inconsistent"<<endl;
+    cerr<<" No power added, addInBandPower skipped"<<endl;
+    return;
+  }
+    
+  long PixelNumber= Map.NbPixels();
+  double out;
+  T temp;
+  if(Tool.getOption()==AllSeparable)
+  {		// Fast !
+  	double FreqIntFactor=Tool.getIntegSpectOverFreq();
+  	for(long k=0; k<PixelNumber; k++)
+  	{  Map.PixThetaPhi(k,AIBtheta,AIBphi);
+	   out= AddInBandPowerpLSrc->powerDensAmpli(AIBtheta,AIBphi)*FreqIntFactor;
+    		 // Lobe weigth do no enters here
+	   temp= (T) out;
+       Map(k)+= temp;
+   // if((k%200)==0) cout<<"200 points calculs "<<"NbPoint Total= "<<k<<endl;
+  	}
+  
+  }
+  else
+  {
+    if(AddInBandPowerpLSrc->IsFreqSep()) 
+    { double FreqMax=Tool.getFreqMax();
+      double FreqMin=Tool.getFreqMin();
+      double out;
+      GLInteg Integrale(AddInBandPowerFreqFunc1,FreqMin,FreqMax);
+      Integrale.NStep(10);	// Serieux!
+  	  for(long k=0; k<PixelNumber; k++)
+  	  {  
+  	    Map.PixThetaPhi(k,AIBtheta,AIBphi);
+                		 // Lobe weigth do no enters here
+        out=Integrale.Value();
+    		 // Lobe weigth do no enters here
+	    temp= (T) out;
+        Map(k)+= temp;
+      }
+  	}
+  }
+   return;
+}
+
+template void addInInBandPowerMap(PixelMap<float>& Map, SigCalcTool& tool);
+template void addInInBandPowerMap(PixelMap<double>& Map, SigCalcTool& tool);
Index: /tags/SigPred_V1/SigPredictor/sigcalctools.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/sigcalctools.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/sigcalctools.h	(revision 799)
@@ -0,0 +1,136 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef SigCalcTools_H
+#define SigCalcTools_H
+
+#include "pixelmap.h" 
+#include "abslightsource.h"
+#include "abslobenopolar.h"
+#include "specresp.h"
+#include "integ.h"
+
+enum CalcOptions {NonSeparable,AllSeparable,IsLightMapPowerInband};
+
+class SigCalcTool {
+friend double SigCalGLFreqFunc1(double freq);
+friend double SigCalGLFreqFunc2(double freq);
+
+public:
+  SigCalcTool(AbsLightSource* pLSrc, AbsLobeNoPolar* pLobe, SpectralResponse* pFilter);
+  ~SigCalcTool() {	delete pIntegrale;}
+  
+ inline void SetLightScr(AbsLightSource* pLightScr) 
+ 	// require to rerun SigCalcToolInit
+ {	pLSrc=pLightScr;
+	delete pIntegrale;
+	SigCalcToolInit();
+ }
+ double compPixel(double theta, double phi);
+    // Return power on detector (Watt/m2) for this Lobe and filter
+    // Pointing at these ccordinates.
+    // Exact if lobes are symmetrical by rotation around pointed axe
+  
+ double compPixel(UnitVector& VP, UnitVector& VY);
+  	// Return power on detector (Watt/m2) for this Lobe and filter
+    // Pointing at these ccordinates. with this telescope orientation
+    // Exact
+ 
+ double CalcLobeSize(double frequency= -10.);
+ // Calcule l'extension spatiale du lobe de cet outil
+ // par integration numerique, calquee sur 
+ // L'integration spatiale de powerInteg();
+ // VALEUR RETOUR EN STERADIAN
+ // Doit EVOLUER SI powerInteg CHANGE;
+ 
+ double NormKelvinRayleighJeans();
+ 	// Compute te normalisation factor to go from Watt/m2 to KelvinRaleighJeans
+ 	// same computation option than the map options
+ 	// KelvinRJ/(Watt/m2)
+ 	
+ double NormKelvinCMB();
+ 	// Compute te normalisation factor to go from Watt/m2 to KelvinCMB
+ 	// same computation option than the map options
+ 	// KelvinCMB/(Watt/m2)
+ 	
+// double CalcInBandPower(double theta, double phi);
+ 	// return pwwer density (W/m2/st) integred over frequency response
+ 	// of LSource, Filter and Lobe
+ 	
+// Ces fonctions servent au hackingtools. elle palient le fait que je n'ai pas reussi
+// A declarer friend la fonction (refu du compilo de donner acces aux variables protected.
+// friend void addInInBandPowerMap(PixelMap<double>& Map, SigCalcTool& Tool);
+// Ne pas utiliser! 	
+ inline AbsLightSource* getpLSrc() const {return pLSrc;}
+ inline AbsLobeNoPolar* getpLobe() const {return pLobe;}
+ inline SpectralResponse* getpFilter() const {return pFilter;}
+ inline CalcOptions getOption() const {return Option;}
+ inline double getFreqMin() {return FreqMin;}
+ inline double getFreqMax() {return FreqMax;}
+ inline double getIntegSpectOverFreq() {return IntegSpectOverFreq;}
+ 
+protected:
+ AbsLightSource* pLSrc;
+ AbsLobeNoPolar* pLobe;
+ SpectralResponse* pFilter;
+ CalcOptions Option;	// Used to drives computations optimisations
+ GLInteg* pIntegrale;
+ 
+ bool emptySignal;
+ double RAngComp;		// angular resolution of computations
+ 
+ double FreqMax;		// Hz
+ double FreqMin;		// Hz
+ double IntegSpectOverFreq;
+
+ UnitVector VPointe;	// Boresigth du telescope (axe "Z" in Ticra focal plan)
+ UnitVector VY;  		// Axe "Y" of telescope ("horizontal") usefull when lobe are not symetrical
+ UnitVector VX;			// By rotation around Z axis.
+ UnitVector VCur;		// VCur, vecteur courant du calcul.
+ void SigCalcToolInit();
+ inline double AngResComp(double angle) const;
+    // Return current angle resolution of lobe computation 
+    // From map resolution lobe parameters. Unit: Radian 
+    
+ UnitVector VecOrigin(UnitVector VPointe);
+    // Define where high resolution computations starts. 
+    // Usually, it's where the horns point to
+    // Usefull with QuasiPtSources. You then want to compute around the source.
+    // Usefull when lobe as offsets in direction, 
+    	// Because he is not in the central part of focal plane. (TTicra data)
+ 
+ // Tools dealing with integration over space or frequency or both.  
+ double calcPowerDens() const;
+  	// compute the power integrated on frequency dependance, (Lobe and LightSource and filter)
+    // on pourra en avoir plusieur versions
+ 
+ double powerInteg();
+ 	// Compute the power recieved by a bolometer (Integration over space)
+ 	// Pointing towards XX
+ 	// From a source LSrc, with a lobe parametrised in Lobe, and an additionnal filter Filter
+	// WeigthVsFreq
+
+ // Computation tools       
+    // Math tools !!!!!!!-------------------------------------------------------
+ inline double max(double a, double b)const;
+ inline double min(double a, double b)const;
+            //compute solid angle between cones of aperture angles ang1 and ang2
+ 
+ inline double diffSolidAng(double ang1,double ang2) const;
+};
+
+// Tools for computing Maps
+template <class T> void addToSkyMap(PixelMap<T>& Map, SigCalcTool& Tool);
+template <class T> void compSkyMap(PixelMap<T>& Map, SigCalcTool& Tool);
+template <class T> void addInInBandPowerMap(PixelMap<T>& Map, SigCalcTool& tool);
+  	// Exact if lobes are symmetrical by rotation around pointed axe
+	// Return a map at the requested resolution and frequency
+
+// Tools for manipulating maps	
+template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2);
+template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2);
+template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2);
+template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map);
+template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min, 
+	double& Max, double& Moy, double& Var);
+
+#endif
Index: /tags/SigPred_V1/SigPredictor/tabulatedfilter.h
===================================================================
--- /tags/SigPred_V1/SigPredictor/tabulatedfilter.h	(revision 799)
+++ /tags/SigPred_V1/SigPredictor/tabulatedfilter.h	(revision 799)
@@ -0,0 +1,29 @@
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#ifndef TabulatedFilter_H
+#define TabulatedFilter_H
+
+// Filter class with unequally spaced (in freq) samples 
+#include <stdio.h>
+#include <stdlib.h>
+//#include <fstream)
+
+#include "specresp.h"
+#include "numrecipes.h"
+
+
+class TabulatedFilter:public SpectralResponse
+{	public:
+		TabulatedFilter(char filename[]);
+		virtual ~TabulatedFilter();	
+		virtual double transmission(double) const;
+		void TraceFiltre(ostream& f, long nbPoint=1000);
+	protected:
+		static NumRecipes NR;
+		unsigned long nbSample;
+		float *Transmi;
+		float *frequency;	
+};
+
+
+#endif
