source: Sophya/trunk/SigPredictor/makebolotimeline.cc@ 3395

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

#include <iomanip.h>
#include <new>

# define ERROR_LABEL -32768

#include "makebolotimeline.h"
//#include "fitsioserver.h"


MakeBoloTimeline::MakeBoloTimeline(char file[], AbsLightSource& LSrc, 
        AbsLobeNoPolar& Lobe, SpectralResponse& Filter, Normalisation Norm)
{       IsSigCalctool=true;
        IsPtSrcInBandtool=false;
        pSigCalcTool=new SigCalcTool(&LSrc,&Lobe,&Filter);
        MakeBoloTimeline(file,*pSigCalcTool,Norm);
}

MakeBoloTimeline::MakeBoloTimeline(char file[], LightPtSrcLevSInBand& LSrc, 
        MeanFreqLobe& Lobe, LevSPanckBand band, Normalisation Norm)
{       IsSigCalctool=false;
        IsPtSrcInBandtool=true;
        
        pPtSrcInBandCalTools=new PtSrcInBandCalTools(&LSrc,&Lobe,band);
        MakeBoloTimeline(file,*pPtSrcInBandCalTools,Norm);
}

MakeBoloTimeline::MakeBoloTimeline(char file[], AbsCalcTool& Tool, Normalisation Norm)
        :NormMethode(Norm), ptool(&Tool)
{       

#ifndef __MWERKS__
  char* PATHResults=getenv("PATHResults"); 
#endif

        
        NormRJ=ptool->NormKelvinRayleighJeans();
        NormTCMB=ptool->NormKelvinCMB();
        saveReadMe(file);
        
        char filecur[150];
        sprintf(filecur,"%s.dat",file);
        
        pMystr= new ofstream(filecur,ios::out|ios::trunc);
//      pMystr->setf(ios::scientific);
        (*pMystr)<<setprecision(5);
        
#ifdef VisuIsActive
        // Visu de controle
        try
        {       pSphere=new SphereHEALPix<r_4>(64);
                pSphereInt=new SphereHEALPix<uint_2>(64);
        }
        catch (bad_alloc)
        {       cerr<<"bad alloc in MakeBoloTimeline:"<<file<<endl;
        }

        sprintf(Name, "%s%s",file,".visu");
#endif
}


MakeBoloTimeline::~MakeBoloTimeline()
{
        (*pMystr)<<flush;
        pMystr->close();
        delete pMystr;

        if(IsSigCalctool) delete pSigCalcTool;
        if(IsPtSrcInBandtool) delete pPtSrcInBandCalTools;

#ifdef VisuIsActive
        for(long index=0; index<(*pSphere).NbPixels();index++)
                if ((*pSphereInt).PixVal(index)!=0) 
                        (*pSphere).PixVal(index)=(*pSphere).PixVal(index)/(*pSphereInt).PixVal(index);
        
//      FitsIoServer FitsServer;
        #ifndef __MWERKS__
                 char* PATHResults=getenv("PATHResults"); 
        #endif

//      FitsServer.Mollweide_picture_projection((*pSphere),Name);  // BUGG XXXXXXXXX
        
        delete pSphere;
        delete pSphereInt;
#endif
}

void MakeBoloTimeline::addToStreamArchTOI(double theta,double phi)
{       //       Theta en heures, Phi en degres, systeme RA Dec
        if((theta==ERROR_LABEL)||(phi==ERROR_LABEL)) return;
        
        // On passe en Radian
        double thetaRad=theta*15./180*M_PI;
        double phiRad=phi/180.*M_PI;
        
        // Go to Galactic ccordinates
        double lgalax=0.;               
        double bgalax=0.;
        kmg_eulerRad(thetaRad,phiRad,1,&lgalax,&bgalax);
        
        // Add to time line
        addToStream(lgalax,bgalax);
}

void MakeBoloTimeline::addToStream(double lgalax,double bgalax)
{       // lgalax et bgalax en Radian, Systeme coord Galactique J2000

        double theta=M_PI/2.-bgalax;    // Coordonnees Spheriques galactiques
        double phi=lgalax;      
        double value;
        
        switch (NormMethode)
        {
                case RayleighJeans:     // Value en TempeRJ
                {       value= NormRJ*ptool->compPixelQD(theta,phi);
                        break;
                }
                case TempeCMB:  // Value en tempeCMB
                {       value= NormTCMB*ptool->compPixelQD(theta,phi);
                        break;
                }
                default:                // Value en power
                {       value   =ptool->compPixelQD(theta,phi);
                        break;
                }
        }
                (*pMystr)<<lgalax<<'\t'<<bgalax<<'\t'<<value<<'\n';
                
                long PixelNumb=(*pSphere).PixIndexSph(theta,phi);
                (*pSphere).PixVal(PixelNumb)+=value;
                (*pSphereInt).PixVal(PixelNumb)++;
                
        return;
}

void  MakeBoloTimeline::saveReadMe(char File[])
{       char filecur[150];
        // Fill Header file for ascii output
        sprintf(filecur,"%s.ReadMe",File);
        cout<<filecur<<endl;
        ofstream ReadMe(filecur,ios::out|ios::trunc);
        
        ReadMe<<" Fichier de timeline pour:"<<endl;
        ptool->print(ReadMe);

        
        ReadMe<< setprecision(3);
        ReadMe<<" Constantes de normalisation"<<endl;
        ReadMe<<" RayleighJeans: "<<NormRJ<< " KelvinRJ/(W/m2)"<<endl;
        ReadMe<<" TempeCMB: "<<NormTCMB<< " Kelvin CMB/(W/m2)"<<endl;
        
        ReadMe<<"Cette carte est en unitŽ: ";
        switch(NormMethode)
        {       
                case TempeCMB: 
                {       ReadMe<< "Kelvin CMB"<<endl;
                        break;
                }
                case RayleighJeans:
                {       ReadMe<< "Kelvin RayleighJeans"<<endl;
                        break;
                }
                default: ReadMe<< "Watt/m2 effectif de mirroir"<<endl;
        }
        
        ReadMe<<"Fichier coordonŽes galactique J2000ecrit selon le format"<<endl;
        ReadMe<<"l_galax(Radian)\tb_galax(Radian)\tAbove Unit"<<endl;
        
        // J'espere que c'est assez 
        ReadMe.close();
        
        return;
}