source: Sophya/trunk/SophyaLib/TArray/utilarr.cc@ 906

//     Utility classes for template numerical arrays
//                     R. Ansari, C.Magneville   03/2000

#include "machdefs.h"
#include "utilarr.h"
#include "srandgen.h"

// Classe utilitaires

//////////////////////////////////////////////////////////
//! Constructor
/*!
  \param typ : generator type
  \param m : mean parameter of the generator (if needed)
  \param s : sigma parameter of the generator (if needed)
 */
RandomSequence::RandomSequence(int typ, double m, double s)
{
  typ_ = (typ == Flat) ? Flat : Gaussian;
  mean_ = m;
  sig_ = s;
}

//! Return random sequence values.
double RandomSequence::Rand()
{
  if (typ_ == Flat) 
    return(drandpm1()*sig_ + mean_);
  else return(GauRnd(mean_, sig_));
}


//////////////////////////////////////////////////////////
//! Constructor
/*!
  \param start : start value
  \param step : step value
  \param f : pointer to the sequence function

  See \ref SequenceOperat "operator()"
 */
Sequence::Sequence(double start, double step, Arr_DoubleFunctionOfX f)
  : rseq_(RandomSequence::Gaussian)
{
  start_ = start;
  step_ = step;
  myf_ = f;
  fgrseq_ = false;
}

//! Constructor
/*!
  \param rseq : RandomSequence

  See \ref SequenceOperat "operator()"
 */
Sequence::Sequence(RandomSequence rseq)
{
  start_ = 0.;
  step_ = 1.;
  myf_ = NULL;  
  rseq_ = rseq;
  fgrseq_ = true;
}

//! Get the \b k th value
/*!
  \param k : index of the value
  \anchor SequenceOperat

  If the constructor was done with RandomSequence, return a RandomSequence
  and \b k doesn't matter.

  If the constructor has a NULL Arr_DoubleFunctionOfX, return start+k*step

  If the constructor has a not NULL Arr_DoubleFunctionOfX, return f(start+k*step)
  \return the \b k th value
 */
double Sequence::operator () (uint_4 k)
{
  if (fgrseq_) return(rseq_.Rand());
  else {
    double x = start_+(double)k*step_;
    if (myf_)  return(myf_(x));
    else return x;
  }
}

//////////////////////////////////////////////////////////
//! Constructor
/*!
  Define a range of indexes
  \param start : start index
  \param end : start end
  \param size : size
  \param step : step

  \warning If \b end \> \b start, \b size is computed automatically
  \warning If not \b size is fixed and \b end recomputed
 */
Range::Range(uint_4 start, uint_4 end, uint_4 size, uint_4 step)
{
  start_ = start;
  step_ = (step > 0) ? step : 1;
  if (end > start) {  // Taille calcule automatiquement 
    end_ = end;
    if (step_ > ((end_-start_)+1))  size_ = 1;
    else size_ = ((end-start)+1)/step_;
  }
  else {     // Taille fixee
    size_ = size;
    end_ = start_+size_*step_;
  }
}

/*
Range & Range::operator = (uint_4 start)
{
  start_ = start;
  size_ = 1;
  step_ = 1;
  return (*this);
}
*/


//////////////////////////////////////////////////////////
//! Constructor of a (n,n) diagonal matrix with value diag on the diagonal
IdentityMatrix::IdentityMatrix(double diag, uint_4 n)
{
  size_ = n;
  diag_ = diag;
}