source: Sophya/trunk/SigPredictor/lightptsrclevsinband.cc@ 3516

#include <fstream>
#include "unitvector.h"
#include "lightptsrclevsinband.h"


LightPtSrcLevSInBand::LightPtSrcLevSInBand(char FileOfPtSrc[])
{       // le fichier de pointsource ordonnŽ en Theta
        sprintf(Name, "Sources Ponctuelles du Level S");
        readLevSPtSrc(FileOfPtSrc);
}

void LightPtSrcLevSInBand::readLevSPtSrc(char PtSrcfile[])
{       char Line[200];
        
        double theta, phi; 
        float f30, f44, f70, f100, f143, f217, f352, f545, f857;
        LevSPtSrcData UneSource;
        UnitVector UnCoord;
        
        ifstream HFI(PtSrcfile);
        if (!HFI)
        {       cout<< "Erreur a l'ouverture du fichier de donnees :"<<PtSrcfile<< endl;
                exit(-1);
        }
        
        long compteur=0;
        long compteur2=0;
        while (HFI)
        {       // Data in file in radian and Jansky
        //      cout<<compteur<<'\t'<<HFI.fail()<<'\n'<<endl;;
                HFI>>theta;
                HFI>>phi;
                HFI>>f30;
                HFI>>f44;
                HFI>>f70;
                HFI>>f100;
                HFI>>f143;
                HFI>>f217;
                HFI>>f352;
                HFI>>f545;
                HFI>>f857;      
                
                compteur++;
                if(HFI.fail())
                {       HFI.clear();
                        HFI.getline(Line, 199); 
                        if(HFI.eof())           // Fin de lecture du fichier
                        {       cerr<< "Fin de lecture du fichier: "<<PtSrcfile<<endl;
                                cerr<< "Point source number : "<<MapOfPtSrc.size()<<endl;
                                HFI.close();
                                break;
                        }
                        else cerr<<"erreur lecture ligne: "<<Line<<endl;
                }
                else
                {       // cout<<theta<<'\t'<<phi<<'\t'<<f30<<'\t'<<f44<<'\t'<<f70<<'\t'<<f100;
                        // cout<<'\t'<<f143<<'\t'<<f217<<'\t'<<f352<<'\t'<<f545<<'\t'<<f857<<endl<<endl;
                        
                        // In Memory Power is sored in Watt, angles in Radian
                        UnitVector UnCoord(theta, phi);
                        UneSource.flux30G=f30*1.e-26;
                        UneSource.flux44G=f44*1.e-26;
                        UneSource.flux70G=f70*1.e-26;
                        UneSource.flux100G=f100*1.e-26;
                        UneSource.flux143G=f143*1.e-26;
                        UneSource.flux217G=f217*1.e-26;
                        UneSource.flux352G=f352*1.e-26;
                        UneSource.flux545G=f545*1.e-26;
                        UneSource.flux857G=f857*1.e-26;
                        
                        const pair<const UnitVector,LevSPtSrcData>UneValeur=make_pair(UnCoord, UneSource);
                        
                        MapOfPtSrc.insert(MapOfPtSrc.end(),UneValeur);
                        compteur2++;
                }
        }
        cerr<< compteur<<'\t'<<compteur2<<endl;
        return;
}

#include "fitsspherehealpix.h"

void LightPtSrcLevSInBand::DumpMap(LevSPanckBand band)
{       // Write a Mollweide projection of the Point Sources Sky
        // power in the LevSPanckBand band
        SphereHEALPix <r_4> Sphere(128);

        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iter;
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterBegin=MapOfPtSrc.begin();
        multimap <UnitVector,LevSPtSrcData,LevSPtSrcCmpUnitVec>::iterator iterEnd=MapOfPtSrc.end();
        
        long indexSphere;
        double theta;
        double phi;
        double power;
        
        for(iter=iterBegin; !(iter==iterEnd); iter++)
        {       theta=(*iter).first.Theta();
                phi=(*iter).first.Phi();
                power=(*iter).second.getPower(band);
                
                indexSphere=Sphere.PixIndexSph(theta,phi);
                Sphere(indexSphere) += power;
        }
        
        FITS_SphereHEALPix <r_4> FSphere(Sphere);
        FSphere.Mollweide_picture_projection("SkyOfPtSrcInBandMollweide.fits");

        return;
}