source: Sophya/trunk/Poubelle/DPC:FitsIOServer/Samba/scan.cc@ 3488

#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "scan.h"
#include <complex>
#include "piocmplx.h"

// valeurs de Pi, 2*Pi, etc                                              
#include "nbmath.h"


//
//++
// Class        Scan
//
// include      scan.h  dvlist.h math.h nbmath.h 
//
//    Storage and treatment of data for a scanning of a part of the sky
//    with a set of given values for parameters (see constructor)
//--
//++
Scan::Scan(float Ouv,float* Omega,float Fech,float T, 
           float t0=0., float phi0=0.)

//  
//    * Ouv      = aperture angle (rad)
//    * Omega[3] = direction of rotation axis of the satellite (teta,phi)
//    and rotation velocity (rad/s)
//    * Fech     = sampling frequency  (Hz)
//    * T        = total time of data acquistion (s)
//    * t0       = starting time (s)
//    * phi0     = offset of antenna (rad)
//--
{
  mInfo_=NULL;
  // On memorise les arguments d'appel
  Ouverture_= Ouv;
  OmegaTeta_= Omega[0];
  //
  // je ne comprends pas ce qui se passe ci-apres (GLM)
  if( (Omega[0]== 0.0) || (Omega[0]==(float)Pi) ) 
    OmegaPhi_ = Omega[1];
  else
    OmegaPhi_ = Omega[1]+(float)Pi/2.;
  OmegaRad_ = Omega[2];
  FrequenceEch_  = Fech;
  TempsFinal_    = T;
  TempsInitial_  = t0;
  PhiZero_ = phi0;

  float aux= 0.0;
  // Nombre total de points
  aux= (TempsFinal_-TempsInitial_)*FrequenceEch_;
  NmaxPts_= (int_4)(aux);

  // Nombre total de tours
  aux= OmegaRad_*(TempsFinal_-TempsInitial_)/(2.*Pi);
  NmaxTrs_= (int_4)(aux);

  // Nombre de points par tour
  aux= 2.*Pi*FrequenceEch_/OmegaRad_;
  NPts1Tr_= (int_4)(aux);

  // Creation et initialisation du vecteur des mesures aux points
  //  sPix_ = new r_8[NmaxPts_];
  sPix_.ReSize(NmaxPts_);
  sPix_.Reset();
  //  for( int_4 i=0; i<NmaxPts_; i++ )  sPix_[i]= 0.;
    for( int_4 i=0; i<NmaxPts_; i++ )  sPix_(i)= 0.;

  // matrice de passage du systeme lie au satellite au systeme fixe

  Rota_[0]=  cos((double)OmegaPhi_);
  Rota_[1]= -sin((double)OmegaPhi_)*cos((double)OmegaTeta_);
  Rota_[2]=  sin((double)OmegaPhi_)*sin((double)OmegaTeta_);
  Rota_[3]=  sin((double)OmegaPhi_);
  Rota_[4]=  cos((double)OmegaPhi_)*cos((double)OmegaTeta_);
  Rota_[5]= -cos((double)OmegaPhi_)*sin((double)OmegaTeta_);
  Rota_[6]=  0.0;
  Rota_[7]=  sin((double)OmegaTeta_);
  Rota_[8]=  cos((double)OmegaTeta_);


  // printf("%f  %f  %f \n",Rota_[0],Rota_[1],Rota_[2]);
  // printf("%f  %f  %f \n",Rota_[3],Rota_[4],Rota_[5]);
  // printf("%f  %f  %f \n",Rota_[6],Rota_[7],Rota_[8]);

}


//++
Scan::Scan(const Scan& s)

//    copy  constructor
//--
{
  NmaxPts_      = s.NmaxPts_;
  Ouverture_    = s.Ouverture_;
  OmegaTeta_    = s.OmegaTeta_;
  OmegaPhi_     = s.OmegaPhi_;
  OmegaRad_     = s.OmegaRad_;
  FrequenceEch_ = s.FrequenceEch_;
  TempsFinal_   = s.TempsFinal_;
  TempsInitial_ = s.TempsInitial_;
  PhiZero_      = s.PhiZero_;
  sPix_=s.sPix_;
  for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k];
}
//++
// Titre        Destructor
//--
//++
Scan::~Scan()

//     
//--
{
  // delete[] sPix_ ;
}


