Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 785)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 787)
@@ -10,26 +10,4 @@
 
 
-// Variables statiques globales
-template <class T>
-char * TArray<T>::ck_op_msg_[6] = {"???", "Size(int )", "IsPacked(int )", 
-				 "Stride(int )", "ElemCheckBound()", "operator()" };
-template <class T>
-uint_4 TArray<T>::max_nprt_ = 50;
-
-// Methodes statiques globales
-template <class T>
-void TArray<T>::SetMaxPrint(uint_4 nprt)
-{
-  max_nprt_ = nprt;
-}
-
-
-template <class T>
-uint_8 TArray<T>::ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset)
-{
-  uint_8 rs = step;
-  for(int k=0; k<ndim; k++) rs *= siz[k];
-  return(rs+offset);
-}
 
 // -------------------------------------------------------
@@ -40,21 +18,12 @@
 template <class T>
 TArray<T>::TArray()
-  : mNDBlock() , mInfo(NULL)
+  : mNDBlock() , BaseArray()
 // Default constructor
 {
-  ndim_ = 0;
-  for(int k=0; k<TARRAY_MAXNDIMS; k++) step_[k] = size_[k] = 0;
-  totsize_ = 0;
-  minstep_ = 0;
-  moystep_ = 0;
-  offset_ = 0;
-  // Default for matrices : Lines = Y , Columns = X
-  macoli_ = 0;
-  marowi_ = 1;
 }
 
 template <class T>
 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
-  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , mInfo(NULL)
+  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , BaseArray()
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
@@ -64,7 +33,7 @@
 template <class T>
 TArray<T>::TArray(uint_4 nx, uint_4 ny=1, uint_4 nz=1)
-  : mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)) , mInfo(NULL)
-{
-  uint_4 size[TARRAY_MAXNDIMS];
+  : mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)) , BaseArray()
+{
+  uint_4 size[BASEARRAY_MAXNDIMS];
   size[0] = nx;   size[1] = ny;   size[2] = nz;
   int ndim = 1;
@@ -78,5 +47,5 @@
 template <class T>
 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
-  : mNDBlock(db, share) , mInfo(NULL)
+  : mNDBlock(db, share) , BaseArray()
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
@@ -87,5 +56,5 @@
 template <class T>
 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
-  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) , mInfo(NULL)
+  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) , BaseArray()
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
@@ -95,5 +64,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a)
-  : mNDBlock(a.mNDBlock) , mInfo(NULL)
+  : mNDBlock(a.mNDBlock) , BaseArray()
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&)";
@@ -104,5 +73,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a, bool share)
-  : mNDBlock(a.mNDBlock, share) , mInfo(NULL)
+  : mNDBlock(a.mNDBlock, share) , BaseArray()
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
@@ -154,80 +123,60 @@
 }
 
