source: Sophya/trunk/SophyaLib/Samba/localmap.cc@ 413

#include "localmap.h"
#include <iostream.h>
#include "nbmath.h"
//*****************************************************************************
//++ 
// Class        LocalMap
//
// include      localmap.h nbmath.h
//
//    A local map of a region of the sky, in cartesian coordinates.
//    It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
//    (x0, y0 might be outside of this local map)
//    default value of (x0, y0) is middle of the map, center of 
//    pixel(nx/2, ny/2)
//
//    la carte est consideree comme un tableau a deux indices i et j, i etant 
//    indice de colonne et j indice de ligne. La carte est supposee resider 
//    dans un plan tangent, dont le point de tangence est repere (x0,y0) dans 
//    la carte et (theta0, phi0) sur la sphere celeste. L'extension de la 
//    carte est definie par les valeurs de deux angles couverts respectivement
//    par la totalite des pixels en x de la carte et la totalite des pixels 
//    en y. (SetSize()).
//    On considere un "plan de reference" : plan tangent a la sphere celeste 
//    aux angles theta=Pi/2 et phi=0. Dans ce plan L'origine des coordonnees 
//    est le point de tangence. L'axe Ox est la tangente parallele a 
//    l'equateur, dirige vers les phi croissants, l'axe Oy est parallele 
//    au meridien, dirige vers le pole nord.
//    De maniere interne a la classe une carte est definie dans ce plan de 
//    reference et transportee  jusqu'au point (theta0, phi0) de sorte que 
//    les axes restent paralleles aux meridiens et paralleles. L'utilisateur 
//    peut definir sa carte selon un repere en rotation par rapport au repere 
//    de reference (par l'angle entre le parallele et l'axe Ox souhaite -- 
//    methode SetOrigin(...))

//
//-- 
//++
//
// Links        Parents
//
//    PixelMap
//    
//-- 
//++
//
// Links        Descendants
//
//     
//-- 
//++
// Titre        Constructeurs
//--
//++

LocalMap::LocalMap()

//--

{InitNul();}

//++

LocalMap::LocalMap(int_4 nx, int_4 ny) : mSzX_(nx), mSzY_(ny)

//--
{
  InitNul();
  mNPix_=nx*ny;
  mPix_ = new r_8[mNPix_];
}


//++
// Titre        Destructeur
//--
//++
LocalMap::~LocalMap()

//--
{
  Clear();
}



void        LocalMap::WriteSelf(POutPersist& s) const
{
  char strg[256];
  if (mInfo_) {sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d HasInfo",  mSzX_, mSzY_);} else  {
    sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d ",  mSzX_, mSzY_);
  }         
  s.PutLine(strg);
  if (mInfo_)  s << (*mInfo_);
  s.PutI4(mSzX_);
  s.PutI4(mSzY_);
  s.PutI4(mNPix_);
  s.PutI4(x0_);
  s.PutI4(y0_);
  s.PutI4(originFlag_);
  s.PutI4(SzFlag_);
  s.PutR4(cos_angle_);
  s.PutR4(sin_angle_);
  s.PutR8(theta0_);
  s.PutR8( phi0_);
  s.PutR8( tgAngleX_);
  s.PutR8(tgAngleY_);
  s.PutR8s(mPix_, mNPix_);
}
void        LocalMap::ReadSelf(PInPersist& s)
{
  Clear();
  char strg[256];
  s.GetLine(strg, 255);
// Pour savoir s'il y avait un DVList Info associe
  bool hadinfo = false;
  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
  if (hadinfo) {    // Lecture eventuelle du DVList Info
    if (mInfo_ == NULL)  mInfo_ = new DVList;
    s >> (*mInfo_);
  }
  s.GetI4(mSzX_);
  s.GetI4(mSzY_);
  s.GetI4(mNPix_);
  s.GetI4(x0_);
  s.GetI4(y0_);
  s.GetI4(originFlag_);
  s.GetI4(SzFlag_);
  s.GetR4(cos_angle_);
  s.GetR4(sin_angle_);
  s.GetR8(theta0_);
  s.GetR8( phi0_);
  s.GetR8( tgAngleX_);
  s.GetR8(tgAngleY_);
  mPix_ = new r_8[mNPix_+1];
  s.GetR8s(mPix_, mNPix_);
}