//++
// Titre        Public Methods
//--

//++
int_4 Scan::NbPoints() const

//    Return the number of points in the scan
//--
{
  return(NmaxPts_);
}

/* --Methode-- */
//++
int_4 Scan::NbTours() const

//    Return total nomber of turns
//--
{
  return(NmaxTrs_);
}

/* --Methode-- */
//++
int_4 Scan::NbPts1Tr() const

//    Return  nomber of points for 1 turn
//--
{
  return(NPts1Tr_);
}

/* --Methode-- */
//++
int_4 Scan::ValueIndex(float t) const

//    Return index of pixel associated to time t
//--
{
  int_4 k;
  float eps= 1.0E-06;

  // verification si t est dans [TempsInitial_,TempsFinal_]
  if( (t< TempsInitial_) || (t> TempsFinal_) ) {
    printf("\n ValueIndex:: probleme sur le temps t= %f",t);
    return(-1);
  }

  k= (int_4)((t-TempsInitial_)*FrequenceEch_+eps);
  if ( (k< 0) || (k >= NmaxPts_) ) {
    printf("\n ValueIndex:: probleme sur l'indice du point k= %d",k); 
    return(-1);
  }
  return(k);
}

/* --Methode-- */
//++
void Scan::Direction(float t, float&  teta , float&  phi)

//    Return  (teta,phi) coordinate of pixel related to time t
//--
{
  r_8 alfa,xs,ys,zs,x,y,z;

  // coordonnees dans le systeme du satellite
  alfa= OmegaRad_*(t-TempsInitial_)+PhiZero_;
  xs = sin((double)Ouverture_)*cos(alfa);
  ys = sin((double)Ouverture_)*sin(alfa);
  zs = cos((double)Ouverture_);

  // coordonnees dans le systeme fixe
  x = Rota_[0]*xs+Rota_[1]*ys+Rota_[2]*zs;
  y = Rota_[3]*xs+Rota_[4]*ys+Rota_[5]*zs;
  z = Rota_[6]*xs+Rota_[7]*ys+Rota_[8]*zs;

  // angles teta,phi
  teta = acos(z);
  phi  = atan2(y,x);
  if( phi< 0. ) phi= DeuxPi+phi;
}

/* --Methode-- */
//++
void Scan::DirectionIndex(int_4 k, float&  teta, float&  phi)

//    Return  (teta,phi) coordinates of  pixel with index k
//--
{
  float t;

  //recupere le temps associe a l'indice k du pixel
  t= TempsInitial_+(float)k/FrequenceEch_;

  // angles teta,phi 
  Direction(t, teta, phi);
}

static r_8 dummy_pixel = 0;
/* --Methode-- */
//++
r_8 const& Scan::PixelValue(int_4 k) const

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >=  NmaxPts_) ) {
    printf("\n ValueIndex::indice du pixel errone k= %d",k); 
    return(dummy_pixel);
  }
  //  return(sPix_[k]);
  //return(sPix_(k));
  return *(sPix_.Data()+k);
}
//++
r_8 & Scan::PixelValue(int_4 k) 

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >=  NmaxPts_) ) {
    printf("\n ValueIndex::indice du pixel errone k= %d",k); 
    return(dummy_pixel);
  }
  //  return(sPix_[k]);
  //return(sPix_(k));
  return *(sPix_.Data()+k);
}

//void Scan::Clear() {
//  if (sPix_)  delete[] sPix_;
//}         
void Scan::InitNull() {
  //  sPix_=NULL;
  sPix_.Reset();
  mInfo_=NULL;
}         
Scan&  Scan::operator = (const Scan& s)
{
  NmaxPts_      = s.NmaxPts_;
  Ouverture_    = s.Ouverture_;
  OmegaTeta_    = s.OmegaTeta_;
  OmegaPhi_     = s.OmegaPhi_;
  OmegaRad_     = s.OmegaRad_;
  FrequenceEch_ = s.FrequenceEch_;
  TempsFinal_   = s.TempsFinal_;
  TempsInitial_ = s.TempsInitial_;
  PhiZero_      = s.PhiZero_;
  sPix_=s.sPix_;
  for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k];
  return *this;
}


//++
// Titre        Operators
//--
//++
// Links        double&        operator()(int_4 k)
//--
//++
//    fill and/or return value of the pixel with index k
//--





