// 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;
}