source: Sophya/trunk/ArchTOIPipe/ProcWSophya/toi2map.cc@ 1501

#include "machdefs.h"
#include "toimanager.h"
#include "pexceptions.h"
#include "ctimer.h"
#include "toi2map.h"
#include "xastropack.h"

////////////////////////////////////////////////////////////////////////
TOI2Map::TOI2Map(SphereHEALPix<r_8>* sph,SphereHEALPix<r_8>* wsph)
  : mSph(sph), mWSph(wsph), mWSphInternal(false), mTypCoor(false), fTypCoor(false), mActualYear(2001.)
{
 if(mSph->NbPixels()<1) {
  cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mSph "
      <<mSph->NbPixels()<<endl;
  throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
 }
 mSph->SetPixels(0.);
 int nlat = mSph->SizeIndex();

 if(mWSph==NULL) {
   mWSph = new SphereHEALPix<r_8>(nlat);
   mWSphInternal = true;
 } else {
   mWSphInternal = false;
   if(nlat != mWSph->SizeIndex()) mWSph->Resize(nlat);
 }
 if(mWSph->NbPixels()<1) {
   cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mWSph "
       <<mWSph->NbPixels()<<endl;
   throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
 }
 mWSph->SetPixels(0);

}

TOI2Map::~TOI2Map()
{
 if(mWSph && !mWSphInternal) delete mWSph;
}

////////////////////////////////////////////////////////////////////////
void TOI2Map::init() {
  cout << "TOI2Map::init" << endl;
  declareInput("Coord1In");     // input index 0
  declareInput("Coord2In");     // input index 1
  declareInput("BoloIn");      // input index 2
}

////////////////////////////////////////////////////////////////////////
void TOI2Map::run()
{
long snb = getMinIn();
long sne = getMaxIn();

if(snb>sne) {
  cout<<"TOI2Map::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
  throw ParmError("TOI2Map::run() - Bad sample interval");
}
if(!checkInputTOIIndex(0) || !checkInputTOIIndex(1) || !checkInputTOIIndex(2)) {
  cout<<"TOI2Map::run() - Input TOI (Coord1In or Coord2In or BoloIn) not connected! "<<endl;
  throw ParmError("TOI2Map::run() Output TOI (Coord1In or Coord2In or BoloIn) not connected!");
}

//---------------------------------------------------------
#define NFILL 25
try {

int ii;
uint_4 mNSnFill=0, mNpixFill=0, NFill[NFILL];
for(ii=0;ii<NFILL;ii++) NFill[ii]=0;

// Remplissage des spheres
for(int s=snb;s<=sne;s++) {
  int_8 fgbolo = 0; 
  double bolo;

  double coord1 = getData(0,s); // gLat ou delta entre [-90,90] en degres
  double coord2 = getData(1,s); // gLon entre [0,360[ en degres ou alpha entre [0,24[ en heures
  getData(2,s,bolo,fgbolo);

  if(coord2<-90. || coord2>90.) fgbolo=1; 
  if((coord1<0.) || (!mTypCoor && coord1>=24.) || (mTypCoor && coord1>=360.) ) {fgbolo=1; 
  cout << "!!!!!!!!" <<coord1 << endl;
  }
  
  
  
  if(bolo<-32767.) fgbolo=1;   // Bidouille Archeops
  
  if(!fgbolo) {
    // sphere phi entre [0,2*Pi] en radians
    // sphere phi entre [0,2*Pi] en radians
    // sphere theta entre [0,Pi] en radians
    double phi,theta;
    if(fTypCoor && !mTypCoor) { //on a alpha,delta et on veut l,b
      double mjd = MJDfrYear(mActualYear);
      EqtoGal(mjd,coord1,coord2,&coord1,&coord2);
    } 

    if(!fTypCoor && mTypCoor) { //on a l,b et on veut alpha,delta
      double mjd = MJDfrYear(mActualYear);
      GaltoEq(mjd,coord1,coord2,&coord1,&coord2);
    }
    
    if(fTypCoor) phi = coord1 * M_PI/180.;
    else  phi = coord1 * M_PI/12.;

    theta = (90.-coord2)*M_PI/180.;

    int_4 ipix = mSph->PixIndexSph(theta,phi);
    (*mSph)(ipix) += bolo;    
    ((*mWSph)(ipix))++;
    mNSnFill++;
  }
}

// Remplissage des spheres
 for(int_4 i=0;i<mSph->NbPixels();i++) {
   r_8 wf = (*mWSph)(i);
   if( wf > 0. ) {
     mNpixFill++;
     (*mSph)(i) /= wf;
   }
   int_4 nf = int_4(wf);
   if(nf>=NFILL) nf=NFILL-1; NFill[nf]++;
 }

 cout<<"TOI2Map::run(): mNpixTot="<<mSph->NbPixels()
     <<"  mNpixFill="<<mNpixFill
     <<"  mNSnFill="<<mNSnFill<<endl
     <<" --> FracSky="<<mNpixFill*100./(double)mSph->NbPixels()<<"%"<<endl
     <<"NFill["<<NFILL<<"] = "<<endl;
 for(ii=0;ii<NFILL;ii++) cout<<NFill[ii]<<" ";
 cout<<endl;

//---------------------------------------------------------
} catch (PException & exc) {
  cout<<"TOI2Map: Catched Exception "<<(string)typeid(exc).name()
      <<"\n .... Msg= "<<exc.Msg()<<endl;
}

return;                                                                            
}