//++
// Titre        class FIO_Scan
//    Delegated objects for persitance management
//--

//*******************************************************************
// class FIO_Scan
//  Les objets delegues pour la gestion de persistance 
//*******************************************************************

//++
FIO_Scan::FIO_Scan()
//
//--
{
  dobj= new Scan;
  ownobj= true;
}
//++
FIO_Scan::FIO_Scan(string const& filename)
//
//--
{
  dobj= new Scan;
  ownobj= true;
  Read(filename);
}
//++
FIO_Scan::FIO_Scan(const Scan& obj)
//
//--
{
  dobj= new Scan(obj);
  ownobj= true;
}
FIO_Scan::FIO_Scan(Scan* obj)
{
  dobj= obj;
  ownobj= false;
}
//++
FIO_Scan::~FIO_Scan()
//
//--
{
  if (ownobj && dobj) delete dobj;
}
//++
AnyDataObj* FIO_Scan::DataObj()
//
//--
{
  return(dobj);
}

//++
void FIO_Scan::ReadSelf(PInPersist& is)
//
//--
{

  int NmaxPts, NmaxTrs, NPts1Tr;
  r_4 Ouverture, OmegaTeta, OmegaPhi, OmegaRad, FrequenceEch; 
  r_4 TempsFinal,  TempsInitial, PhiZero; 
  r_8 Rota[9];

  if(dobj == NULL) 
    {
      dobj= new Scan;
    }

  // Pour savoir s'il y avait un DVList Info associe
  char strg[256];
  is.GetLine(strg, 255);
  bool hadinfo= false;
  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
  if(hadinfo) 
    {    // Lecture eventuelle du DVList Info
      is >> dobj->Info(); 
    }


  is.GetI4(NmaxPts);
  is.GetI4(NmaxTrs);
  is.GetI4(NPts1Tr);
  dobj->SetIntParams(NmaxPts,NmaxTrs,NPts1Tr);

  is.GetR4(Ouverture);
  is.GetR4(OmegaTeta);
  is.GetR4(OmegaPhi);
  is.GetR4(OmegaRad);
  is.GetR4(FrequenceEch);
  is.GetR4(TempsFinal);
  is.GetR4(TempsInitial);
  is.GetR4(PhiZero);
  is.GetR8s(Rota, 9);
  dobj->SetFloatParams(Ouverture,OmegaTeta,OmegaPhi,OmegaRad,
                       FrequenceEch,TempsFinal,TempsInitial,PhiZero,Rota);
  r_8* sPix=new r_8[NmaxPts];
  PIOSReadArray(is, sPix, NmaxPts);
  dobj->setDataBlock(sPix, NmaxPts);
  delete [] sPix;
}
//++
void FIO_Scan::WriteSelf(POutPersist& os) const
//
//--
{
  r_4 Ouverture, OmegaTeta, OmegaPhi, OmegaRad, FrequenceEch; 
  r_4 TempsFinal,  TempsInitial, PhiZero; 
  r_8 Rota[9];
  char strg[256];

  if(dobj == NULL) 
    {
      cout << " FIO_Scan::WriteSelf:: dobj= null " << endl;
      return;
    }
  dobj->GetFloatParams(Ouverture,OmegaTeta,OmegaPhi,OmegaRad,
                       FrequenceEch,TempsFinal,TempsInitial,PhiZero,Rota);
  if (dobj->ptrInfo()) 
    {sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f HasInfo", 
             (float)OmegaTeta, (float)OmegaPhi, (float)OmegaRad);
    os.PutLine(strg);
    os << dobj->Info();
    }
  else 
    {
      sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f ", 
              (float)OmegaTeta, (float)OmegaPhi, (float)OmegaRad); 
      os.PutLine(strg);
    }

  os.PutI4(dobj->NbPoints());
  os.PutI4(dobj->NbTours());
  os.PutI4(dobj->NbPts1Tr());

  os.PutR4(Ouverture);
  os.PutR4(OmegaTeta);
  os.PutR4(OmegaPhi);
  os.PutR4(OmegaRad);
  os.PutR4(FrequenceEch);
  os.PutR4(TempsFinal);
  os.PutR4(TempsInitial);
  os.PutR4(PhiZero);
  os.PutR8s(Rota, 9);

  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), dobj->NbPoints());

}