source: Sophya/trunk/SigPredictor/sigcalctools.cc@ 3264

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

#include <math.h>
#include <iostream>
#include <iostream>
#include <fstream>
#ifdef __MWERKS__
    
   #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 sŽrieuse
                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(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 pointŽe
         }
                // Maintenant on intgre
         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.);
    // AngleSteps are not necessarily constant.
  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 dŽcalages 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::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::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.);
    // AngleSteps are not necessarily constant.
  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 dŽcalages 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;
}

// 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 calculŽs "<<"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);