source: Sophya/trunk/SigPredictor/ptsrcinbandcalctools.cc@ 3937

#include <iostream>
#include "ptsrcinbandcalctools.h"
#include "lightptsrclevsinband.h"
#include "spherehealpix.h"

PtSrcInBandCalTools::PtSrcInBandCalTools(LightPtSrcLevSInBand* pLightSource,MeanFreqLobe* pLob, LevSPanckBand thisBand)
        :pLightSrc(pLightSource), Band(thisBand)
{       pLobe=pLob;     
        FreqMax=pLobe->maxFreq();
        FreqMin=pLobe->minFreq();
        pFilter=pLob->pRespShape;
        RAngComp =  pLobe->lobeResol();
}

double PtSrcInBandCalTools::compPixel(UnitVector& VP, UnitVector& VY)
{       double OuvAng=pLobe->AngleMax();
        double thetaPointe=VP.Theta();
        double PhiPointe=VP.Phi();

        double ThetaMin=thetaPointe-OuvAng; if(ThetaMin<0.) ThetaMin=0.;
        double ThetaMax=thetaPointe+OuvAng;     if(ThetaMax>M_PI) ThetaMax=M_PI;
        
/*
        cout<< "test Passage"<<endl;
        cout<<VP.Theta()<<'\t'<<VP.Phi()<<'\t'<<VY.Theta()<<'\t'<<VY.Phi()<<endl;
*/      
        double deltaPhi;
        {       // Compute limits in Phi
                double deltaPhi1;
                if(ThetaMin!=0.) deltaPhi1=OuvAng/sin(ThetaMin);
                else    deltaPhi1=M_PI;
                double deltaPhi2;
                if(ThetaMax!=0.) deltaPhi2=OuvAng/sin(ThetaMax);
                else deltaPhi2=M_PI;
                
                // Compute deltaPhi
                if(deltaPhi1>deltaPhi2) deltaPhi=deltaPhi1;
                else deltaPhi=deltaPhi2;
        }
        double PhiMin=PhiPointe-deltaPhi;
        double PhiMax=PhiPointe+deltaPhi;

/*      // Useless!
        if((PhiMin<-M_PI)||(PhiMax>M_PI))
        {       PhiMin=-M_PI;
                PhiMax=M_PI;
        }
*/
        double PowerTot=0.;
        
        const UnitVector VecBidon(ThetaMin,0.);
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterlow=pLightSrc->MapOfPtSrc.lower_bound(VecBidon);
        const UnitVector VecBidon2(ThetaMax,0.);
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterhigh=pLightSrc->MapOfPtSrc.upper_bound(VecBidon2);

        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
        double PhiSrc;
                
//      long compteur=0;
//      long compteur2=0;

        for(iter=iterlow; iter!=iterhigh; iter++)
        {       PhiSrc=(*iter).first.Phi();
                if( (PhiSrc>PhiMin) && (PhiSrc<PhiMax) )
                {       double poid=pLobe->weigthAmpl((*iter).first,VP,VY);
                        PowerTot+=(*iter).second.getPower(Band)*poid;
//                      compteur++;
//                      if (poid!=0.) compteur2++;
                }
        }
        return PowerTot;
}

double PtSrcInBandCalTools::CalcLobeSize(double frequency)
        // I KNOW! Awfully heavy. But It's the cost of generality.
        // This function should be called only once!
{       double AngMax=pLobe->AngleMax();

        // On cherche la resolution healpix correspondante
        double x= 0.5*M_PI/RAngComp;
        long nlat=1;
        while (x>1.)
        {x/=2.; nlat*=2;}
        
        // On reserve une sphere pour en connaitres les coordonnŽes
        SphereHEALPix<uint_2> dummSphere(nlat);
        
        double SolidAngPixel = dummSphere.PixSolAngle();
        long Nbtranch = dummSphere.NbThetaSlices()*AngMax/M_PI;
        double thetaStep=M_PI/Nbtranch+1.;
        
        TVector<uint_2> VecValues;
        TVector<int_4> VecIndice;
        r_8 dumTheta;
        r_8 dumPhi;
        
        double AngSolid=0.;
        UnitVector VZ(0.,0.,1.);
        
        for(double ThetaAnn=0.5*thetaStep;ThetaAnn<AngMax; ThetaAnn+=thetaStep)
        {       // Obtient l'indice d'un pixel sur l'anneau theta=ThetaAnn
                int_4 index=dummSphere.PixIndexSph(ThetaAnn,0.);
                dummSphere.GetThetaSlice(index,dumTheta,dumPhi,VecIndice,VecValues);
                
                double theta,phi;
                // Somme la valeur du lobe sur les pointe d'un anneau
                for(long index=0; index<VecIndice.Size(); index++) 
                {       int_4 ind=VecIndice(index);
                        dummSphere.PixThetaPhi(ind,theta,phi);
                        UnitVector VInteg(theta, phi);
                        UnitVector VY=VZ^VInteg;
                        
                        AngSolid+= pLobe->weigthAmpl(VInteg,VZ,VY)*SolidAngPixel;
                }
        }
        cout<<AngSolid<<endl;
        
        return AngSolid;
}

/*
double PtSrcInBandCalTools::powerInteg(double ThetaMin, double ThetaMax, 
                                double PhiMin, double PhiMax, UnitVector& VP, UnitVector& VY, bool InvLog)
{       
        const UnitVector VecBidon(ThetaMin,0.);
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterlow=pLightSrc->MapOfPtSrc.lower_bound(VecBidon);
        const UnitVector VecBidon2(ThetaMax,0.);
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterhigh=pLightSrc->MapOfPtSrc.upper_bound(VecBidon2);

        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
        double Power=0.;
        double PhiSrc;
        bool test;
                
//      long compteur=0;
//      long compteur2=0;

        for(iter=iterlow; iter!=iterhigh; iter++)
        {       PhiSrc=(*iter).first.Phi();
                test=(PhiSrc>PhiMin) && (PhiSrc<PhiMax);
                if (InvLog) test= (!test);      
                if(test)
                {       double poid=pLobe->weigthAmpl((*iter).first,VP,VY);
                        Power+=(*iter).second.getPower(Band)*poid;
//                      compteur++;
//                      if (poid!=0.) compteur2++;
                }
        }
        return Power;
}
*/