source: Sophya/trunk/SigPredictor/lightgalaxresol.cc@ 800

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

#include <stdio.h>
#include <stdlib.h>

#include <math.h>
#include <iostream>
#include <fstream>

#ifdef __MWERKS__
   #include "mwerksmath.h"
   #include "unixmac.h"
   #include "macenvvariables.h"
#endif

#include "nbrandom.h"
#include "lightgalaxresol.h"
#include "numrecipes.h"
#include "integ.h"


LightGalaxResol::LightGalaxResol(int_4 nside) {
  sprintf(Name, "Galaxies Resolues");
  nlat=nside;
  resolution=2*M_PI/4./nlat;            // Radians
  nbFreq=31;
  nbPointSource=2711;
  PtSourceS=true;
  
  pNR=new NumRecipes;
  LastFreqIndex=1;
  
// Calcul des frŽquences des cartes guiderdonni
  dataFreq=new  r_8[nbFreq];
  dataFreqDegueux= new float[nbFreq];

  double fmin=15.e9;
  double fmax=15.e12;
  double dLogFreq=log10(fmax/fmin)/(nbFreq-1);

  for(int i=0; i<nbFreq; i++) 
  {
    dataFreq[i]=fmin*pow(10,i*dLogFreq);
    dataFreqDegueux[i]=(float)dataFreq[i];                              // Pour la NR
  }

// Booking of skymap of index
  try { pLightMap= new SphereGorski<uint_2> (nlat); }
  catch(bad_alloc){
    cerr<<"Erreur Reservation de memoire classe LightGalacResol, nlat= "<<nlat<<endl;
    exit(-1);
  }
  
  nbPixelLight=pLightMap->NbPixels();
  int size=sizeof(uint_2)*nbPixelLight;
  cout<<"Objet LightGalaxResol: "<<endl;
  cout<<"Vous avez reservŽ: "<<size<<"Octets de mŽmoire vive"<<endl;
  
  
// La sphere est initialisŽe ˆ nbPointSource+1, cad pas de sources
  int_2 bidon=nbPointSource+1;
  for(int i=0; i<nbPixelLight; i++)  (*pLightMap)(i)=bidon;
  
// On distribue les sources uniformement sur le ciel
  int_4 skyIndex;
  for(int_2 i=0; i<nbPointSource; i++) 
  {
    skyIndex=(int)(MyRan()*nbPixelLight);
    if ( (*pLightMap)(skyIndex)==(nbPointSource+1)) (*pLightMap)(skyIndex)=i;
    else 
    { 
      while(! ((*pLightMap)(skyIndex) == (nbPointSource+1)) ) 
      { 
//      float test;        
//      test= MyRan();          // frand01() ne fonctionne pas au 05/12/99 sur Mac
        skyIndex=(int)(MyRan()*nbPixelLight);
      }
           (*pLightMap)(skyIndex)=i;
    } 
  }


// _____________ On lit les parametres des sources _______________________________________

        // Reservation de l'espace memoire de stockage
  try { 
        ppPointSourceData= new float* [nbPointSource];
        for(int i=0; i<nbPointSource; i++) ppPointSourceData[i]= new float[nbFreq];
  }
  catch(bad_alloc) {
        cerr<<" Memory booking error in LightGalaxResol: ppskyData "<<endl;
        exit(-1);
  }

        // Lecture des fichiers
  char filename[150];
  double DumStarNb;
  char DumChaine[20];
  r_4 thisdata;


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

  for(int NoFreq=0; NoFreq<nbFreq; NoFreq++) 
  {      
     // cout<<"NoFreq :"<<NoFreq<<"Frequence: "<<dataFreq[NoFreq]<<endl;
     // On genere le bon nom de fichier
         sprintf(filename,"%seran102_441_%08imhz.s_rare",PATHDataLScr,(int)(dataFreq[NoFreq]/1.e6+.5));
 
         ifstream InFile(filename,ios::in);             //      Fichier Ouvert!
         
#ifdef __MWERKS__
         if (!InFile.is_open())         // cxx DEC est trop con pour comprendre is_open()
#else
        if (!InFile)
#endif
         {      cerr<<"erreur a l'ouverture du fichier de donnees: "<<filename<<endl;
                exit(-1);
         }


         InFile>>DumStarNb; 
         if(DumStarNb!=nbPointSource) {
                cerr<<"ya une bogue dans le nb de galaxptsources."<<endl;
                cerr<<"fichier: "<< filename<<"nb sources annoncees:"<< DumStarNb<<endl;
         }
         InFile>>DumChaine;
         
         for(long noPtSource=0; noPtSource<nbPointSource; noPtSource++) 
         {
                InFile>> thisdata;
                ppPointSourceData [noPtSource] [NoFreq] =thisdata;
//              cout<<noPtSource<<'\t'<<ppPointSourceData [noPtSource] [NoFreq]<<endl;
         }
         InFile.close();
  }
/*     for (int NoFreq=0; NoFreq<nbFreq; NoFreq ++) 
           cout<< NoFreq <<'\t'<< ppPointSourceData [] [NoFreq]<< endl;
*/
}


LightGalaxResol::~LightGalaxResol() 
{
  // LibŽrons la mŽmoire!! Bient™t mai!
  for(int i=0; i<nbPointSource; i++) delete[] ppPointSourceData[i];  
  delete[] ppPointSourceData;
  delete pLightMap;
  delete[] dataFreqDegueux;
  delete[] dataFreq;
  delete pNR;
}

double LightGalaxResol::powSpecDens(double theta, double phi, double freq)
         //(Watt/m2/Hz)
{     
  long pixelNo= pLightMap->PixIndexSph(theta, phi);
  int_2 NoSource=(*pLightMap)(pixelNo);
 
  if ( NoSource==(nbPointSource+1)) return 0.;  // Pas de source
  else 
  {  
    float* binSpectre=(ppPointSourceData[NoSource]); 
    float InterpRes;
    float InterpResErr;
    
    float thisFreq=freq;

        //Ou est la frequence dans le tableau?
        pNR->hunt(dataFreqDegueux-1,(unsigned long) nbFreq, thisFreq, &LastFreqIndex);
        if((LastFreqIndex==0)||(LastFreqIndex==nbFreq)) 
        {       cerr<<"Appel a powSpecDens dans LightGalaxResol hors limites de frequence"<<endl;
                LastFreqIndex=1;
                return 0;
        }
                        
        
      // On interpole 3 eme ordre. Num Reciepies fonction.
    pNR->polint(dataFreqDegueux-1+LastFreqIndex,binSpectre-1+LastFreqIndex,
        2,thisFreq,&InterpRes,&InterpResErr);
      // InterpRes en Jansky:10-26 W/m2/Hz           
    return (double) (InterpRes*10.e-26);        // W/m2/Hz
      // RMQ: Le Jansky est une unite de source ponctuelle !
      // Pas de steradian dans les homogeneites.
 }
}

static double randMax=RAND_MAX;
double LightGalaxResol::MyRan() {
        return rand()/randMax;
}