Index: /trunk/SigPredictor/abscalctool.cc
===================================================================
--- /trunk/SigPredictor/abscalctool.cc	(revision 1149)
+++ /trunk/SigPredictor/abscalctool.cc	(revision 1149)
@@ -0,0 +1,295 @@
+
+
+			// Dominique YVON, CEA/DAPNIA/SPP 02/2000
+
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#ifdef __MWERKS__
+    
+   #include "unixmac.h"
+#endif
+
+//#include "integ.h"
+
+#include "abscalctool.h"
+
+#include "lightdipole.h"
+#include "sigcalctools.h"
+
+
+double AbsCalcTool::compPixelQD(double theta,double phi) 
+{	
+//	cout<<"theta: "<<theta<<"\tphi: "<<phi<<endl;
+	UnitVector VP(theta,phi);
+	UnitVector VYbidon=VP.VperpPhi();	  
+// Compute unit vector perpendicular to Vpoin at same theta
+	return compPixel(VP,VYbidon);
+}
+
+double AbsCalcTool::NormKelvinRayleighJeans()
+{	
+ 	double tempeCNoir=10000.;	
+ 					// Kelvin
+ 	double CutFreq=1.380662e-23*tempeCNoir/6.626176e-34/5.;
+ 	if(FreqMax>1.380662e-23*tempeCNoir/6.626176e-34/5.) 
+ 	{	cerr<< "RaleighJeans approximation is not valid for this frequency"<<endl;
+ 		cerr<< "Frequency: "<< FreqMax<<" in SigCalcTool::NormRayleighJeans"<<endl;
+ 	}
+ 	
+ 	LightBlackBody CorpsNoir(tempeCNoir, RAngComp);
+ 	SigCalcTool ToolRJ(&CorpsNoir,pLobe,pFilter);
+ 	double puissNorm = ToolRJ.compPixelQD(M_PI/2.,M_PI);	// Un pixel au hasard
+ 	return tempeCNoir/puissNorm;				// Kelvin RaleighJeans/(Watt/m2)
+ 
+}
+ 
+double AbsCalcTool::NormKelvinCMB()
+{	
+ 	double deltatempeCNoir=1.;					// Kelvin
+ 	LightNormTCMB DeltaCorpsNoir(deltatempeCNoir, RAngComp);
+ 	SigCalcTool ToolDeltaCMB(&DeltaCorpsNoir,pLobe,pFilter);
+ 	double puissNorm = ToolDeltaCMB.compPixelQD(M_PI/2.,M_PI);	// Un pixel au hasard
+ 	return deltatempeCNoir/puissNorm;			// KelvinCMB/(Watt/m2)
+}
+
+double AbsCalcTool::diffSolidAng(double ang1,double ang2) const
+{	return fabs(2*M_PI*(cos(ang1)-cos(ang2)));		// Steradians
+ 		// Cas d'une source tendue.
+}
+
+double AbsCalcTool::max(double a, double b) const{
+  if(a>b) return a;
+  else return b;
+}
+
+double AbsCalcTool::min(double a, double b) const{
+  if(a<b) return a;
+  else return b;
+}
+
+int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo) {
+  // All coordinates are in Radian.
+  /* kmg_euler.c
+   *
+   * Converts between different coordinate systems. 
+   *
+   * SELECT   From           To         |   SELECT    From          To
+   *    1     RA-Dec (2000)  Galactic   |     4       Ecliptic      RA-Dec    
+   *    2     Galactic       RA-DEC     |     5       Ecliptic
+Galactic  
+   *    3     RA-Dec         Ecliptic   |     6       Galactic
+Ecliptic  
+   *
+   *    Date     Programmer Remarks
+   * ----------- ---------- -------
+   * 08-Aug-1999 K. Ganga   First version. Copied and modified EULER from 
+   *                        the IDL Astrolib. 
+   * May 2000, D. Yvon		Change coordinates units to Radians
+   */
+
+  /* Local Declarations */
+  double       a, b, sb, cb, cbsa;
+  long         i;
+  const double TWOPI   = 2.0*M_PI;
+  const double FOURPI  = 4.0*M_PI;
+  const double DEG2RAD = 180.0/M_PI;
+
+  /* J2000 coordinate conversions are based on the following constants
+   * eps    =  23.4392911111  Obliquity of the ecliptic
+   * alphaG = 192.85948       Right Ascension of Galactic North Pole
+   * deltaG =  27.12825       Declination of Galactic North Pole
+   * lomega =  32.93192       Galactic longitude of celestial equator
+   * alphaE = 180.02322       Ecliptic longitude of Galactic North Pole
+   * deltaE =  29.811438523   Ecliptic latitude of Galactic North Pole
+   * Eomega =   6.3839743     Galactic longitude of ecliptic equator
+   */
+
+  const double psi[6]    = {0.57477043300, 4.9368292465 ,
+                            0.00000000000, 0.0000000000 , 
+                            0.11142137093, 4.71279419371};
+  const double stheta[6] = {0.88998808748,-0.88998808748,
+                            0.39777715593,-0.39777715593,
+                            0.86766622025,-0.86766622025};    
+  const double ctheta[6] = {0.45598377618, 0.45598377618, 
+                            0.91748206207, 0.91748206207,
+                            0.49714719172, 0.49714719172};
+  const double phi[6]    = {4.9368292465 , 0.57477043300,
+                            0.0000000000 , 0.00000000000, 
+                            4.71279419371, 0.11142137093};
+
+  i    = select - 1;
+  a    = ai - phi[i];
+  b    = bi;
+  sb   = sin(b);
+  cb   = cos(b);
+  cbsa = cb*sin(a);
+  b    = -stheta[i]*cbsa + ctheta[i]*sb;
+  b    = ( (b >  1.0) ?  1.0 : b );
+  b    = ( (b < -1.0) ? -1.0 : b );
+  *bo  = asin(b);
+
+  a    = atan2( ctheta[i] * cbsa + stheta[i] * sb, cb * cos(a) );
+  *ao  = fmod(a + psi[i] + FOURPI, TWOPI);
+
+  /* Later */
+  return 0;
+}
+
+template <class T> void addToSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool) 
+{
+  double theta, phi;
+  long PixelNumber= Map.NbPixels();
+  cout<<"Nbre de pIxel a calculer dans addToSkyMap: "<<PixelNumber<<endl;
+  T temp;
+  for(long k=0; k<PixelNumber; k++)
+  { 
+    Map.PixThetaPhi(k,theta,phi);
+    temp=(T) Tool.compPixelQD(theta,phi); 
+//    if (temp<1.e-35) temp=0.;		// BUGG XXXXXX
+    								// pour vivre avec les pb de fichiers fits D.Y. 
+     Map(k)+=temp;
+     if((k%500)==0) cout<<"500 points calculs "<<"NbPoint Total= "<<k<<endl;
+  }
+   return;
+}
+
+template <class T> void compSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool)
+{
+  double theta, phi;
+  long PixelNumber= Map.NbPixels();
+  cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
+  T temp;
+  for(long k=0; k<PixelNumber; k++)
+  { 
+    Map.PixThetaPhi(k,theta,phi);
+	temp= (T) Tool.compPixelQD(theta,phi);
+//    if (temp<1.e-35) temp=0.;
+				// BUGG XXXXXX
+     			// pour vivre avec les pb de fichiers fits D.Y.
+	//if (out!=0.) cout<< out <<'\t'<< temp<<endl;		 
+    Map(k)= temp;
+    if((k%500)==0) cout<<"500 points 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, AbsCalcTool& Tool);
+template void addToSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
+template void compSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
+template void compSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);	
+
+// Idem pour les outils de cartes
+
+template void scaleMap(double scalefactor, PixelMap<double>& Map2);
+template void scaleMap(double scalefactor, PixelMap<float>& Map2);
+
+template int MinMaxSigMap(PixelMap<float>& Map, double& Min, 
+	double& Max, double& Moy, double& sigma);
+template int MinMaxSigMap(PixelMap<double>& Map, double& Min, 
+	double& Max, double& Moy, double& sigma);
+
+
+template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
+template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
+template  void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);
+
+template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
+template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
+template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);
+
+template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
+template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
+template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);
+
+template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
+template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
+template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);
+
Index: /trunk/SigPredictor/abscalctool.h
===================================================================
--- /trunk/SigPredictor/abscalctool.h	(revision 1149)
+++ /trunk/SigPredictor/abscalctool.h	(revision 1149)
@@ -0,0 +1,103 @@
+#ifndef AbsCalTool_H
+#define AbsCalTool_H
+
+#include <ostream>
+#include "unitvector.h"
+#include "pixelmap.h" 
+#include "specresp.h"
+#include "abslobenopolar.h"
+
+class AbsCalcTool
+{
+ public:
+ 	AbsCalcTool(double CompResol=0.,double freqMin=0.,double freqMax=-1.) 
+ 		:FreqMax(freqMax), FreqMin(freqMin), RAngComp(CompResol) {	}
+ 	//	Dummy values. Have to be initialised in offspring classes
+ 	
+ 	~AbsCalcTool() {	}
+ 	
+ 	virtual double compPixel(UnitVector& VP, UnitVector& VY)=0;
+ 	 // Return power on detector (Watt/m2) for this Lobe and filter
+    // Pointing at these ccordinates and orientation. Exact
+    
+   virtual void print(ostream& ReadMe)=0;
+   //	Main infos about calcTool
+    
+   virtual double CalcLobeSize(double frequency= -10.)=0;
+	// Calcule l'extension spatiale du lobe de cet outil
+	// par integration numerique, calquee sur 
+	// L'integration spatiale ddu signal physique;
+	// VALEUR RETOUR EN STERADIAN
+	// Doit EVOLUER SI l'integration spatiale du signal physique CHANGE;
+	
+ 	double compPixelQD(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
+   
+  // Normalisation tools
+	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)
+ 	
+ protected:
+ 	double RAngComp;		// angular resolution of computations
+ 	double FreqMax;		// Hz
+ 	double FreqMin;
+ 
+ 	AbsLobeNoPolar* pLobe;
+ 	SpectralResponse* pFilter;		// Hz
+ 	
+ 	 // 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
+ 
+	virtual inline double diffSolidAng(double ang1,double ang2) const;
+
+};
+ 	
+// Tools for computing Maps
+template <class T> void addToSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool);
+template <class T> void compSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool);
+template <class T> void addInInBandPowerMap(PixelMap<T>& Map, AbsCalcTool& 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);
+	
+// Tool for changing ccordinates
+// All coordinates are in Radian.
+  /* kmg_euler.c
+   *
+   * Converts between different coordinate systems. 
+   *
+   * SELECT   From           To         |   SELECT    From          To
+   *    1     RA-Dec (2000)  Galactic   |     4       Ecliptic      RA-Dec    
+   *    2     Galactic       RA-DEC     |     5       Ecliptic
+Galactic  
+   *    3     RA-Dec         Ecliptic   |     6       Galactic
+Ecliptic  
+   *
+   *    Date     Programmer Remarks
+   * ----------- ---------- -------
+   * 08-Aug-1999 K. Ganga   First version. Copied and modified EULER from 
+   *                        the IDL Astrolib. 
+   * May 2000, D. Yvon		 Change coordinates units to Radians
+   */
+int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo);
+
+#endif
