Index: trunk/SophyaLib/TArray/basarr.cc
===================================================================
--- trunk/SophyaLib/TArray/basarr.cc	(revision 1493)
+++ trunk/SophyaLib/TArray/basarr.cc	(revision 1517)
@@ -290,5 +290,5 @@
   in this case.
 */
-bool BaseArray::CompareSizes(const BaseArray& a, bool& smo)
+bool BaseArray::CompareSizes(const BaseArray& a, bool& smo) const
 {
   if (ndim_ != a.ndim_)  return(false);
@@ -454,5 +454,5 @@
 //! return various parameters for double loop operations on two arrays.
 void BaseArray::GetOpeParams(const BaseArray& a, bool smo, int_4& ax, int_4& axa, sa_size_t& step,
-			     sa_size_t& stepa, sa_size_t& gpas, sa_size_t& naxa)
+			     sa_size_t& stepa, sa_size_t& gpas, sa_size_t& naxa) const
 {
   if (smo) { // Same memory organisation
Index: trunk/SophyaLib/TArray/basarr.h
===================================================================
--- trunk/SophyaLib/TArray/basarr.h	(revision 1493)
+++ trunk/SophyaLib/TArray/basarr.h	(revision 1517)
@@ -62,5 +62,5 @@
 
   // Returns true if ndim and sizes are equal
-  virtual bool CompareSizes(const BaseArray& a, bool& smo);
+  virtual bool CompareSizes(const BaseArray& a, bool& smo) const;
 
   // Compacts \b size=1 array dimensions
@@ -143,5 +143,5 @@
   // Pour recuperer pas et numero d'axe pour operations sur deux arrays
   void   GetOpeParams(const BaseArray& a, bool smo, int_4& ax, int_4& axa, sa_size_t& step,
-		      sa_size_t& stepa, sa_size_t& gpas, sa_size_t& naxa);
+		      sa_size_t& stepa, sa_size_t& gpas, sa_size_t& naxa) const;
   // Impression, I/O, ...
   void           Show(ostream& os, bool si=false) const;
Index: trunk/SophyaLib/TArray/matharr.cc
===================================================================
--- trunk/SophyaLib/TArray/matharr.cc	(revision 1493)
+++ trunk/SophyaLib/TArray/matharr.cc	(revision 1517)
@@ -152,8 +152,334 @@
 }
 
