// 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 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(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 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.)*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 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::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 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;
}
 
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 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);