// 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 "integ.h"
#include "lightdiffdust.h"
#include "fitsioserver.h"

LightDiffDust::LightDiffDust(int_4 nside=256) {
  nlat=nside;
  resolution=2*3.14159/4./nlat; // Resolution angulaire de la carte lue pour ce fond.
                      			// Radians
  sprintf(Name,"Poussi¸re Diffuse Schlegel");                  			
   
  // Mod¸le Physique                  			
  TDust=18.0;
  index=2.;
  double lref=240.e-6;		// Longueur d'onde carte schlegel
  fRef=3.e8/lref;			// Frequence de la carte Schelgel
							// Car nous parlons en frequence
  index2=-0.15;	
  Hfraction=0.05;
  double kbol=1.38e-23;
  double freeAtPole=5.3e-6;
  double freeAtPoleMJy_Sr=freeAtPole*2.*kbol
  	* pow(53.e9/3.e8,2) * pow(fRef/53.e9,index2) / 1.e-20;
  double dustAtPoleMJy_Sr=3*0.9/1.1;
  CorrelatedfreeWeight=freeAtPoleMJy_Sr/dustAtPoleMJy_Sr;
//  nuPlanckfRef=nuPlanck(fRef);
  
  cout<<"Objet LightDiffDust :"<<endl;
  cout <<"Resolution de travail nlat= "<<nlat<<endl;
  
  try { pMap=new SphereGorski<float> (nlat);	}
  catch (bad_alloc) {
    cerr<<"memory booking error Lightdiffdust constructor nlat= "<<nlat<<endl;
    exit(-1); 
  }
  cerr<<"Vous avez rˇservˇ: "<<12*nlat*nlat*sizeof(float)<<" Octets"<<endl;
  char filename[150];
  
#ifndef __MWERKS__  
  char *PATHDataLScr;
  PATHDataLScr=getenv("PATHDataLScr");
#endif

  sprintf(filename,"%sdust_res%04i.fits",PATHDataLScr,nlat);
	
  cout<<"On lit la carte "<<filename<<endl;
  if(nlat<1024) cerr<<"la carte avec nlat=1024 est dispo"<<endl;
  
  FitsIoServer FitsServer;
  FitsServer.load((*pMap),filename);
		// unite: MJy/st. Dixit Romain
}

LightDiffDust::~LightDiffDust(){
  delete pMap;
}