+
+//-------------------------------------------------------------------------------
+//      Definition utilitaire d'application de fonction 
+inline complex<r_8> ApplyComplexDoubleFunction(complex<r_8> z, 
+					       Arr_ComplexDoubleFunctionOfX f)
+{
+  return(f(z));
+}
+
+inline complex<r_4> ApplyComplexDoubleFunction(complex<r_4> z, 
+					       Arr_ComplexDoubleFunctionOfX f)
+{
+  complex<r_8> zd((r_8)z.real(), (r_8)z.imag());
+  zd = f(zd);
+  complex<r_8> zr((r_4)zd.real(), (r_4)zd.imag());
+  return(zr);
+}
+
+//-------------------------------------------------------------------------------
+
+/*!
+  \class SOPHYA::ComplexMathArray
+  \ingroup TArray
+  Class for simple mathematical operation on arrays 
+  \warning Instanciated only for \b real and \b double (r_4, r_8) complex arrays
+*/
+
+//! Apply Function In Place (complex arrays)
+/*!
+  \param a : complex array to be replaced in place
+  \param f : function for replacement
+  \return Return an array \b a filled with function f(a(i,j))
+*/
+template <class T>
+TArray< complex<T> >& ComplexMathArray<T>::ApplyFunctionInPlace(TArray< complex<T> > & a, Arr_ComplexDoubleFunctionOfX f)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::ApplyFunctionInPlace(TArray< complex<T> > & a..) Not Allocated Array a !");
+  complex<T> * pe;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    for(k=0; k<maxx; k+=step )  pe[k] = ApplyComplexDoubleFunction(pe[k],f); 
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  pe[k] = ApplyComplexDoubleFunction(pe[k],f);
+    }
+  }
+  return(a);
+}
+
+
+
+//! Apply Function (complex arrays)
+/*!
+  \param a : argument array of the function
+  \param f : function for replacement
+  \return Return a new array filled with function f(a(i,j))
+*/
+template <class T>
+TArray< complex<T> > ComplexMathArray<T>::ApplyFunction(TArray< complex<T> > const & a, Arr_ComplexDoubleFunctionOfX f)
+{
+  TArray< complex<T> > ra;
+  ra = a;
+  ApplyFunctionInPlace(ra, f);
+  return(ra);
+}
+
+//! Create a complex array, from a real and an imaginary arrays
+/*!
+  \param p_real : array containing the real part of the complex output array
+  \param p_imag : array containing the imaginary part of the complex output array
+  \return Return a new complex array build from \b p_real and \b p_imag
+*/
+template <class T>
+TArray< complex<T> > ComplexMathArray<T>::FillFrom(TArray<T> const & p_real,
+						   TArray<T> const & p_imag)
+{
+  if (p_real.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray<T>::FillFrom() - Not Allocated Array ! ");
+  bool smo;
+  if (!p_real.CompareSizes(p_imag, smo)) 
+    throw(SzMismatchError("ComplexMathArray<T>::FillFrom() SizeMismatch")) ;
+
+  TArray< complex<T> > ra;
+  ra.ReSize(p_real);
+
+  complex<T> * pe;
+  const T * per;
+  const T * pei;
+  sa_size_t j,k,ka;
+  if (smo && (p_real.AvgStep() > 0) && (p_imag.AvgStep() > 0))   {  // regularly spaced elements
+    sa_size_t step = p_real.AvgStep();
+    sa_size_t stepa = p_imag.AvgStep();
+    sa_size_t maxx = p_real.Size()*step;
+    per = p_real.Data();
+    pei = p_imag.Data();
+    pe = ra.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  
+      pe[k] = complex<T>(per[k], pei[ka]) ;
+  }
+  else {    // Non regular data spacing ...
+    int_4 ax,axa;
+    sa_size_t step, stepa;
+    sa_size_t gpas, naxa;
+    p_real.GetOpeParams(p_imag, smo, ax, axa, step, stepa, gpas, naxa);
+    for(j=0; j<naxa; j++)  {
+      per = p_real.Data()+p_real.Offset(ax,j);
+      pei = p_imag.Data()+p_imag.Offset(axa,j);
+      pe = ra.Data()+ra.Offset(ax,j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  
+	pe[k] = complex<T>(per[k], pei[ka]) ;
+    }
+  }
+  return(ra);
+}
+
+
+//! Returns the real part of the complex input array.
+/*!
+  \param a : input complex array
+  \return Return a new array filled with the real part of the input complex array elements
+*/
+
+template <class T>
+TArray<T> ComplexMathArray<T>::real(TArray< complex<T> > const & a)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::real(TArray< complex<T> >& a) Not Allocated Array a !");
+  TArray<T> ra;
+  ra.ReSize(a);
+
+  const complex<T> * pe;
+  T * po;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    po = ra.Data();
+    for(k=0; k<maxx; k+=step )  po[k] = pe[k].real();
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      po = ra.DataBlock().Begin()+ra.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  po[k] = pe[k].real();
+    }
+  }
+  return(ra);
+}
+
+//! Returns the imaginary part of the complex input array.
+/*!
+  \param a : input complex array
+  \return Return a new array filled with the imaginary part of the input complex array elements
+*/
+
+template <class T>
+TArray<T> ComplexMathArray<T>::imag(TArray< complex<T> > const & a)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::imag(TArray< complex<T> >& a) Not Allocated Array a !");
+  TArray<T> ra;
+  ra.ReSize(a);
+
+  const complex<T> * pe;
+  T * po;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    po = ra.Data();
+    for(k=0; k<maxx; k+=step )  po[k] = pe[k].imag();
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      po = ra.DataBlock().Begin()+ra.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  po[k] = pe[k].imag();
+    }
+  }
+  return(ra);
+}
+
+//! Returns the module squared of the complex input array.
+/*!
+  \param a : input complex array
+  \return Return a new array filled with the module squared of the input complex array elements 
+*/
+
+template <class T>
+TArray<T> ComplexMathArray<T>::module2(TArray< complex<T> > const & a)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::module2(TArray< complex<T> >& a) Not Allocated Array a !");
+  TArray<T> ra;
+  ra.ReSize(a);
+
+  const complex<T> * pe;
+  T * po;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    po = ra.Data();
+    for(k=0; k<maxx; k+=step )  
+      po[k] = (pe[k].real()*pe[k].real()+pe[k].imag()*pe[k].imag());
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      po = ra.DataBlock().Begin()+ra.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  
+	po[k] = (pe[k].real()*pe[k].real()+pe[k].imag()*pe[k].imag());
+    }
+  }
+  return(ra);
+}
+
+//! Returns the module of the complex input array.
+/*!
+  \param a : input complex array
+  \return Return a new array filled with the module of the input complex array elements 
+*/
+
+template <class T>
+TArray<T> ComplexMathArray<T>::module(TArray< complex<T> > const & a)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::module(TArray< complex<T> >& a) Not Allocated Array a !");
+  TArray<T> ra;
+  ra.ReSize(a);
+
+  const complex<T> * pe;
+  T * po;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    po = ra.Data();
+    for(k=0; k<maxx; k+=step )  
+      po[k] = sqrt((double)(pe[k].real()*pe[k].real()+pe[k].imag()*pe[k].imag()));
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      po = ra.DataBlock().Begin()+ra.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  
+	po[k] = sqrt((double)(pe[k].real()*pe[k].real()+pe[k].imag()*pe[k].imag()));
+    }
+  }
+  return(ra);
+}
+
+
+//! Returns the phase of the complex input array.
+/*!
+  \param a : input complex array
+  \return Return a new array filled with the phase of the input complex array elements 
+*/
+
+template <class T>
+TArray<T> ComplexMathArray<T>::phase(TArray< complex<T> > const & a)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("ComplexMathArray< complex<T> >::phase(TArray< complex<T> >& a) Not Allocated Array a !");
+  TArray<T> ra;
+  ra.ReSize(a);
+
+  const complex<T> * pe;
+  T * po;
+  sa_size_t j,k;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    sa_size_t step = a.AvgStep(); 
+    sa_size_t maxx = a.Size()*step;
+    pe = a.Data();
+    po = ra.Data();
+    for(k=0; k<maxx; k+=step )  
+      po[k] = atan2((double)pe[k].imag(), (double)pe[k].real());
+  }
+  else {    // Non regular data spacing ...
+    int_4 ka = a.MaxSizeKA();
+    sa_size_t step = a.Step(ka);
+    sa_size_t gpas = a.Size(ka)*step;
+    sa_size_t naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
+      po = ra.DataBlock().Begin()+ra.Offset(ka,j);
+      for(k=0; k<gpas; k+=step)  
+	po[k] = atan2((double)pe[k].imag(), (double)pe[k].real());
+    }
+  }
+  return(ra);
+}
+
+
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
 #pragma define_template MathArray<r_4>
 #pragma define_template MathArray<r_8>