//++
// Titre        Methodes
//--


//++
int_4 LocalMap::NbPixels() const

//    Retourne le nombre de pixels du découpage
//--
{
  return(mNPix_);
}

//++
r_8& LocalMap::PixVal(int_4 k)

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >= mNPix_) ) {
    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
    //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
  THROW(rangeCheckErr);

    //throw "LocalMap::PIxVal Pixel index out of range";
  }
  //return (*(pixel_->Data()+k));
  return(mPix_[k]);
}

//++

r_8 const& LocalMap::PixVal(int_4 k) const

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >= mNPix_) ) {
    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));

    throw "LocalMap::PIxVal Pixel index out of range"; 
  } 
  //return *(pixel_->Data()+k);
  return(mPix_[k]);
}

//++
int_4 LocalMap::PixIndexSph(float theta, float phi) const

//    Retourne l'indice du pixel vers lequel pointe une direction définie par 
//    ses coordonnées sphériques
//--
{
  int i,j;
  if (!(originFlag_==1) || !(SzFlag_==1) ) {
    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    exit(0);
  }
  // theta et phi en coordonnees relatives (on se ramene a une situation 
  // par rapport au plan de reference)
  float theta_aux=theta;
  float phi_aux=phi;
  UserToReference(theta_aux, phi_aux);
  // coordonnees dans le plan local en unites de pixels
  float x,y;
  AngleProjToPix(theta_aux, phi_aux, x,y);

  float xmin=-x0_-0.5;
  float xmax=xmin+mSzX_;
  if( (x >  xmax) || (x < xmin ) ) return(-1);
  float xcurrent=xmin;
  for( i=0; i < mSzX_; i++ ) {
    xcurrent += 1.;
    if( x < xcurrent ) break;
  }
  float ymin=-y0_-0.5;
  float ymax=ymin+mSzY_;
  if( (y >  ymax) || (y < ymin ) ) return(-1);
  float ycurrent=ymin;
  for( j=0; j < mSzY_; j++ ) {
    ycurrent += 1.;
    if( y < ycurrent ) break;
  }


  return (j*mSzX_+i);
}

//++
void  LocalMap::PixThetaPhi(int_4 k, float& theta, float& phi) const

//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
//--
{
  if (!(originFlag_==1) || !(SzFlag_==1) ) {
    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    exit(0);
  }
  int i,j;
  Getij(k,i,j);
  float X=float(i-x0_);
  float Y=float(j-y0_);
  // situation de ce pixel dans le plan de reference
  float x=X*cos_angle_-Y*sin_angle_;
  float y=X*sin_angle_+Y* cos_angle_;
  // projection sur la sphere
  PixProjToAngle(x, y,theta, phi);
  // retour au plan utilisateur
  ReferenceToUser(theta, phi);
}
 
//++
r_8    LocalMap::PixSolAngle(int_4 k) const
//    Pixel Solid angle  (steradians)
//    All the pixels have not necessarly the same size in (theta, phi)
//    because of the projection scheme which is not yet fixed. 
//    
//--
{
  int i,j;
  Getij(k,i,j);
  float X=float(i-x0_);
  float Y=float(j-y0_);
  float XR=X+float(i)*0.5;
  float XL=X-float(i)*0.5;
  float YU=Y+float(j)*0.5;
  float YL=Y-float(j)*0.5;
  // situation  dans le plan de reference
  float x0=XL*cos_angle_-YL*sin_angle_;
  float y0=XL*sin_angle_+YL* cos_angle_;
  float xa=XR*cos_angle_-YL*sin_angle_;
  float ya=XR*sin_angle_+YL* cos_angle_;
  float xb=XL*cos_angle_-YU*sin_angle_;
  float yb=XL*sin_angle_+YU* cos_angle_;
  // projection sur la sphere
  float tet0,phi0,teta,phia,tetb,phib;
  PixProjToAngle(x0, y0,tet0, phi0);
  PixProjToAngle(xa, ya,teta, phia);
  PixProjToAngle(xb, yb,tetb, phib);
  // angle solide
  float sol=fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
  return r_8(sol);
}


