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

// ArchTOIPipe           (C)     CEA/DAPNIA/SPP IN2P3/LAL
//                               Eric Aubourg
//                               Christophe Magneville
//                               Reza Ansari
// $Id: toi2map.cc,v 1.12 2001-11-10 00:14:55 aubourg Exp $

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

////////////////////////////////////////////////////////////////////////
TOI2Map::TOI2Map(SphereHEALPix<r_8>* sph,SphereHEALPix<r_8>* wsph)
  : mSph(sph), mWSph(wsph), mWSphInternal(false)
{
 SetEquinox();
 SetCoorIn();
 SetCoorOut();
 SetTestFlag();
 SetTestMin();
 SetTestMax();

 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::Print(ostream & os)
{
  os<<"TOI2Map::Print -- Sphere NLat = "<<mSph->SizeIndex()<<endl

    <<"   - Equinoxe="<<mActualYear<<endl

    <<"   - TypCoorIn: Gal("<<((mTypCoorIn&TypCoordGal)?1:0)
    <<") Eq("<<((mTypCoorIn&TypCoordEq)?1:0)
    <<") Deg("<<((mTypCoorIn&TypCoordDD)?1:0)
    <<") Hour("<<((mTypCoorIn&TypCoordHD)?1:0)
    <<") Rad("<<((mTypCoorIn&TypCoordRR)?1:0)<<")"<<endl

    <<"   - TypCoorOut: Gal("<<((mTypCoorOut&TypCoordGal)?1:0)
    <<") Eq("<<((mTypCoorOut&TypCoordEq)?1:0)
    <<") Deg("<<((mTypCoorOut&TypCoordDD)?1:0)
    <<") Hour("<<((mTypCoorOut&TypCoordHD)?1:0)
    <<") Rad("<<((mTypCoorOut&TypCoordRR)?1:0)<<")"<<endl

    <<"   - Tests: Flag("<<mTFlag<<") bad="<<mBadFlag
    <<" / Value Min("<<mTMin<<")="<<mValMin
    <<" , Max("<<mTMax<<")="<<mValMax<<endl;
}

////////////////////////////////////////////////////////////////////////
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!");
}
if( !(mTypCoorIn&TypCoordEq || mTypCoorIn&TypCoordGal) ) {
  cout<<"TOI2Map::run() - Input Coordinates not Eq or Gal! "<<endl;
  throw ParmError("TOI2Map::run() - Input Coordinates not Eq or Gal!");
}
if( !(mTypCoorOut&TypCoordEq || mTypCoorOut&TypCoordGal) ) {
  cout<<"TOI2Map::run() - Output Coordinates not Eq or Gal! "<<endl;
  throw ParmError("TOI2Map::run() - Output Coordinates not Eq or Gal!");
}

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

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

 int bufsz = 100;
 uint_8* bfgbolo = new uint_8[bufsz];
 double* bbolo   = new double[bufsz];
 double* bc1     = new double[bufsz];
 double* bc2     = new double[bufsz];
 int i0 = -1;

// Remplissage des spheres
for(int s=snb;s<=sne;s++) {
  if (i0<0 ||
      s<i0 ||
      s>=i0+bufsz) {
    i0 = s;
    int nget = bufsz;
    if (sne-s+1 < bufsz) {
      nget = sne-s+1;
    }
    getData(0, i0, nget, bc1);
    getData(1, i0, nget, bc2);
    getData(2, i0, nget, bbolo, bfgbolo);
  }

  
  uint_8 fgbolo = 0;
  double bolo;
  //              Equatoriales   /   Galactiques
  // coord1,2 =   alpha,delta    /   gLon,gLat
  //double coord1 = getData(0,s);
  //double coord2 = getData(1,s);

  //getData(2,s,bolo,fgbolo);
  bolo = bbolo[s-i0];
  fgbolo = bfgbolo[s-i0];
  double coord1 = bc1[s-i0];
  double coord2 = bc2[s-i0];
  if (s%100 == 0) wontNeedBefore(s-1);
  if(mTFlag && fgbolo&mBadFlag) continue;
  if(mTMin && bolo<mValMin) continue;
  if(mTMax && bolo>mValMax) continue;

  // sphere phi   entre [0,2*Pi] en radians
  // sphere theta entre [0,Pi]   en radians
  double phi=-1.,theta;
  CoordConvertToStd(mTypCoorIn,coord1,coord2);
  if(mTypCoorIn&TypCoordEq && mTypCoorOut&TypCoordGal) { // Eq -> Gal
    EqtoGal(mjd,coord1,coord2,&coord1,&coord2);
    phi   = coord1 * M_PI/180.;
  } else if(mTypCoorIn&TypCoordGal && mTypCoorOut&TypCoordEq) { // Gal -> Eq
    GaltoEq(mjd,coord1,coord2,&coord1,&coord2);
    phi   = coord1 * M_PI/12.;
  } else if(mTypCoorOut&TypCoordGal) { // Gal -> Gal
    phi   = coord1 * M_PI/180.;
  } else if(mTypCoorOut&TypCoordEq) { // Eq -> Eq
    phi   = coord1 * M_PI/12.;
  }
  theta = (90.-coord2) * M_PI/180.;
  if(phi<0.   || phi>=2*M_PI) continue;
  if(theta<0. || theta>=M_PI) continue;

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

 cout<<"TOI2Map::run(): Fin de boucle sur les sampleNum"<<endl;

// 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()<<"%"
     <<"  NFill["<<NFILL<<"] ="<<endl;
 for(ii=0;ii<NFILL;ii++) {cout<<NFill[ii]<<" "; if(ii%10==9) cout<<endl;}
 cout<<endl;

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

return;                                                                            
}