+#pragma define_template ComplexMathArray<r_4>
+#pragma define_template ComplexMathArray<r_8>
 #endif
 
@@ -161,3 +487,5 @@
 template class MathArray<r_4>;
 template class MathArray<r_8>;
+template class ComplexMathArray<r_4>;
+template class ComplexMathArray<r_8>;
 #endif
Index: trunk/SophyaLib/TArray/matharr.h
===================================================================
--- trunk/SophyaLib/TArray/matharr.h	(revision 1493)
+++ trunk/SophyaLib/TArray/matharr.h	(revision 1517)
@@ -110,4 +110,69 @@
   { MathArray<T> ma;  double mean, sig;  return( ma.MeanSigma(a, mean, sig) ); }
 
+
+//! Class for simple mathematical operation on complex arrays 
+template <class T>
+class ComplexMathArray {
+public:
+// Applying a function 
+  // Replaces the input array content with the result f(x)
+  virtual TArray< complex<T> >& ApplyFunctionInPlace(TArray< complex<T> >& a, 
+						     Arr_ComplexDoubleFunctionOfX f);
+  // Creates a new array and fills it with f(x)
+  virtual TArray< complex<T> > ApplyFunction(TArray< complex<T> > const & a, 
+					     Arr_ComplexDoubleFunctionOfX f);
+
+  // Creates a new array and fills it with f(x)
+  virtual TArray< complex<T> > FillFrom(TArray<T> const & p_real,
+					TArray<T> const & p_imag);
+
+  // Returns real-imaginary part of the complex array
+  virtual TArray<T> real(TArray< complex<T> >  const & a);
+  virtual TArray<T> imag(TArray< complex<T> >  const & a);
+  
+  // Returns module and phase of the complex input array 
+  virtual TArray<T> module2(TArray< complex<T> >  const & a);
+  virtual TArray<T> module(TArray< complex<T> >  const & a);
+  virtual TArray<T> phase(TArray< complex<T> >  const & a);  
+};
+
+/*! \ingroup TArray \fn real(const TArray< complex<T> >&)
+  \brief Return the \b real part of the input complex array \b a */
+template <class T>
+inline TArray<T> real(const TArray< complex<T> >& a)
+  { ComplexMathArray<T> cma;  return( cma.real(a) ); }
+
+/*! \ingroup TArray \fn imag(const TArray< complex<T> >&)
+  \brief Return the \b imaginary part of the input complex array \b a */
+template <class T>
+inline TArray<T> imag(const TArray< complex<T> >& a)
+  { ComplexMathArray<T> cma;  return( cma.imag(a) ); }
+
+/*! \ingroup TArray \fn module2(const TArray< complex<T> >&)
+  \brief Return the \b module squared of the input complex array \b a */
+template <class T>
+inline TArray<T> module2(const TArray< complex<T> >& a)
+  { ComplexMathArray<T> cma;  return( cma.module2(a) ); }
+
+/*! \ingroup TArray \fn module(const TArray< complex<T> >&)
+  \brief Return the \b module of the input complex array \b a */
+template <class T>
+inline TArray<T> module(const TArray< complex<T> >& a)
+  { ComplexMathArray<T> cma;  return( cma.module(a) ); }
+
+/*! \ingroup TArray \fn phase(const TArray< complex<T> >&)
+  \brief Return the \b phase of the input complex array \b a */
+template <class T>
+inline TArray<T> phase(const TArray< complex<T> >& a)
+  { ComplexMathArray<T> cma;  return( cma.phase(a) ); }
+
+/*! \ingroup TArray \fn ComplexArray(const TArray<T> &, const TArray<T> &)
+  \brief Return a complex array, with real and imaginary parts filled from the arguments  */
+template <class T>
+inline TArray< complex<T> > ComplexArray(const TArray<T> & p_real,
+					 const TArray<T> & p_imag)
+  { ComplexMathArray<T> cma;  return( cma.FillFrom(p_real, p_imag) ); }
+
+
 } // Fin du namespace
 
Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 1493)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 1517)
@@ -131,4 +131,5 @@
   : BaseArray() , mNDBlock(a.mNDBlock, share) 
 {
+  if (a.NbDimensions() == 0)  return;  
   string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
   if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
@@ -141,4 +142,5 @@
   : BaseArray() , mNDBlock() 
 {
+  if (a.NbDimensions() == 0)  return;  
   string exmsg = "TArray<T>::TArray(const BaseArray&)";
   if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
@@ -256,5 +258,5 @@
  */
 template <class T>
-void TArray<T>::ReSize(const TArray<T>& a)
+void TArray<T>::ReSize(const BaseArray& a)
 {
   if (arrtype_ != 0) {
@@ -1031,4 +1033,24 @@
 }
 
+//! Fill the array, decoding the ASCII input stream
+/*!
+  \param is : input stream (ASCII)
+ */
+template <class T>
+void TArray<T>::ReadASCII(istream& is)
+{
+  cerr << " TArray<T>::ReadASCII() - Pas encore implemente - Reza 12/6/2001 " << endl;
+}
+
+//! Writes the array content to the output stream, (in ASCII) 
+/*!
+  \param os : output stream (ASCII)
+ */
+template <class T>
+void TArray<T>::WriteASCII(ostream& os) const
+{
+  Print(os, Size(), true);
+}
+
 
 
Index: trunk/SophyaLib/TArray/tarray.h
===================================================================
--- trunk/SophyaLib/TArray/tarray.h	(revision 1493)
+++ trunk/SophyaLib/TArray/tarray.h	(revision 1517)
@@ -58,10 +58,10 @@
 
   void ReSize(int_4 ndim, sa_size_t * siz, sa_size_t step=1);
-  void ReSize(const TArray<T>& a);
+  void ReSize(const BaseArray& a);
   //! a synonym (alias) for method ReSize(int_4, ...)
   inline void SetSize(int_4 ndim, sa_size_t * siz, sa_size_t step=1)
 		{ ReSize(ndim, siz, step); }
-  //! a synonym (alias) for method ReSize(const TArray<T>&)
-  inline void SetSize(const TArray<T>& a)
+  //! a synonym (alias) for method ReSize(const BaseArray&)
+  inline void SetSize(const BaseArray& a)
 		{ ReSize(a); }
   void Realloc(int_4 ndim, sa_size_t * siz, sa_size_t step=1, bool force=false);
@@ -180,4 +180,8 @@
   virtual void   Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
 
+// Lecture,Ecriture sur fichier ASCII
+  virtual void   ReadASCII(istream& is);
+  virtual void   WriteASCII(ostream& os) const;
+
 //  Pour la gestion de persistance
   friend class  FIO_TArray<T>;
@@ -195,4 +199,11 @@
 inline ostream& operator << (ostream& os, const TArray<T>& a)
                             { a.Print(os);    return(os);    }
+
+//   Surcharge d'operateur >> 
+//! Decodes the ASCII input stream \b is , filling TArray \b a elements
+template <class T>
+inline istream& operator >> (istream& is, TArray<T>& a)
+                            { a.ReadASCII(os);    return(is);    }
+
 
 ////////////////////////////////////////////////////////////////
Index: trunk/SophyaLib/TArray/utilarr.h
===================================================================
--- trunk/SophyaLib/TArray/utilarr.h	(revision 1493)
+++ trunk/SophyaLib/TArray/utilarr.h	(revision 1517)
@@ -18,12 +18,23 @@
 /*! \ingroup TArray
   \typedef Arr_DoubleFunctionOfX
-  \brief define a function of float which returns a double
+  \brief define a function of double which returns a double
 */
 typedef double (* Arr_DoubleFunctionOfX) (double x);
 /*! \ingroup TArray
   \typedef Arr_FloatFunctionOfX
-  \brief define a function of float which returns a double
+  \brief define a function of float which returns a float
 */
 typedef float  (* Arr_FloatFunctionOfX)  (float x);
+
+/*! \ingroup TArray
+  \typedef Arr_ComplexDoubleFunctionOfX
+  \brief define a function of a complex<double> which returns a complex<double>
+*/
+typedef complex<double> (* Arr_ComplexDoubleFunctionOfX) (complex<double> x);
+/*! \ingroup TArray
+  \typedef Arr_ComplexFloatFunctionOfX
+  \brief define a function of complex<float> which returns a complex<float>
+*/
+typedef float  (* Arr_ComplexFloatFunctionOfX)  (complex<float> x);
 
 //////////////////////////////////////////////////////////
