Changeset 682 in Sophya for trunk/SophyaLib/Samba/localmap.cc


Ignore:
Timestamp:
Dec 10, 1999, 5:56:03 PM (26 years ago)
Author:
ansari
Message:

Compilation Mac pour CodeWarrior PRO 5

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaLib/Samba/localmap.cc

    r604 r682  
    66#include <string.h>
    77#include <iostream.h>
     8
     9
     10//*****************************************************************************
     11//++
     12// Class        LocalMap
     13//
     14// include      localmap.h nbmath.h
     15//
     16//    A local map of a region of the sky, in cartesian coordinates.
     17//    It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
     18//    (x0, y0 might be outside of this local map)
     19//    default value of (x0, y0) is middle of the map, center of
     20//    pixel(nx/2, ny/2)
     21//
     22//    A local map is a 2 dimensional array, with i as column index and j
     23//    as row index. The map is supposed to lie on a plan tangent to the
     24//    celestial sphere in a point whose coordinates are (x0,y0) on the local
     25//    map and (theta0, phi0) on the sphere. The range of the map is defined
     26//    by two values of angles covered respectively by all the pixels in
     27//    x direction and all the pixels in y direction (SetSize()).
     28//
     29//    A "reference plane" is considered : this plane is tangent to the
     30//    celestial sphere in a point with angles theta=Pi/2 and phi=0. This
     31//    point is the origine of coordinates is of the reference plane. The
     32//    x-axis is the tangent parallel to the equatorial line and oriented
     33//    toward the increasing phi's ; the y-axis is parallel to the meridian
     34//    line and oriented toward the north pole.
     35//   
     36//    Internally, a map is first defined within this reference plane and
     37//    tranported until the point (theta0, phi0) in such a way that both
     38//    axes are kept parallel to meridian and parallel lines of the sphere.
     39//    The user can define its own map with axes rotated with respect to
     40//    reference axes (this rotation is characterized by angle between
     41//    the local parallel line and the wanted x-axis-- see method
     42//    SetOrigin(...))
     43//   
     44//
     45
     46//
     47//--
     48//++
     49//
     50// Links        Parents
     51//
     52//    PixelMap
     53//   
     54//--
     55//++
     56//
     57// Links        Childs
     58//
     59//     
     60//--
     61//++
     62// Titre        Constructors
     63//--
     64//++
     65template<class T>
     66LocalMap<T>::LocalMap()
     67//
     68//     
     69//--
     70{
     71  InitNul();
     72  pixels_.Reset();
     73}
     74
     75//++
     76template<class T>
     77LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
     78//
     79//   
     80//--
     81{
     82  InitNul();
     83  nPix_= nx*ny;
     84  pixels_.ReSize(nPix_);
     85  pixels_.Reset();
     86}
     87
     88//++
     89template<class T>
     90LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)
     91  : pixels_(lm.pixels_, share)
     92
     93//
     94//    copy constructor
     95//--
     96{
     97  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
     98
     99  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
     100  nSzX_= lm.nSzX_;
     101  nSzY_= lm.nSzY_;
     102  nPix_= lm.nPix_;
     103  originFlag_= lm.originFlag_;
     104  extensFlag_= lm.extensFlag_;
     105  x0_= lm.x0_;
     106  y0_= lm.y0_;
     107  theta0_= lm.theta0_;
     108  phi0_= lm.phi0_;
     109  angle_= lm.angle_;
     110  cos_angle_= lm.cos_angle_;
     111  sin_angle_= lm.sin_angle_;
     112  angleX_= lm.angleX_;
     113  angleY_= lm.angleY_;
     114  tgAngleX_= lm.tgAngleX_;
     115  tgAngleY_= lm.tgAngleY_;
     116}
     117
     118//++
     119// Titre        Destructor
     120//--
     121//++
     122template<class T>
     123LocalMap<T>::~LocalMap()
     124//
     125//--
     126{
     127  InitNul();
     128}
     129
     130
     131
     132//++
     133// Titre        Public Methods
     134//--
     135
     136//++
     137template<class T>
     138void LocalMap<T>::ReSize(int_4 nx, int_4 ny)
     139//
     140//    Resize storage area for pixels
     141//--
     142{
     143  InitNul();
     144  nSzX_ = nx;
     145  nSzY_ = ny;
     146  nPix_= nx*ny;
     147  pixels_.ReSize(nPix_);
     148  pixels_.Reset();
     149}
     150
     151//++
     152template<class T>
     153int_4 LocalMap<T>::NbPixels() const
     154//
     155//    Return number of pixels
     156//--
     157{
     158  return(nPix_);
     159}
     160
     161//++
     162template<class T>
     163T& LocalMap<T>::PixVal(int_4 k)
     164//
     165//    Return value of pixel with index k
     166//--
     167{
     168  if((k < 0) || (k >= nPix_))
     169    {
     170      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
     171      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     172      THROW(rangeCheckErr);
     173      //throw "LocalMap::PIxVal Pixel index out of range";
     174    }
     175  return(pixels_(k));
     176}
     177
     178//++
     179
     180template<class T>
     181T const& LocalMap<T>::PixVal(int_4 k) const
     182//
     183//   const version of previous method
     184//--
     185{
     186  if((k < 0) || (k >= nPix_))
     187    {
     188      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
     189      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     190     
     191      throw "LocalMap::PIxVal Pixel index out of range";
     192    }
     193  return *(pixels_.Data()+k);
     194}
     195
     196template<class T>
     197bool LocalMap<T>::ContainsSph(double theta, double phi) const
     198{
     199// $CHECK$  Reza 11/11/99 - A modifier
     200return(true);
     201}
     202
     203//++
     204template<class T>
     205int_4 LocalMap<T>::PixIndexSph(double theta,double phi) const
     206//
     207//    Return index of the pixel with spherical coordinates (theta,phi)
     208//--
     209{
     210  int_4 i,j;
     211  if(!(originFlag_) || !(extensFlag_))
     212    {
     213      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
     214      exit(0);
     215    }
     216
     217  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
     218  double theta_aux= theta;
     219  double phi_aux  = phi;
     220  UserToReference(theta_aux, phi_aux);
     221
     222  // coordonnees dans le plan local en unites de pixels
     223  double x,y;
     224  AngleProjToPix(theta_aux,phi_aux, x, y);
     225
     226  double xmin= -x0_-0.5;
     227  double xmax= xmin+nSzX_;
     228  if((x > xmax) || (x < xmin)) return(-1);
     229  double xcurrent= xmin;
     230  for(i = 0; i < nSzX_; i++ )
     231    {
     232      xcurrent += 1.;
     233      if( x < xcurrent ) break;
     234    }
     235  double ymin= -y0_-0.5;
     236  double ymax= ymin+nSzY_;
     237  if((y >  ymax) || (y < ymin)) return(-1);
     238  double ycurrent= ymin;
     239  for(j = 0; j < nSzY_; j++ )
     240    {
     241      ycurrent += 1.;
     242      if( y < ycurrent ) break;
     243    }
     244  return (j*nSzX_+i);
     245}
     246
     247//++
     248template<class T>
     249void LocalMap<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
     250//
     251//    Return (theta, phi) coordinates of pixel with index k
     252//--
     253{
     254  if(!(originFlag_) || !(extensFlag_))
     255    {
     256      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
     257      exit(0);
     258    }
     259
     260  int_4 i,j;
     261  Getij(k,i,j);
     262
     263  double X= double(i-x0_);
     264  double Y= double(j-y0_);
     265  // situation de ce pixel dans le plan de reference
     266  double x= X*cos_angle_-Y*sin_angle_;
     267  double y= X*sin_angle_+Y* cos_angle_;
     268  // projection sur la sphere
     269  PixProjToAngle(x, y, theta, phi);
     270  // retour au plan utilisateur
     271  ReferenceToUser(theta, phi);
     272}
     273
     274template <class T>
     275T LocalMap<T>::SetPixels(T v)
     276{
     277pixels_.Reset(v);
     278return(v);
     279}
     280 
     281//++
     282template<class T>
     283double LocalMap<T>::PixSolAngle(int_4 k) const
     284//
     285//    Pixel Solid angle  (steradians)
     286//    All the pixels have not necessarly the same size in (theta, phi)
     287//    because of the projection scheme which is not yet fixed.   
     288//--
     289{
     290  int_4 i,j;
     291  Getij(k,i,j);
     292  double X= double(i-x0_);
     293  double Y= double(j-y0_);
     294  double XR= X+double(i)*0.5;
     295  double XL= X-double(i)*0.5;
     296  double YU= Y+double(j)*0.5;
     297  double YL= Y-double(j)*0.5;
     298
     299  // situation  dans le plan de reference
     300  double x0= XL*cos_angle_-YL*sin_angle_;
     301  double y0= XL*sin_angle_+YL*cos_angle_;
     302  double xa= XR*cos_angle_-YL*sin_angle_;
     303  double ya= XR*sin_angle_+YL*cos_angle_;
     304  double xb= XL*cos_angle_-YU*sin_angle_;
     305  double yb= XL*sin_angle_+YU*cos_angle_;
     306
     307  // projection sur la sphere
     308  double thet0,phi0,theta,phia,thetb,phib;
     309  PixProjToAngle(x0, y0, thet0, phi0);
     310  PixProjToAngle(xa, ya, theta, phia);
     311  PixProjToAngle(xb, yb, thetb, phib);
     312
     313  // angle solide
     314  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
     315  return sol;
     316}
     317
     318//++
     319template<class T>
     320void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
     321//
     322//    set the referential of the map (angles in degrees)
     323//    (default x0=siz_x/2,  y0=siz_y/2)
     324//--
     325{
     326  theta0_= theta0;
     327  phi0_  = phi0;
     328  angle_ = angle;
     329  x0_= nSzX_/2;
     330  y0_= nSzY_/2;
     331  cos_angle_= cos(angle*Pi/180.);
     332  sin_angle_= sin(angle*Pi/180.);
     333  originFlag_= true;
     334  cout << " LocalMap:: set origin 1 done" << endl;
     335}
     336 
     337//++
     338template<class T>
     339void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)
     340//
     341//    set the referential of the map (angles in degrees)
     342//--
     343{
     344  theta0_= theta0;
     345  phi0_  = phi0;
     346  angle_ = angle;
     347  x0_= x0;
     348  y0_= y0;
     349  cos_angle_= cos(angle*Pi/180.);
     350  sin_angle_= sin(angle*Pi/180.);
     351  originFlag_= true;
     352  cout << " LocalMap:: set origin 2 done" << endl;
     353}
     354
     355//++
     356template<class T>
     357void LocalMap<T>::SetSize(double angleX,double angleY)
     358//
     359//    angle range of tthe map (angles in degrees)
     360//--
     361{
     362  angleX_= angleX;
     363  angleY_= angleY;
     364
     365  // tangente de la moitie de l'ouverture angulaire totale
     366  tgAngleX_= tan(0.5*angleX_*Pi/180.);
     367  tgAngleY_= tan(0.5*angleY_*Pi/180.);
     368
     369  extensFlag_= true;
     370  cout << " LocalMap:: set extension done" << endl;
     371}
     372
     373//++
     374template<class T>
     375void LocalMap<T>::Project(SphericalMap<T>& sphere) const
     376//
     377//    Projection to a spherical map   
     378//--
     379{
     380  for(int m = 0; m < nPix_; m++)
     381    {
     382      double theta,phi;
     383      PixThetaPhi(m,theta,phi);
     384      sphere(theta,phi)= pixels_(m);
     385      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
     386    }
     387}
     388// Titre        Private Methods
     389//++
     390template<class T>
     391void LocalMap<T>::InitNul()
     392//
     393//    set some attributes to zero
     394//--
     395{
     396  originFlag_= false;
     397  extensFlag_= false;
     398  cos_angle_= 1.0;
     399  sin_angle_= 0.0;
     400//  pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
     401}
     402
     403//++
     404template<class T>
     405void LocalMap<T>::Getij(int_4 k,int_4& i,int_4& j) const
     406//
     407//    Return 2 indices corresponding to the pixel number k
     408//--
     409{
     410  i= (k+1)%nSzX_-1;
     411  if(i == -1) i= nSzX_-1;
     412  j= (k-i+2)/nSzX_;
     413}
     414
     415//++
     416template<class T>
     417void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
     418//
     419//    Transform a pair of coordinates (theta, phi) given in
     420//    reference coordinates into map coordinates
     421//--     
     422{
     423  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     424    {
     425      throw "LocalMap::ReferenceToUser (theta,phi) out of range";
     426    }
     427
     428  theta= theta0_*Pi/180.+theta-Pi*0.5;
     429  if(theta < 0.)
     430    {
     431      theta= -theta;
     432      phi += Pi;
     433    }
     434  else
     435    {
     436      if(theta > Pi)
     437        {
     438          theta= 2.*Pi-theta;
     439          phi += Pi;
     440        }
     441    }
     442
     443  phi= phi0_*Pi/180.+phi;
     444  while(phi >= 2.*Pi) phi-= 2.*Pi;
     445
     446  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     447    {
     448      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
     449      cout << " theta= " << theta << " phi= " << phi << endl;
     450      exit(0);
     451    }
     452}
     453
     454//++
     455template<class T>
     456void  LocalMap<T>::UserToReference(double& theta,double& phi) const
     457//
     458//    Transform a pair of coordinates (theta, phi) given in
     459//    map coordinates into reference coordinates
     460//--     
     461{
     462  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     463    {
     464      cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;
     465      // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     466      throw "LocalMap::UserToReference (theta,phi) out of range";
     467    }
     468
     469  double phi1= phi-phi0_*Pi/180.;
     470  if(phi1 < 0.) phi1+= 2.*Pi;
     471
     472  double theta1= theta-theta0_*Pi/180.+Pi*0.5;
     473  if(theta1 < 0.)
     474    {
     475      theta= -theta1;
     476      phi1+= Pi;
     477    }
     478  else
     479    {
     480      if(theta1 > Pi)
     481        {
     482          theta= 2.*Pi-theta1;
     483          phi1+= Pi;
     484        }
     485    }
     486
     487  while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
     488  phi= phi1;
     489  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     490    {
     491      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
     492      cout << " theta= " << theta << " phi= " << phi << endl;
     493      exit(0);
     494    }
     495}
     496
     497//++
     498template<class T>
     499void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const
     500//
     501// 
     502//    Given coordinates in pixel units in the REFERENCE PLANE, return
     503//    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
     504//--
     505{
     506  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
     507  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
     508  if(phi < 0.) phi += DeuxPi;
     509}
     510
     511//++
     512template<class T>
     513void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
     514//
     515//    Given coordinates  (theta, phi) in "absolute" referential
     516//    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
     517//--
     518{
     519  if(phi >= Pi) phi-= DeuxPi;
     520  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
     521  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ? 
     522  x= 0.5*nSzX_*tan(phi)/tgAngleX_;
     523}
     524
     525template<class T>
     526void LocalMap<T>::print(ostream& os) const
     527{
     528  os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
     529  if(LocalMap_isDone())
     530    {
     531      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
     532      os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
     533      os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
     534      os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
     535      os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
     536    }
     537
     538  os << " contenu de pixels : ";
     539  for(int i=0; i < nPix_; i++)
     540    {
     541      if(i%5 == 0) os << endl;
     542      os <<  pixels_(i) <<", ";
     543    }
     544  os << endl;
     545}
     546//++
     547// Titre        class FIO_LocalMap
     548//    Delegated objects for persitance management
     549//--
     550
     551//*******************************************************************
     552// class FIO_LocalMap<T>
     553//  Les objets delegues pour la gestion de persistance
     554//*******************************************************************
     555
     556//++
     557template <class T>
     558FIO_LocalMap<T>::FIO_LocalMap()
     559//
     560//--
     561{
     562  dobj= new LocalMap<T>;
     563  ownobj= true;
     564}
     565//++
     566template <class T>
     567FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
     568//
     569//--
     570{
     571  dobj= new LocalMap<T>;
     572  dobj->DataBlock().SetTemp(true);
     573  ownobj= true;
     574  Read(filename);
     575}
     576
     577//++
     578template <class T>
     579FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
     580//
     581//--
     582{
     583  dobj= new LocalMap<T>(obj, true);
     584  dobj->DataBlock().SetTemp(true);
     585  ownobj= true;
     586}
     587
     588template <class T>
     589FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
     590{
     591  dobj= obj;
     592  ownobj= false;
     593}
     594
     595//++
     596template <class T>
     597FIO_LocalMap<T>::~FIO_LocalMap()
     598//
     599//--
     600{
     601  if (ownobj && dobj) delete dobj;
     602}
     603//++
     604template <class T>
     605AnyDataObj* FIO_LocalMap<T>::DataObj()
     606//
     607//--
     608{
     609  return(dobj);
     610}
     611
     612//++
     613template <class T>
     614void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
     615//
     616//--
     617{
     618
     619  if(dobj == NULL)
     620    {
     621      dobj= new LocalMap<T>;
     622      dobj->DataBlock().SetTemp(true);
     623      ownobj= true;
     624    }
     625
     626  // Pour savoir s'il y avait un DVList Info associe
     627  char strg[256];
     628  is.GetLine(strg, 255);
     629  bool hadinfo= false;
     630  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
     631  if(hadinfo)
     632    {    // Lecture eventuelle du DVList Info
     633      is >> dobj->Info();
     634    }
     635
     636  int_4 nSzX;
     637  is.GetI4(nSzX);
     638  dobj->setSize_x(nSzX);
     639
     640  int_4 nSzY;
     641  is.GetI4(nSzY);
     642  dobj->setSize_y(nSzY);
     643
     644  int_4 nPix;
     645  is.GetI4(nPix);
     646  dobj->setNbPixels(nPix);
     647
     648  string ss("local mapping is done");
     649  string sso;
     650  is.GetStr(sso);
     651  if(sso == ss)
     652    {
     653      cout<<" ReadSelf:: local mapping"<<endl;
     654      int_4 x0, y0;
     655      double theta, phi, angle;
     656      is.GetI4(x0);
     657      is.GetI4(y0);
     658      is.GetR8(theta);
     659      is.GetR8(phi);
     660      is.GetR8(angle);
     661      dobj->SetOrigin(theta, phi, x0, y0, angle);
     662
     663      double angleX, angleY;
     664      is.GetR8(angleX);
     665      is.GetR8(angleY);
     666      dobj->SetSize(angleX, angleY);
     667    }
     668
     669// On lit le DataBlock;
     670  is >> dobj->DataBlock();
     671}
     672
     673//++
     674template <class T>
     675void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
     676//
     677//--
     678{
     679  if(dobj == NULL)
     680    {
     681      cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;
     682      return;
     683    }
     684
     685  char strg[256];
     686  int_4 nSzX= dobj->Size_x();
     687  int_4 nSzY= dobj->Size_y();
     688  int_4 nPix= dobj->NbPixels();
     689 
     690  if(dobj->ptrInfo())
     691    {
     692      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
     693      os.PutLine(strg);
     694      os << dobj->Info();
     695    }
     696  else
     697    {
     698      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
     699      os.PutLine(strg); 
     700    }
     701
     702  os.PutI4(nSzX);
     703  os.PutI4(nSzY);
     704  os.PutI4(nPix);
     705
     706  if(dobj->LocalMap_isDone())
     707    {
     708      string ss("local mapping is done");
     709      os.PutStr(ss);
     710      int_4 x0, y0;
     711      double theta, phi, angle;
     712      dobj->Origin(theta, phi, x0, y0, angle);
     713      os.PutI4(x0);
     714      os.PutI4(y0);
     715      os.PutR8(theta);
     716      os.PutR8(phi);
     717      os.PutR8(angle);
     718
     719      double angleX, angleY;
     720      dobj->Aperture(angleX, angleY);
     721      os.PutR8(angleX);
     722      os.PutR8(angleY);
     723    }
     724  else
     725    {
     726      string ss("no local mapping");
     727      os.PutStr(ss);
     728    }
     729
     730// On ecrit le dataBlock
     731  os << dobj->DataBlock();
     732}
    8733
    9734#ifdef __CXX_PRAGMA_TEMPLATES__
     
    27752template class FIO_LocalMap< complex<double> >;
    28753#endif
    29 
    30 //*****************************************************************************
    31 //++
    32 // Class        LocalMap
    33 //
    34 // include      localmap.h nbmath.h
    35 //
    36 //    A local map of a region of the sky, in cartesian coordinates.
    37 //    It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
    38 //    (x0, y0 might be outside of this local map)
    39 //    default value of (x0, y0) is middle of the map, center of
    40 //    pixel(nx/2, ny/2)
    41 //
    42 //    A local map is a 2 dimensional array, with i as column index and j
    43 //    as row index. The map is supposed to lie on a plan tangent to the
    44 //    celestial sphere in a point whose coordinates are (x0,y0) on the local
    45 //    map and (theta0, phi0) on the sphere. The range of the map is defined
    46 //    by two values of angles covered respectively by all the pixels in
    47 //    x direction and all the pixels in y direction (SetSize()).
    48 //
    49 //    A "reference plane" is considered : this plane is tangent to the
    50 //    celestial sphere in a point with angles theta=Pi/2 and phi=0. This
    51 //    point is the origine of coordinates is of the reference plane. The
    52 //    x-axis is the tangent parallel to the equatorial line and oriented
    53 //    toward the increasing phi's ; the y-axis is parallel to the meridian
    54 //    line and oriented toward the north pole.
    55 //   
    56 //    Internally, a map is first defined within this reference plane and
    57 //    tranported until the point (theta0, phi0) in such a way that both
    58 //    axes are kept parallel to meridian and parallel lines of the sphere.
    59 //    The user can define its own map with axes rotated with respect to
    60 //    reference axes (this rotation is characterized by angle between
    61 //    the local parallel line and the wanted x-axis-- see method
    62 //    SetOrigin(...))
    63 //   
    64 //
    65 
    66 //
    67 //--
    68 //++
    69 //
    70 // Links        Parents
    71 //
    72 //    PixelMap
    73 //   
    74 //--
    75 //++
    76 //
    77 // Links        Childs
    78 //
    79 //     
    80 //--
    81 //++
    82 // Titre        Constructors
    83 //--
    84 //++
    85 template<class T>
    86 LocalMap<T>::LocalMap()
    87 //
    88 //     
    89 //--
    90 {
    91   InitNul();
    92   pixels_.Reset();
    93 }
    94 
    95 //++
    96 template<class T>
    97 LocalMap<T>::LocalMap(int nx, int ny) : nSzX_(nx), nSzY_(ny)
    98 //
    99 //   
    100 //--
    101 {
    102   InitNul();
    103   nPix_= nx*ny;
    104   pixels_.ReSize(nPix_);
    105   pixels_.Reset();
    106 }
    107 
    108 //++
    109 template<class T>
    110 LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)
    111   : pixels_(lm.pixels_, share)
    112 
    113 //
    114 //    copy constructor
    115 //--
    116 {
    117   cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
    118 
    119   if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
    120   nSzX_= lm.nSzX_;
    121   nSzY_= lm.nSzY_;
    122   nPix_= lm.nPix_;
    123   originFlag_= lm.originFlag_;
    124   extensFlag_= lm.extensFlag_;
    125   x0_= lm.x0_;
    126   y0_= lm.y0_;
    127   theta0_= lm.theta0_;
    128   phi0_= lm.phi0_;
    129   angle_= lm.angle_;
    130   cos_angle_= lm.cos_angle_;
    131   sin_angle_= lm.sin_angle_;
    132   angleX_= lm.angleX_;
    133   angleY_= lm.angleY_;
    134   tgAngleX_= lm.tgAngleX_;
    135   tgAngleY_= lm.tgAngleY_;
    136 }
    137 
    138 //++
    139 // Titre        Destructor
    140 //--
    141 //++
    142 template<class T>
    143 LocalMap<T>::~LocalMap()
    144 //
    145 //--
    146 {
    147   InitNul();
    148 }
    149 
    150 
    151 
    152 //++
    153 // Titre        Public Methods
    154 //--
    155 
    156 //++
    157 template<class T>
    158 void LocalMap<T>::ReSize(int nx, int ny)
    159 //
    160 //    Resize storage area for pixels
    161 //--
    162 {
    163   InitNul();
    164   nSzX_ = nx;
    165   nSzY_ = ny;
    166   nPix_= nx*ny;
    167   pixels_.ReSize(nPix_);
    168   pixels_.Reset();
    169 }
    170 
    171 //++
    172 template<class T>
    173 int LocalMap<T>::NbPixels() const
    174 //
    175 //    Return number of pixels
    176 //--
    177 {
    178   return(nPix_);
    179 }
    180 
    181 //++
    182 template<class T>
    183 T& LocalMap<T>::PixVal(int k)
    184 //
    185 //    Return value of pixel with index k
    186 //--
    187 {
    188   if((k < 0) || (k >= nPix_))
    189     {
    190       cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
    191       //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    192       THROW(rangeCheckErr);
    193       //throw "LocalMap::PIxVal Pixel index out of range";
    194     }
    195   return(pixels_(k));
    196 }
    197 
    198 //++
    199 
    200 template<class T>
    201 T const& LocalMap<T>::PixVal(int k) const
    202 //
    203 //   const version of previous method
    204 //--
    205 {
    206   if((k < 0) || (k >= nPix_))
    207     {
    208       cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
    209       //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    210      
    211       throw "LocalMap::PIxVal Pixel index out of range";
    212     }
    213   return *(pixels_.Data()+k);
    214 }
    215 
    216 template<class T>
    217 bool LocalMap<T>::ContainsSph(double theta, double phi) const
    218 {
    219 // $CHECK$  Reza 11/11/99 - A modifier
    220 return(true);
    221 }
    222 
    223 //++
    224 template<class T>
    225 int LocalMap<T>::PixIndexSph(double theta,double phi) const
    226 //
    227 //    Return index of the pixel with spherical coordinates (theta,phi)
    228 //--
    229 {
    230   int i,j;
    231   if(!(originFlag_) || !(extensFlag_))
    232     {
    233       cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    234       exit(0);
    235     }
    236 
    237   // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
    238   double theta_aux= theta;
    239   double phi_aux  = phi;
    240   UserToReference(theta_aux, phi_aux);
    241 
    242   // coordonnees dans le plan local en unites de pixels
    243   double x,y;
    244   AngleProjToPix(theta_aux,phi_aux, x, y);
    245 
    246   double xmin= -x0_-0.5;
    247   double xmax= xmin+nSzX_;
    248   if((x > xmax) || (x < xmin)) return(-1);
    249   double xcurrent= xmin;
    250   for(i = 0; i < nSzX_; i++ )
    251     {
    252       xcurrent += 1.;
    253       if( x < xcurrent ) break;
    254     }
    255   double ymin= -y0_-0.5;
    256   double ymax= ymin+nSzY_;
    257   if((y >  ymax) || (y < ymin)) return(-1);
    258   double ycurrent= ymin;
    259   for(j = 0; j < nSzY_; j++ )
    260     {
    261       ycurrent += 1.;
    262       if( y < ycurrent ) break;
    263     }
    264   return (j*nSzX_+i);
    265 }
    266 
    267 //++
    268 template<class T>
    269 void LocalMap<T>::PixThetaPhi(int k,double& theta,double& phi) const
    270 //
    271 //    Return (theta, phi) coordinates of pixel with index k
    272 //--
    273 {
    274   if(!(originFlag_) || !(extensFlag_))
    275     {
    276       cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    277       exit(0);
    278     }
    279 
    280   int i,j;
    281   Getij(k,i,j);
    282 
    283   double X= double(i-x0_);
    284   double Y= double(j-y0_);
    285   // situation de ce pixel dans le plan de reference
    286   double x= X*cos_angle_-Y*sin_angle_;
    287   double y= X*sin_angle_+Y* cos_angle_;
    288   // projection sur la sphere
    289   PixProjToAngle(x, y, theta, phi);
    290   // retour au plan utilisateur
    291   ReferenceToUser(theta, phi);
    292 }
    293 
    294 template <class T>
    295 T LocalMap<T>::SetPixels(T v)
    296 {
    297 pixels_.Reset(v);
    298 return(v);
    299 }
    300  
    301 //++
    302 template<class T>
    303 double LocalMap<T>::PixSolAngle(int k) const
    304 //
    305 //    Pixel Solid angle  (steradians)
    306 //    All the pixels have not necessarly the same size in (theta, phi)
    307 //    because of the projection scheme which is not yet fixed.   
    308 //--
    309 {
    310   int i,j;
    311   Getij(k,i,j);
    312   double X= double(i-x0_);
    313   double Y= double(j-y0_);
    314   double XR= X+double(i)*0.5;
    315   double XL= X-double(i)*0.5;
    316   double YU= Y+double(j)*0.5;
    317   double YL= Y-double(j)*0.5;
    318 
    319   // situation  dans le plan de reference
    320   double x0= XL*cos_angle_-YL*sin_angle_;
    321   double y0= XL*sin_angle_+YL*cos_angle_;
    322   double xa= XR*cos_angle_-YL*sin_angle_;
    323   double ya= XR*sin_angle_+YL*cos_angle_;
    324   double xb= XL*cos_angle_-YU*sin_angle_;
    325   double yb= XL*sin_angle_+YU*cos_angle_;
    326 
    327   // projection sur la sphere
    328   double thet0,phi0,theta,phia,thetb,phib;
    329   PixProjToAngle(x0, y0, thet0, phi0);
    330   PixProjToAngle(xa, ya, theta, phia);
    331   PixProjToAngle(xb, yb, thetb, phib);
    332 
    333   // angle solide
    334   double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
    335   return sol;
    336 }
    337 
    338 //++
    339 template<class T>
    340 void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
    341 //
    342 //    set the referential of the map (angles in degrees)
    343 //    (default x0=siz_x/2,  y0=siz_y/2)
    344 //--
    345 {
    346   theta0_= theta0;
    347   phi0_  = phi0;
    348   angle_ = angle;
    349   x0_= nSzX_/2;
    350   y0_= nSzY_/2;
    351   cos_angle_= cos(angle*Pi/180.);
    352   sin_angle_= sin(angle*Pi/180.);
    353   originFlag_= true;
    354   cout << " LocalMap:: set origin 1 done" << endl;
    355 }
    356  
    357 //++
    358 template<class T>
    359 void LocalMap<T>::SetOrigin(double theta0,double phi0,int x0,int y0,double angle)
    360 //
    361 //    set the referential of the map (angles in degrees)
    362 //--
    363 {
    364   theta0_= theta0;
    365   phi0_  = phi0;
    366   angle_ = angle;
    367   x0_= x0;
    368   y0_= y0;
    369   cos_angle_= cos(angle*Pi/180.);
    370   sin_angle_= sin(angle*Pi/180.);
    371   originFlag_= true;
    372   cout << " LocalMap:: set origin 2 done" << endl;
    373 }
    374 
    375 //++
    376 template<class T>
    377 void LocalMap<T>::SetSize(double angleX,double angleY)
    378 //
    379 //    angle range of tthe map (angles in degrees)
    380 //--
    381 {
    382   angleX_= angleX;
    383   angleY_= angleY;
    384 
    385   // tangente de la moitie de l'ouverture angulaire totale
    386   tgAngleX_= tan(0.5*angleX_*Pi/180.);
    387   tgAngleY_= tan(0.5*angleY_*Pi/180.);
    388 
    389   extensFlag_= true;
    390   cout << " LocalMap:: set extension done" << endl;
    391 }
    392 
    393 //++
    394 template<class T>
    395 void LocalMap<T>::Project(SphericalMap<T>& sphere) const
    396 //
    397 //    Projection to a spherical map   
    398 //--
    399 {
    400   for(int m = 0; m < nPix_; m++)
    401     {
    402       double theta,phi;
    403       PixThetaPhi(m,theta,phi);
    404       sphere(theta,phi)= pixels_(m);
    405       // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
    406     }
    407 }
    408 // Titre        Private Methods
    409 //++
    410 template<class T>
    411 void LocalMap<T>::InitNul()
    412 //
    413 //    set some attributes to zero
    414 //--
    415 {
    416   originFlag_= false;
    417   extensFlag_= false;
    418   cos_angle_= 1.0;
    419   sin_angle_= 0.0;
    420 //  pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
    421 }
    422 
    423 //++
    424 template<class T>
    425 void LocalMap<T>::Getij(int k,int& i,int& j) const
    426 //
    427 //    Return 2 indices corresponding to the pixel number k
    428 //--
    429 {
    430   i= (k+1)%nSzX_-1;
    431   if(i == -1) i= nSzX_-1;
    432   j= (k-i+2)/nSzX_;
    433 }
    434 
    435 //++
    436 template<class T>
    437 void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
    438 //
    439 //    Transform a pair of coordinates (theta, phi) given in
    440 //    reference coordinates into map coordinates
    441 //--     
    442 {
    443   if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    444     {
    445       throw "LocalMap::ReferenceToUser (theta,phi) out of range";
    446     }
    447 
    448   theta= theta0_*Pi/180.+theta-Pi*0.5;
    449   if(theta < 0.)
    450     {
    451       theta= -theta;
    452       phi += Pi;
    453     }
    454   else
    455     {
    456       if(theta > Pi)
    457         {
    458           theta= 2.*Pi-theta;
    459           phi += Pi;
    460         }
    461     }
    462 
    463   phi= phi0_*Pi/180.+phi;
    464   while(phi >= 2.*Pi) phi-= 2.*Pi;
    465 
    466   if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    467     {
    468       cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
    469       cout << " theta= " << theta << " phi= " << phi << endl;
    470       exit(0);
    471     }
    472 }
    473 
    474 //++
    475 template<class T>
    476 void  LocalMap<T>::UserToReference(double& theta,double& phi) const
    477 //
    478 //    Transform a pair of coordinates (theta, phi) given in
    479 //    map coordinates into reference coordinates
    480 //--     
    481 {
    482   if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    483     {
    484       cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;
    485       // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    486       throw "LocalMap::UserToReference (theta,phi) out of range";
    487     }
    488 
    489   double phi1= phi-phi0_*Pi/180.;
    490   if(phi1 < 0.) phi1+= 2.*Pi;
    491 
    492   double theta1= theta-theta0_*Pi/180.+Pi*0.5;
    493   if(theta1 < 0.)
    494     {
    495       theta= -theta1;
    496       phi1+= Pi;
    497     }
    498   else
    499     {
    500       if(theta1 > Pi)
    501         {
    502           theta= 2.*Pi-theta1;
    503           phi1+= Pi;
    504         }
    505     }
    506 
    507   while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
    508   phi= phi1;
    509   if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    510     {
    511       cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
    512       cout << " theta= " << theta << " phi= " << phi << endl;
    513       exit(0);
    514     }
    515 }
    516 
    517 //++
    518 template<class T>
    519 void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const
    520 //
    521 // 
    522 //    Given coordinates in pixel units in the REFERENCE PLANE, return
    523 //    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
    524 //--
    525 {
    526   theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
    527   phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
    528   if(phi < 0.) phi += DeuxPi;
    529 }
    530 
    531 //++
    532 template<class T>
    533 void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
    534 //
    535 //    Given coordinates  (theta, phi) in "absolute" referential
    536 //    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
    537 //--
    538 {
    539   if(phi >= Pi) phi-= DeuxPi;
    540   //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
    541   y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ? 
    542   x= 0.5*nSzX_*tan(phi)/tgAngleX_;
    543 }
    544 
    545 template<class T>
    546 void LocalMap<T>::print(ostream& os) const
    547 {
    548   os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
    549   if(LocalMap_isDone())
    550     {
    551       os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
    552       os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
    553       os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
    554       os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
    555       os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
    556     }
    557 
    558   os << " contenu de pixels : ";
    559   for(int i=0; i < nPix_; i++)
    560     {
    561       if(i%5 == 0) os << endl;
    562       os <<  pixels_(i) <<", ";
    563     }
    564   os << endl;
    565 }
    566 //++
    567 // Titre        class FIO_LocalMap
    568 //    Delegated objects for persitance management
    569 //--
    570 
    571 //*******************************************************************
    572 // class FIO_LocalMap<T>
    573 //  Les objets delegues pour la gestion de persistance
    574 //*******************************************************************
    575 
    576 //++
    577 template <class T>
    578 FIO_LocalMap<T>::FIO_LocalMap()
    579 //
    580 //--
    581 {
    582   dobj= new LocalMap<T>;
    583   ownobj= true;
    584 }
    585 //++
    586 template <class T>
    587 FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
    588 //
    589 //--
    590 {
    591   dobj= new LocalMap<T>;
    592   dobj->DataBlock().SetTemp(true);
    593   ownobj= true;
    594   Read(filename);
    595 }
    596 
    597 //++
    598 template <class T>
    599 FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
    600 //
    601 //--
    602 {
    603   dobj= new LocalMap<T>(obj, true);
    604   dobj->DataBlock().SetTemp(true);
    605   ownobj= true;
    606 }
    607 
    608 template <class T>
    609 FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
    610 {
    611   dobj= obj;
    612   ownobj= false;
    613 }
    614 
    615 //++
    616 template <class T>
    617 FIO_LocalMap<T>::~FIO_LocalMap()
    618 //
    619 //--
    620 {
    621   if (ownobj && dobj) delete dobj;
    622 }
    623 //++
    624 template <class T>
    625 AnyDataObj* FIO_LocalMap<T>::DataObj()
    626 //
    627 //--
    628 {
    629   return(dobj);
    630 }
    631 
    632 //++
    633 template <class T>
    634 void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
    635 //
    636 //--
    637 {
    638 
    639   if(dobj == NULL)
    640     {
    641       dobj= new LocalMap<T>;
    642       dobj->DataBlock().SetTemp(true);
    643       ownobj= true;
    644     }
    645 
    646   // Pour savoir s'il y avait un DVList Info associe
    647   char strg[256];
    648   is.GetLine(strg, 255);
    649   bool hadinfo= false;
    650   if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
    651   if(hadinfo)
    652     {    // Lecture eventuelle du DVList Info
    653       is >> dobj->Info();
    654     }
    655 
    656   int nSzX;
    657   is.GetI4(nSzX);
    658   dobj->setSize_x(nSzX);
    659 
    660   int nSzY;
    661   is.GetI4(nSzY);
    662   dobj->setSize_y(nSzY);
    663 
    664   int nPix;
    665   is.GetI4(nPix);
    666   dobj->setNbPixels(nPix);
    667 
    668   string ss("local mapping is done");
    669   string sso;
    670   is.GetStr(sso);
    671   if(sso == ss)
    672     {
    673       cout<<" ReadSelf:: local mapping"<<endl;
    674       int x0, y0;
    675       double theta, phi, angle;
    676       is.GetI4(x0);
    677       is.GetI4(y0);
    678       is.GetR8(theta);
    679       is.GetR8(phi);
    680       is.GetR8(angle);
    681       dobj->SetOrigin(theta, phi, x0, y0, angle);
    682 
    683       double angleX, angleY;
    684       is.GetR8(angleX);
    685       is.GetR8(angleY);
    686       dobj->SetSize(angleX, angleY);
    687     }
    688 
    689 // On lit le DataBlock;
    690   is >> dobj->DataBlock();
    691 }
    692 
    693 //++
    694 template <class T>
    695 void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
    696 //
    697 //--
    698 {
    699   if(dobj == NULL)
    700     {
    701       cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;
    702       return;
    703     }
    704 
    705   char strg[256];
    706   int nSzX= dobj->Size_x();
    707   int nSzY= dobj->Size_y();
    708   int nPix= dobj->NbPixels();
    709  
    710   if(dobj->ptrInfo())
    711     {
    712       sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
    713       os.PutLine(strg);
    714       os << dobj->Info();
    715     }
    716   else
    717     {
    718       sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
    719       os.PutLine(strg); 
    720     }
    721 
    722   os.PutI4(nSzX);
    723   os.PutI4(nSzY);
    724   os.PutI4(nPix);
    725 
    726   if(dobj->LocalMap_isDone())
    727     {
    728       string ss("local mapping is done");
    729       os.PutStr(ss);
    730       int x0, y0;
    731       double theta, phi, angle;
    732       dobj->Origin(theta, phi, x0, y0, angle);
    733       os.PutI4(x0);
    734       os.PutI4(y0);
    735       os.PutR8(theta);
    736       os.PutR8(phi);
    737       os.PutR8(angle);
    738 
    739       double angleX, angleY;
    740       dobj->Aperture(angleX, angleY);
    741       os.PutR8(angleX);
    742       os.PutR8(angleY);
    743     }
    744   else
    745     {
    746       string ss("no local mapping");
    747       os.PutStr(ss);
    748     }
    749 
    750 // On ecrit le dataBlock
    751   os << dobj->DataBlock();
    752 }
    753 
Note: See TracChangeset for help on using the changeset viewer.