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

#include "scan.h"

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


//
//++
// Class	Scan
//
// include      scan.h  nbmath.h 
//
//    Cette classe permet de stocker et traiter les données pour un balayage 
//    d'une portion de ciel pour un jeu de valeurs fixé des paramètres précisé
//    dans le constructeur.    
//--
//++
//
// Links	Parents
//
//    
//-- 
//++
//
// Links	Descendants
//
//     
//-- 
//++
// Titre	Constructeurs
//--
//++
Scan::Scan(float Ouv,float* Omega,float Fech,float T, 
	   float t0=0., float phi0=0.)

//  
//    * Ouv      = angle d'ouverture (rad)
//    * Omega[3] = direction de l'axe de rotation du satellite (teta,phi)
//    et vitesse de rotation (rad/s)
//    * Fech     = fréquence d'échantillonnage (Hz)
//    * T        = temps total de prise des données (s)
//    * t0       = instant de départ (s)
//    * phi0     = offset antenne (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_];
  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)

//    Constructeur de copie
//--
{
  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_=new r_8[ NmaxPts_];
  for (int_4 k=0; k<NmaxPts_; k++)  sPix_[k]=s.sPix_[k];
  for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k];
}
//++
// Titre	Destructeur
//--
//++
Scan::~Scan()

//     
//--
{
  delete[] sPix_ ;
}
//++
void Scan::WriteSelf(POutPersist& s) const

//    créer un fichier persistant
//--
{
  char strg[256];
  if (mInfo_) {sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f HasInfo", (float)OmegaTeta_, (float)OmegaPhi_, (float)OmegaRad_);
  }
  else {sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f ", (float)OmegaTeta_, (float)OmegaPhi_, (float)OmegaRad_); 
  }
  s.PutLine(strg);
  if (mInfo_)  s << (*mInfo_);
   s.PutI4( NmaxPts_);
  s.PutI4(NmaxTrs_);
  s.PutI4( NPts1Tr_);
  s.PutR4(Ouverture_);
  s.PutR4(OmegaTeta_);
  s.PutR4(OmegaPhi_);
  s.PutR4(OmegaRad_);
  s.PutR4(FrequenceEch_);
  s.PutR4(TempsFinal_);
  s.PutR4( TempsInitial_);
  s.PutR4( PhiZero_);
  s.PutR8s(sPix_,NmaxPts_);
  s.PutR8s(Rota_, 9);
  return;
}


//++
void Scan::ReadSelf(PInPersist& s) 

//    relire un fichier persistant
//--
{
  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;
// Pour savoir s'il y avait un DVList Info associe
  if (hadinfo) {    // Lecture eventuelle du DVList Info
    if (mInfo_ == NULL)  mInfo_ = new DVList;
    s >> (*mInfo_);
  }
  s.GetI4( NmaxPts_);
  s.GetI4(NmaxTrs_);
  s.GetI4( NPts1Tr_);
  s.GetR4(Ouverture_);
  s.GetR4(OmegaTeta_);
  s.GetR4(OmegaPhi_);
  s.GetR4(OmegaRad_);
  s.GetR4(FrequenceEch_);
  s.GetR4(TempsFinal_);
  s.GetR4( TempsInitial_);
  s.GetR4( PhiZero_);
  sPix_=new r_8[NmaxPts_];
  s.GetR8s(sPix_,NmaxPts_);
  s.GetR8s(Rota_, 9);
}


//++
// Titre	Méthodes
//--

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

//    Retourne le nombre de Points du scan
//--
{
  return(NmaxPts_);
}

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

//    Retourne le nombre total de tours
//--
{
  return(NmaxTrs_);
}

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

//    Retourne le nombre de points pour 1 tour
//--
{
  return(NPts1Tr_);
}

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

//    Retourne l'indice du pixel associé au temps 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)

//    Retourne les coordonnées (teta,phi) du pixel lié au temps 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)

//    Retourne les coordonnées (teta,phi) du pixel d'indice 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& 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]);
}

void Scan::Clear() {
  if (sPix_)  delete[] sPix_;
}         
void Scan::InitNull() {
  sPix_=NULL;
  mInfo_=NULL;
}         


//++
// Titre	Opérateurs
//--
//++
// Links	double&        operator()(int_4 k)
//--
//++
//    Surcharge de la parenthese a un indice entier : remplit ou/et renvoie 
//    la valeur du contenu du pixel d'indice k                               
//--