-template <class T>
-bool TArray<T>::CompareSizes(const TArray<T>& a)
-{
-  if (ndim_ != a.ndim_)  return(false);
-  for(int k=0; k<ndim_; k++) 
-    if (size_[k] != a.size_[k])  return(false);
-  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) )  return(false);
-  return(true);
-}
-
-template <class T>
-TArray<T>& TArray<T>::CompactDim()
-{
-  if (ndim_ < 2)  return(*this);
+
+template <class T>
+TArray<T>& TArray<T>::CompactAllDimensions()
+{
+  CompactAllDim();
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::CompactTrailingDimensions()
+{
+  CompactTrailingDim();
+  return(*this);
+}
+
+template <class T>
+double TArray<T>::ValueAtPosition(uint_8 ip) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)  throw( ParmError("TArray<T>::ValueAtPosition(uint_8 ip) Out-of-bound Error") );
+#endif
+  return( (double)(*(mNDBlock.Begin()+Offset(ip))) );
+}
+
+// For complex values, we return the module of the complex number
+double TArray< complex<r_4> >::ValueAtPosition(uint_8 ip) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)  throw( ParmError("TArray<T>::ValueAtPosition(uint_8 ip) Out-of-bound Error") );
+#endif
+  complex<r_4> c = *(mNDBlock.Begin()+Offset(ip));
+  double cr = (double)(c.real());
+  double ci = (double)(c.imag());
+  return( sqrt(cr*cr+ci*ci) );
+}
+
+double TArray< complex<r_8> >::ValueAtPosition(uint_8 ip) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)  throw( ParmError("TArray<T>::ValueAtPosition(uint_8 ip) Out-of-bound Error") );
+#endif
+  complex<r_8> c = *(mNDBlock.Begin()+Offset(ip));
+  double cr = (double)(c.real());
+  double ci = (double)(c.imag());
+  return( sqrt(cr*cr+ci*ci) );
+}
+
+
+// SubArrays 
+template <class T>
+TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
+{
   uint_4 ndim = 0;
-  uint_4 size[TARRAY_MAXNDIMS];
-  uint_4 step[TARRAY_MAXNDIMS];
-  for(int k=0; k<ndim_; k++) {
-    if (size_[k] < 2)  continue;
-    size[ndim] = size_[k];
-    step[ndim] = step_[k];
-    ndim++;
-  }
-  if (ndim == 0)  {
-    size[0] = size_[0];
-    step[0] = step_[0];
-    ndim = 1;
-  }
-  string exmsg = "TArray<T>::CompactDim() ";
-  if (!UpdateSizes(ndim, size, step, offset_, exmsg))  throw( ParmError(exmsg) );
-  return(*this);
-}
-
-
-template <class T>
-uint_4 TArray<T>::MinStepKA() const
-{
-  for(int ka=0; ka<ndim_; ka++)
-    if (step_[ka] == minstep_) return(ka);
-  return(0);
-}
-
-template <class T>
-uint_4 TArray<T>::MaxSizeKA() const
-{
-  uint_4 ka = 0;
-  uint_4 mx = size_[0];
-  for(int k=0; k<ndim_; k++)  
-    if (size_[k] > mx) {  ka = k;  mx = size_[k];  }
-  return(ka);
-}
-
-
-//  Acces lineaire aux elements ....  Calcul d'offset
-
-template <class T>
-uint_8 TArray<T>::Offset(uint_8 ip) const
-{
-if (ip == 0)  return(offset_);
-uint_4 idx[TARRAY_MAXNDIMS];
-int k;
-uint_8 rest = ip;
-for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
-uint_8 off = offset_;
-for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
-return (off);
-}
-
-// SubArrays 
-template <class T>
-TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
-{
-  uint_4 ndim = 0;
-  uint_4 size[TARRAY_MAXNDIMS];
-  uint_4 step[TARRAY_MAXNDIMS];
-  uint_4 pos[TARRAY_MAXNDIMS];
+  uint_4 size[BASEARRAY_MAXNDIMS];
+  uint_4 step[BASEARRAY_MAXNDIMS];
+  uint_4 pos[BASEARRAY_MAXNDIMS];
   size[0] = rx.Size();
   size[1] = ry.Size();
@@ -251,4 +200,5 @@
   TArray<T> ra;
   SubArray(ra, ndim, size, pos, step); 
+  ra.DataBlock().Share(this->DataBlock());
   ra.SetTemp(true);
   return(ra);
@@ -523,92 +473,12 @@
 
 // ----------------------------------------------------
-//          Application d'une fonction
-// ----------------------------------------------------
-template <class T>
-TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
-{
-  T * pe;
-  uint_8 j,k;
-  if (AvgStep() > 0)   {  // regularly spaced elements
-    uint_8 step = AvgStep();
-    uint_8 maxx = totsize_*step;
-    pe = Data();
-    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((double)pe[k]));
-  }
-  else {    // Non regular data spacing ...
-    uint_4 ka = MaxSizeKA();
-    uint_8 step = Step(ka);
-    uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
-      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
-    }
-  }
-  return(*this);
-}
-
-template <class T>
-TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
-{
-  T * pe;
-  uint_8 j,k;
-  if (AvgStep() > 0)   {  // regularly spaced elements
-    uint_8 step = AvgStep();
-    uint_8 maxx = totsize_*step;
-    pe = Data();
-    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((float)pe[k]));
-  }
-  else {    // Non regular data spacing ...
-    uint_4 ka = MaxSizeKA();
-    uint_8 step = Step(ka);
-    uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
-      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
-    }
-  }
-  return(*this);
-}
-
-TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_DoubleFunctionOfX f)
-{
-  throw( NotAvailableOperation("TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
-}
-
-TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_FloatFunctionOfX f)
-{
-  throw(NotAvailableOperation( "TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
-}
-
-TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_DoubleFunctionOfX f)
-{
-  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
-}
-
-TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_FloatFunctionOfX f)
-{
-  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
-}
-
-// ----------------------------------------------------
 //       Impression, etc ...
 // ----------------------------------------------------
 
 template <class T>
-void TArray<T>::Show(ostream& os, bool si) const
-{
-  os << " TArray< " << typeid(T).name() << " NDim= " << ndim_
-     << "(" << typeid(*this).name() << " )" << endl;
-  os << "TotSize=" << totsize_ << " Size(X*Y*...)= " ;
-  for(int k=0; k<ndim_; k++) { 
-    os << size_[k];
-    if (k<ndim_-1)  os << " x ";
-  }
-  os << endl;
-  os <<  " Step(X Y ...)= " ;
-  for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
-  os << " Offset= " << offset_  << "\n " << endl;
-  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
- 
+string TArray<T>::DataType() const
+{
+  string rs = typeid(T).name();
+  return(rs);
 }
 
@@ -643,132 +513,4 @@
 }
 
-template <class T>
-DVList& TArray<T>::Info()
-{
-if (mInfo == NULL)  mInfo = new DVList;
-return(*mInfo);
-}
-
-
-template <class T>
-bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
-{
-  if (ndim >= TARRAY_MAXNDIMS) {
-    exmsg += " NDim Error";  return false;
-  }
-  if (step < 1) {
-    exmsg += " Step(=0) Error";  return false;
-  }
-
-  minstep_ = moystep_ = step;
-
-  // Flagging bad updates ...
-  ndim_ = 0;
-
-  totsize_ = 1;
-  int k;
-  for(k=0; k<TARRAY_MAXNDIMS; k++) {
-    size_[k] = 1;
-    step_[k] = 0;
-  }
-  for(k=0; k<ndim; k++) {
-    size_[k] = siz[k] ;
-    step_[k] = totsize_*step;
-    totsize_ *= size_[k];
-  }
-  if (totsize_ < 1) {
-    exmsg += " Size Error";  return false;
-  }
-  offset_ = offset;
-  // Default for matrices : Lines = Y , Columns = X
-  macoli_ = 0;
-  marowi_ = 1;
-  // Update OK 
-  ndim_ = ndim;
-  return true;
-}
-
-template <class T>
-bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
-{
-  if (ndim >= TARRAY_MAXNDIMS) {
-    exmsg += " NDim Error";  return false;
-  }
-
-  // Flagging bad updates ...
-  ndim_ = 0;
-
-  totsize_ = 1;
-  int k;
-  for(k=0; k<TARRAY_MAXNDIMS; k++) {
-    size_[k] = 1;
-    step_[k] = 0;
-  }
-  uint_4 minstep = step[0];
-  for(k=0; k<ndim; k++) {
-    size_[k] = siz[k] ;
-    step_[k] = step[k];
-    totsize_ *= size_[k];
-    if (step_[k] < minstep)  minstep = step_[k];
-  }
-  if (minstep < 1) {
-    exmsg += " Step(=0) Error";  return false;
-  }
-  if (totsize_ < 1) {
-    exmsg += " Size Error";  return false;
-  }
-  uint_8 plast = 0;
-  for(k=0; k<ndim; k++)   plast += (siz[k]-1)*step[k];
-  if (plast == minstep*totsize_ )  moystep_ = minstep;
-  else moystep_ = 0;
-  minstep_ = minstep;
-  offset_ = offset;
-  // Default for matrices : Lines = Y , Columns = X
-  macoli_ = 0;
-  marowi_ = 1;
-  // Update OK 
-  ndim_ = ndim;
-  return true;
-}
-
-template <class T>
-bool TArray<T>::UpdateSizes(const TArray<T>& a, string & exmsg)
-{
-  if (a.ndim_ >= TARRAY_MAXNDIMS) {
-    exmsg += " NDim Error";  return false;
-  }
-
-  // Flagging bad updates ...
-  ndim_ = 0;
-
-  totsize_ = 1;
-  int k;
-  for(k=0; k<TARRAY_MAXNDIMS; k++) {
-    size_[k] = 1;
-    step_[k] = 0;
-  }
-  uint_4 minstep = a.step_[0];
-  for(k=0; k<a.ndim_; k++) {
-    size_[k] = a.size_[k] ;
-    step_[k] = a.step_[k];
-    totsize_ *= size_[k];
-    if (step_[k] < minstep)  minstep = step_[k];
-  }
-  if (minstep < 1) {
-    exmsg += " Step(=0) Error";  return false;
-  }
-  if (totsize_ < 1) {
-    exmsg += " Size Error";  return false;
-  }
-
-  minstep_ = a.minstep_;
-  moystep_ = a.moystep_;
-  offset_ = a.offset_;
-  macoli_ = a.macoli_;
-  marowi_ = a.marowi_;
-  // Update OK 
-  ndim_ = a.ndim_;
-  return true;
-}
 
 template <class T>
@@ -789,24 +531,4 @@
 
 
-template <class T>
-void TArray<T>::SubArray(TArray<T> & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step)
-{
-  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
-  int k;
-  for(k=0; k<ndim; k++) 
-    if ( (siz[k]*step[k]+pos[k]) > size_[k] )  
-      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error") );
-  (ra.mNDBlock).Share(mNDBlock);
-  uint_8 offset = offset_;
-  for(k=0; k<ndim_; k++) { 
-    offset += pos[k]*step_[k]; 
-    step[k] *= step_[k];
-  }
-  string exm = "TArray<T>::SubArray() ";
-  if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
-     throw( ParmError(exm) );
-  (ra.mNDBlock).Share(mNDBlock);
-  return;
-}
 
 ///////////////////////////////////////////////////////////////
