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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef __MWERKS__
   #include "mwerksmath.h"
   #include "unixmac.h"
   #include "macenvvariables.h"
#endif
#include "nbrandom.h"
#include "tmatrix.h"
#include "numrecipes.h"
#include "integ.h"
#include "fitsioserver.h" 

#include "lightgalaxnoresol.h"


LightGalaxNoResol::LightGalaxNoResol(int_4 nside=512)
{
 
 /* cerr<< "LightGalaxNoResol en travaux, ne pas utiliser"<<endl;
    exit (-1);
 */
  sprintf(Name,"Galaxies non resolues");
  nlat=nside;              
  if(nside!=2048) 
  {  cerr<<" Warning : Donnees Guiderdonni compatible avec nlat=2048"<<endl;
  	 cerr<<" Vous avez choisi nlat= "<<nlat<<endl;
  }
  nbFreq=31;			   // Nbre de frequences calculŽes 
  NbPixUnCote=441;
  nbPixel=NbPixUnCote*NbPixUnCote;    // Nbre de pixel cartes Guiderdonni
  nbPixelKept=65535;       // Ce que lon peut stocker dans un U2
  pNR= new NumRecipes;

  // Tableau des frŽquences auxquelles sont calculŽes les cartes guiderdonni
  dataFreq=new r_8[nbFreq]; 
  dataFreqDegueux= new float[nbFreq];
  // Calcul des frŽquences des cartes guiderdonni
  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];
  }
  // Chaque pixel du ciel est caracterise par un spectre ŽchantillonnŽ
  // ˆ nbfreq frequences extrait des cartes Guiderdonni.
  // On stocke le no de pixel de la carte guiderdonni correspondant au pixel ciel
  
  try {pSphereIndex=new SphereGorski<uint_2>(nlat); }
  catch(bad_alloc) {
  	cerr<<"Memory booking error in LightGalaxNoResol: pSphereIndex"<<endl;
  	exit(-1);
  }
    // Sphre faisant la correspondance ciel, no de pixel de carte guiderdonni
  nbPixelLight=pSphereIndex->NbPixels();
  	// On attribue a chaque pixel ciel un pixel carte Guiderdonni.
  for( long i=0; i< nbPixelLight; i++) 
  { 
    // double test= MyRan()*nbPixelKept;
    // bogue sur Mac frand01() retourne toujours la meme valeur.
    // On utilise ran() standard.
    
    (*pSphereIndex)(i)= (int_2) (MyRan()*nbPixelKept);
  }
  
    // Stockage des donnŽes Guiderdonni petites cartes 
        // RŽserve l'espace mŽmoire nŽcessaire au stockage des cartes partielles
  ppskyData= new r_4* [nbPixelKept];	
  		// Un tableau pour chaque point du ciel
  try {
  	for(int i=0; i<nbPixelKept; i++) ppskyData[i]= new r_4 [nbFreq];
  }
  catch(bad_alloc) {
  	cerr<<"Memory booking error in LightGalaxNoResol: ppskyData"<<endl;
  	exit(-1);
  }

  
/*// ca ne marche plus. se renseigner de comment faire!.
  if(pSphereIndex==NULL){
    cerr<<"Memory booking error LightGalaxNoResol constructor"<<endl;
    exit(-1);
  }  
*/

  int_4 size=sizeof(r_4)*(nbFreq*nbPixelKept) + sizeof(uint_2)*nbPixelLight;
  cout<<"Objet LightGalaxNoResol: "<<endl;
  cout<<"Vous avez reservŽ: "<<size<<" Octets de mŽmoire vive"<<endl;

  // Lit les cartes Guiderdonni:les charge dans dataTampon,puis dans skyData
  #ifndef __MWERKS__
    char* PATHDataLScr=getenv("PATHDataLScr");
  #endif
  
  FitsIoServer FitsServer;
  
  TMatrix<float> UneCarte(NbPixUnCote, NbPixUnCote);
  
  long cartefreqMHz; 
  char filename[150];
  for(int noFreq=0; noFreq<nbFreq; noFreq++)
  {
  	  // GŽnre le nom des cartes
	cartefreqMHz=(long)(dataFreq[noFreq]/1000000.+0.5);
    sprintf(filename,"%seran102_441_%08imhz.fits",PATHDataLScr,cartefreqMHz);

	FitsServer.load(UneCarte,filename);
//	cout<<"Carte Lue: "<<filename<<'\t'<<"FrŽquence: "<<dataFreq[noFreq]<<endl;
	
	long indexUneCarte1=0;
	long indexUneCarte2=0;
	
    // On remplit LightData 
    for(int NoPixel=0;NoPixel<nbPixelKept;NoPixel++)
    { 
      indexUneCarte1= (long) (NoPixel/NbPixUnCote); 		// Division entiere
      indexUneCarte2= NoPixel%NbPixUnCote;
      ppskyData[NoPixel][noFreq]= UneCarte(indexUneCarte1,indexUneCarte2);
      cout <<"NoPixel: "<<NoPixel<<" Valeur: ";
      cout << ppskyData[NoPixel][noFreq]<<'\t'<<UneCarte(indexUneCarte1,indexUneCarte2)<<endl;
    }
    
  }  // C'est dans la boite LightData
  
  // C'est fini pour le constructeur
}

LightGalaxNoResol::~LightGalaxNoResol(){
 for(int i=0; i<nbPixelKept; i++) delete[] ppskyData[i];
 delete[] ppskyData;
 delete pNR;
 delete[] dataFreq;
 delete[] dataFreqDegueux;
 delete[] pSphereIndex; 

}

double LightGalaxNoResol::powSpecDens(double theta, double phi, double freq){
  long index=pSphereIndex->PixIndexSph(theta,phi);
  int_2 pixNb=pSphereIndex->PixVal(index);

  float* pbinSpectre= ppskyData[pixNb]; //10-26 W/m2/Hz/st
 
  // On interpole au troisieme ordre.
  float InterpRes=0.;
  float InterpResErr=0.;
  pNR->polint(dataFreqDegueux-1,pbinSpectre-1,3,freq,&InterpRes,&InterpResErr);

  return (double)InterpRes*1.e-26;  // W/m2/Hz/st
}


static double randMax=RAND_MAX;

double LightGalaxNoResol::MyRan() 
{
	return rand()/randMax;
}