source: Sophya/trunk/SigPredictor/sigcalchackingtools.cc@ 808

                        // 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 calculŽs "<<"NbPoint Total= "<<k<<endl;
  }
   return;
}

template <class T> void compSkyMap(PixelMap<T>& Map, SigCalcTool& Tool)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
  T temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
        temp= (T) Tool.compPixel(theta,phi);
//    if (temp<1.e-35) temp=0.;
                                // BUGG XXXXXX
                        // pour vivre avec les pb de fichiers fits D.Y.
        //if (out!=0.) cout<< out <<'\t'<< temp<<endl;           
    Map(k)= temp;
    if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
  }
   return;
}

template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
   // temp=(T1) Map2.PixValSph(theta,phi);
    if(Map2.ContainsSph(theta,phi)) Map(k) += (T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
    Map.PixThetaPhi(k,theta,phi);
 //   temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
    if(Map2.ContainsSph(theta,phi)) Map(k) = Map(k)-(T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2)
{
  double theta, phi;
  long PixelNumber= Map.NbPixels();
  //T1 temp;
  for(long k=0; k<PixelNumber; k++)
  { 
     Map.PixThetaPhi(k,theta,phi);
 //   temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
     if(Map2.ContainsSph(theta,phi)) 
         if (Map2.PixValSph(theta,phi)!=0.) 
                 Map(k) = Map(k)/(T1) Map2.PixValSph(theta,phi);
  }
   return;
}

template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map)
{
        long PixelNumber= Map.NbPixels();
        //T temp;
        for(long k=0; k<PixelNumber; k++) 
        {       
//              temp= (T) ( scalefactor*Map2(k) );
                Map(k) = (T) ( scalefactor*Map(k) );    
        }
}


template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min, 
        double& Max, double& Moy, double& sigma)
{       
        long PixelNumber= Map.NbPixels();
        double val=0.;
        double variance=0;
        Moy=0.;
        Min=1.e36;
        Max=-1.e-36;
        sigma=0.;
        
        for(long k=0; k<PixelNumber; k++) 
        {       val=Map(k); // cout<<k<<'\t'<<val<<'\n';
            if(val>Max) Max=val;
                if(val<Min) Min=val;
                Moy+=val;       
        }
        
        Moy/=PixelNumber;
        for(long k=0; k<PixelNumber; k++) variance+=(val-Moy)*(val-Moy);
        variance/=PixelNumber;
        sigma=sqrt(variance);
        
        if(Min==Max) 
        {       cerr<< "bizarre, une carte vide ou uniforme"<<endl;
                return 1;
        }
        return 0;
}

// On instancie les plus courants
template void addToSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
template void addToSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);
template void compSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
template void compSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);      

// Idem pour les outils de cartes

template void scaleMap(double scalefactor, PixelMap<double>& Map2);
template void scaleMap(double scalefactor, PixelMap<float>& Map2);

template int MinMaxSigMap(PixelMap<float>& Map, double& Min, 
        double& Max, double& Moy, double& sigma);
template int MinMaxSigMap(PixelMap<double>& Map, double& Min, 
        double& Max, double& Moy, double& sigma);


template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
template  void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);

template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);

template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);

template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
template  void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);