//++
void  LocalMap::SetOrigin(float theta0, float phi0, float angle)

//    
//--
{
  theta0_=theta0;
  phi0_=phi0;
  x0_=mSzX_/2;
  y0_=mSzY_/2;
  angle=angle*Pi/180.;
  cos_angle_=cos(angle);
  sin_angle_=sin(angle);
  originFlag_ =1; 
}

 
//++
void  LocalMap::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)

//    
//--
{
  theta0_=theta0;
  phi0_=phi0;
  x0_=x0;
  y0_=y0;
  angle=angle*Pi/180.;
  cos_angle_=cos(angle);
  sin_angle_=sin(angle);
  originFlag_ =1; 
} 

//++
void  LocalMap::SetSize(float angleX, float angleY)

//    
//--
{
 tgAngleX_=r_8(tan(angleX*Pi/360.));
 tgAngleY_=r_8(tan(angleY*Pi/360.));
 SzFlag_=1;
}
//++
void  LocalMap::Project(SphericalMap& sphere) const

//    Projection to spherical map   
//--
{
  for (int m=0; m<mNPix_; m++) {
    float theta,phi;
    PixThetaPhi(m,theta,phi);
    sphere(theta,phi)=mPix_[m];
    //    cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
  }
}

void LocalMap::InitNul()
//
//    initialise à zéro certaines variables internes
{ 
originFlag_=0;
SzFlag_=0;
cos_angle_=1.;
sin_angle_=0.;
mPix_=NULL;
}

void LocalMap::Clear()
//
//    
{ 
  if (mPix_)   delete[] mPix_;
}

void    LocalMap::Getij(int k, int& i, int& j) const
{
  i=(k+1)%mSzX_-1;
  if (i==-1) i=mSzX_-1;
  j=(k-i+2)/mSzX_;
}

void  LocalMap::ReferenceToUser(float &theta, float &phi) const

//       
{
  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));

    throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
  } 
  //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
  theta=(r_4)theta0_+theta-(r_4)Pi*0.5;
  if (theta<0.) {
    theta=-theta;
    phi+=(r_4)Pi;
  } else 
    if (theta>(r_4)Pi) {
      theta=2.*(r_4)Pi-theta;
      phi+=(r_4)Pi;
    } 
  
  phi=(r_4)phi0_+phi;
  while (phi>=2.*(r_4)Pi) phi-=2.*(r_4)Pi;
  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
    cout << " theta= " << theta << " phi= " << phi << endl;
    exit(0);
  }
  //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
}

void  LocalMap::UserToReference(float &theta, float &phi) const

//       
{
  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));

    throw "LocalMap::UserToReference (theta,phi) out of range"; 
  } 
  float phi1=phi-(r_4)phi0_;
  if (phi1<0.) phi1+=2.*(r_4)Pi;

  float theta1=theta-(r_4)theta0_+(r_4)Pi*0.5;
  if (theta1<0.) {
    theta=-theta1;
    phi1+=(r_4)Pi;
  } else 
    if (theta1>(r_4)Pi) {
      theta=2.*(r_4)Pi-theta1;
      phi1+=(r_4)Pi;
    }
  while (phi1>=2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
  phi=phi1;
  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
    cout << " theta= " << theta << " phi= " << phi << endl;
    exit(0);
  }
}

void LocalMap::PixProjToAngle(float x, float y,float &theta, float &phi) const {

// (x,y) representent les coordonnees en unites de pixels d'un point DANS 
// LE PLAN DE REFERENCE. On recupere (theta,phi) par rapport au repere "absolu"
// theta=pi/2,phi=0.

theta=(r_4)Pi*0.5-atan(2*y* tgAngleY_/(float) mSzY_);
phi=atan2(2*x* tgAngleX_,(float) mSzX_);
if( phi< 0. ) phi+= DeuxPi;
}

void  LocalMap::AngleProjToPix(float theta, float phi, float& x, float& y) const {

// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere 
// (i,j)  DANS LE PLAN DE REFERENCE. 
  if (phi>=(r_4)Pi) phi-= DeuxPi;
  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
  y=0.5*mSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
  x=0.5*mSzX_*tan(phi)/tgAngleX_;
}