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


Ignore:
Timestamp:
Oct 18, 1999, 4:37:44 PM (26 years ago)
Author:
ansari
Message:

modifs francois : passage en double etc. 18-OCT-99

File:
1 edited

Legend:

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

    r470 r473  
    8989//++
    9090template<class T>
    91 LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
     91LocalMap<T>::LocalMap(int nx, int ny) : nSzX_(nx), nSzY_(ny)
    9292//
    9393// Constructeur
     
    155155//++
    156156template<class T>
    157 void LocalMap<T>::ReSize(int_4 nx, int_4 ny)
     157void LocalMap<T>::ReSize(int nx, int ny)
    158158//
    159159// Redimensionne l'espace de stokage des pixels
     
    175175//++
    176176template<class T>
    177 int_4 LocalMap<T>::NbPixels() const
     177int LocalMap<T>::NbPixels() const
    178178//
    179179//    Retourne le nombre de pixels du découpage
     
    185185//++
    186186template<class T>
    187 T& LocalMap<T>::PixVal(int_4 k)
     187T& LocalMap<T>::PixVal(int k)
    188188//
    189189//    Retourne la valeur du contenu du pixel d'indice k
     
    203203
    204204template<class T>
    205 T const& LocalMap<T>::PixVal(int_4 k) const
     205T const& LocalMap<T>::PixVal(int k) const
    206206//
    207207//    Retourne la valeur du contenu du pixel d'indice k
     
    220220//++
    221221template<class T>
    222 int_4 LocalMap<T>::PixIndexSph(float theta, float phi) const
     222int LocalMap<T>::PixIndexSph(double theta,double phi) const
    223223//
    224224//    Retourne l'indice du pixel vers lequel pointe une direction définie par
     
    234234
    235235  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
    236   float theta_aux=theta;
    237   float phi_aux=phi;
     236  double theta_aux= theta;
     237  double phi_aux  = phi;
    238238  UserToReference(theta_aux, phi_aux);
    239239
    240240  // coordonnees dans le plan local en unites de pixels
    241   float x,y;
    242   AngleProjToPix(theta_aux, phi_aux, x, y);
    243 
    244   float xmin= -x0_-0.5;
    245   float xmax= xmin+nSzX_;
    246   if((x >  xmax) || (x < xmin)) return(-1);
    247   float xcurrent= xmin;
     241  double x,y;
     242  AngleProjToPix(theta_aux,phi_aux, x, y);
     243
     244  double xmin= -x0_-0.5;
     245  double xmax= xmin+nSzX_;
     246  if((x > xmax) || (x < xmin)) return(-1);
     247  double xcurrent= xmin;
    248248  for(i = 0; i < nSzX_; i++ )
    249249    {
     
    251251      if( x < xcurrent ) break;
    252252    }
    253   float ymin= -y0_-0.5;
    254   float ymax= ymin+nSzY_;
     253  double ymin= -y0_-0.5;
     254  double ymax= ymin+nSzY_;
    255255  if((y >  ymax) || (y < ymin)) return(-1);
    256   float ycurrent= ymin;
     256  double ycurrent= ymin;
    257257  for(j = 0; j < nSzY_; j++ )
    258258    {
     
    265265//++
    266266template<class T>
    267 void LocalMap<T>::PixThetaPhi(int_4 k, float& theta, float& phi) const
     267void LocalMap<T>::PixThetaPhi(int k,double& theta,double& phi) const
    268268//
    269269//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
     
    279279  Getij(k,i,j);
    280280
    281   float X= float(i-x0_);
    282   float Y= float(j-y0_);
     281  double X= double(i-x0_);
     282  double Y= double(j-y0_);
    283283  // situation de ce pixel dans le plan de reference
    284   float x= X*cos_angle_-Y*sin_angle_;
    285   float y= X*sin_angle_+Y* cos_angle_;
     284  double x= X*cos_angle_-Y*sin_angle_;
     285  double y= X*sin_angle_+Y* cos_angle_;
    286286  // projection sur la sphere
    287287  PixProjToAngle(x, y, theta, phi);
     
    292292//++
    293293template<class T>
    294 r_8 LocalMap<T>::PixSolAngle(int_4 k) const
     294double LocalMap<T>::PixSolAngle(int k) const
    295295//
    296296// Pixel Solid angle  (steradians)
     
    301301  int i,j;
    302302  Getij(k,i,j);
    303   float X= float(i-x0_);
    304   float Y= float(j-y0_);
    305   float XR= X+float(i)*0.5;
    306   float XL= X-float(i)*0.5;
    307   float YU= Y+float(j)*0.5;
    308   float YL= Y-float(j)*0.5;
     303  double X= double(i-x0_);
     304  double Y= double(j-y0_);
     305  double XR= X+double(i)*0.5;
     306  double XL= X-double(i)*0.5;
     307  double YU= Y+double(j)*0.5;
     308  double YL= Y-double(j)*0.5;
     309
    309310  // situation  dans le plan de reference
    310   float x0= XL*cos_angle_-YL*sin_angle_;
    311   float y0= XL*sin_angle_+YL*cos_angle_;
    312   float xa= XR*cos_angle_-YL*sin_angle_;
    313   float ya= XR*sin_angle_+YL*cos_angle_;
    314   float xb= XL*cos_angle_-YU*sin_angle_;
    315   float yb= XL*sin_angle_+YU*cos_angle_;
     311  double x0= XL*cos_angle_-YL*sin_angle_;
     312  double y0= XL*sin_angle_+YL*cos_angle_;
     313  double xa= XR*cos_angle_-YL*sin_angle_;
     314  double ya= XR*sin_angle_+YL*cos_angle_;
     315  double xb= XL*cos_angle_-YU*sin_angle_;
     316  double yb= XL*sin_angle_+YU*cos_angle_;
     317
    316318  // projection sur la sphere
    317   float tet0,phi0,teta,phia,tetb,phib;
    318   PixProjToAngle(x0, y0, tet0, phi0);
    319   PixProjToAngle(xa, ya, teta, phia);
    320   PixProjToAngle(xb, yb, tetb, phib);
     319  double thet0,phi0,theta,phia,thetb,phib;
     320  PixProjToAngle(x0, y0, thet0, phi0);
     321  PixProjToAngle(xa, ya, theta, phia);
     322  PixProjToAngle(xb, yb, thetb, phib);
     323
    321324  // angle solide
    322   float sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
    323   return r_8(sol);
    324 }
    325 
    326 //++
    327 template<class T>
    328 void  LocalMap<T>::SetOrigin(float theta0, float phi0, float angle)
     325  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
     326  return sol;
     327}
     328
     329//++
     330template<class T>
     331void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
    329332//
    330333// Fixe la repere de reference ( angles en degres)
    331334//--
    332335{
    333   theta0_= (double)theta0;
    334   phi0_  = (double)phi0;
    335   angle_ = (double)angle;
     336  theta0_= theta0;
     337  phi0_  = phi0;
     338  angle_ = angle;
    336339  x0_= nSzX_/2;
    337340  y0_= nSzY_/2;
    338   cos_angle_= cosdf(angle);
    339   sin_angle_= sindf(angle);
     341  cos_angle_= cos(angle*Pi/180.);
     342  sin_angle_= sin(angle*Pi/180.);
    340343  originFlag_= true;
    341344  cout << " LocalMap:: set origin 1 done" << endl;
     
    344347//++
    345348template<class T>
    346 void  LocalMap<T>::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
     349void LocalMap<T>::SetOrigin(double theta0,double phi0,int x0,int y0,double angle)
    347350//
    348351// Fixe le repere de reference (angles en degres) 
    349352//--
    350353{
    351   theta0_= (double)theta0;
    352   phi0_  = (double)phi0;
    353   angle_ = (double)angle;
     354  theta0_= theta0;
     355  phi0_  = phi0;
     356  angle_ = angle;
    354357  x0_= x0;
    355358  y0_= y0;
    356   cos_angle_= cosdf(angle);
    357   sin_angle_= sindf(angle);
     359  cos_angle_= cos(angle*Pi/180.);
     360  sin_angle_= sin(angle*Pi/180.);
    358361  originFlag_= true;
    359362  cout << " LocalMap:: set origin 2 done" << endl;
     
    362365//++
    363366template<class T>
    364 void LocalMap<T>::SetSize(float angleX, float angleY)
     367void LocalMap<T>::SetSize(double angleX,double angleY)
    365368//
    366369// Fixe l'extension de la carte (angles en degres)
    367370//--
    368371{
    369   angleX_= (double)angleX;
    370   angleY_= (double)angleY;
    371 
    372   //tgAngleX_= T(tan(angleX*Pi/360.));
    373   //tgAngleY_= T(tan(angleY*Pi/360.));
    374 
     372  angleX_= angleX;
     373  angleY_= angleY;
    375374
    376375  // tangente de la moitie de l'ouverture angulaire totale
    377   tgAngleX_= tand(0.5*angleX_);
    378   tgAngleY_= tand(0.5*angleY_);
     376  tgAngleX_= tan(0.5*angleX_*Pi/180.);
     377  tgAngleY_= tan(0.5*angleY_*Pi/180.);
    379378
    380379  extensFlag_= true;
     
    391390  for(int m = 0; m < nPix_; m++)
    392391    {
    393       float theta,phi;
     392      double theta,phi;
    394393      PixThetaPhi(m,theta,phi);
    395394      sphere(theta,phi)= pixels_(m);
     
    400399//++
    401400template<class T>
    402 void LocalMap<T>::Getij(int k, int& i, int& j) const
     401void LocalMap<T>::Getij(int k,int& i,int& j) const
    403402//
    404403//--
     
    411410//++
    412411template<class T>
    413 void  LocalMap<T>::ReferenceToUser(float &theta, float &phi) const
     412void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
    414413//
    415414// --     
    416415{
    417   if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi)
    418     {
    419       //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
    420       //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     416  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     417    {
    421418      throw "LocalMap::ReferenceToUser (theta,phi) out of range";
    422419    }
    423420
    424   //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
    425   theta= (r_4)theta0_*Pi/180.+theta-(r_4)Pi*0.5;
     421  theta= theta0_*Pi/180.+theta-Pi*0.5;
    426422  if(theta < 0.)
    427423    {
    428424      theta= -theta;
    429       phi += (r_4)Pi;
     425      phi += Pi;
    430426    }
    431427  else
    432428    {
    433       if(theta > (r_4)Pi)
     429      if(theta > Pi)
    434430        {
    435           theta= 2.*(r_4)Pi-theta;
    436           phi += (r_4)Pi;
     431          theta= 2.*Pi-theta;
     432          phi += Pi;
    437433        }
    438434    }
    439435
    440   phi= (r_4)phi0_*Pi/180.+phi;
    441   while(phi >= 2.*(r_4)Pi) phi-=2.*(r_4)Pi;
    442 
    443   if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi)
     436  phi= phi0_*Pi/180.+phi;
     437  while(phi >= 2.*Pi) phi-= 2.*Pi;
     438
     439  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    444440    {
    445441      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
     
    447443      exit(0);
    448444    }
    449   //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
    450 }
    451 
    452 //++
    453 template<class T>
    454 void  LocalMap<T>::UserToReference(float &theta, float &phi) const
     445}
     446
     447//++
     448template<class T>
     449void  LocalMap<T>::UserToReference(double& theta,double& phi) const
    455450//
    456451// --     
    457452{
    458   if(theta > (r_4)Pi || theta < 0.  || phi<0. || phi >= 2*(r_4)Pi )
    459     {
    460       cout << " LocalMap::UserToReference : exceptions a mettre en place" <<endl;
    461       //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     453  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
     454    {
     455      cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;
     456      // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    462457      throw "LocalMap::UserToReference (theta,phi) out of range";
    463458    }
    464459
    465   float phi1=phi-(r_4)phi0_*Pi/180.;
    466   if(phi1 < 0.) phi1+=2.*(r_4)Pi;
    467 
    468   float theta1= theta-(r_4)theta0_*Pi/180.+(r_4)Pi*0.5;
     460  double phi1= phi-phi0_*Pi/180.;
     461  if(phi1 < 0.) phi1+= 2.*Pi;
     462
     463  double theta1= theta-theta0_*Pi/180.+Pi*0.5;
    469464  if(theta1 < 0.)
    470465    {
    471466      theta= -theta1;
    472       phi1+= (r_4)Pi;
     467      phi1+= Pi;
    473468    }
    474469  else
    475470    {
    476       if(theta1 > (r_4)Pi)
     471      if(theta1 > Pi)
    477472        {
    478           theta= 2.*(r_4)Pi-theta1;
    479           phi1+= (r_4)Pi;
     473          theta= 2.*Pi-theta1;
     474          phi1+= Pi;
    480475        }
    481476    }
    482477
    483   while(phi1 >= 2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
     478  while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
    484479  phi= phi1;
    485   if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi )
     480  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)
    486481    {
    487482      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
     
    493488//++
    494489template<class T>
    495 void LocalMap<T>::PixProjToAngle(float x, float y, float& theta, float& phi) const
    496 //
    497 // (x,y) representent les coordonnees en unites de pixels d'un point DANS  LE PLAN DE REFERENCE.
     490void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const
     491//
     492// (x,y) representent les coordonnees en unites de pixels d'un point DANS LE PLAN DE REFERENCE.
    498493// On recupere (theta,phi) par rapport au repere "absolu" theta=pi/2 et phi=0.
    499494//--
    500495{
    501   theta= (r_4)Pi*0.5-atan(2*y*tgAngleY_/(float)nSzY_);
    502   phi= atan2(2*x*tgAngleX_,(float)nSzX_);
     496  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
     497  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
    503498  if(phi < 0.) phi += DeuxPi;
    504499}
     
    506501//++
    507502template<class T>
    508 void LocalMap<T>::AngleProjToPix(float theta, float phi, float& x, float& y) const
     503void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
    509504//
    510505// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere
     
    512507//--
    513508{
    514   if(phi >= (r_4)Pi) phi-= DeuxPi;
     509  if(phi >= Pi) phi-= DeuxPi;
    515510  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
    516511  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ? 
     
    606601    }
    607602
    608   int_4 nSzX;
     603  int nSzX;
    609604  is.GetI4(nSzX);
    610605  dobj->setSize_x(nSzX);
    611606
    612   int_4 nSzY;
     607  int nSzY;
    613608  is.GetI4(nSzY);
    614609  dobj->setSize_y(nSzY);
    615610
    616   int_4 nPix;
     611  int nPix;
    617612  is.GetI4(nPix);
    618613  dobj->setNbPixels(nPix);
     
    624619    {
    625620      cout<<" ReadSelf:: local mapping"<<endl;
    626       int_4 x0, y0;
    627       float theta, phi, angle;
     621      int x0, y0;
     622      double theta, phi, angle;
    628623      is.GetI4(x0);
    629624      is.GetI4(y0);
    630       is.GetR4(theta);
    631       is.GetR4(phi);
    632       is.GetR4(angle);
     625      is.GetR8(theta);
     626      is.GetR8(phi);
     627      is.GetR8(angle);
    633628      dobj->SetOrigin(theta, phi, x0, y0, angle);
    634629
    635       float angleX, angleY;
    636       is.GetR4(angleX);
    637       is.GetR4(angleY);
     630      double angleX, angleY;
     631      is.GetR8(angleX);
     632      is.GetR8(angleY);
    638633      dobj->SetSize(angleX, angleY);
    639634    }
     
    657652
    658653  char strg[256];
    659   int_4 nSzX= dobj->Size_x();
    660   int_4 nSzY= dobj->Size_y();
    661   int_4 nPix= dobj->NbPixels();
     654  int nSzX= dobj->Size_x();
     655  int nSzY= dobj->Size_y();
     656  int nPix= dobj->NbPixels();
    662657 
    663658  if(dobj->ptrInfo())
     
    681676      string ss("local mapping is done");
    682677      os.PutStr(ss);
    683       int_4 x0, y0;
    684       float theta, phi, angle;
     678      int x0, y0;
     679      double theta, phi, angle;
    685680      dobj->Origin(theta, phi, x0, y0, angle);
    686681      os.PutI4(x0);
    687682      os.PutI4(y0);
    688       os.PutR4(theta);
    689       os.PutR4(phi);
    690       os.PutR4(angle);
    691 
    692       float angleX, angleY;
     683      os.PutR8(theta);
     684      os.PutR8(phi);
     685      os.PutR8(angle);
     686
     687      double angleX, angleY;
    693688      dobj->Aperture(angleX, angleY);
    694       os.PutR4(angleX);
    695       os.PutR4(angleY);
     689      os.PutR8(angleX);
     690      os.PutR8(angleY);
    696691    }
    697692  else
Note: See TracChangeset for help on using the changeset viewer.