Index: trunk/SophyaLib/TArray/fioarr.cc
===================================================================
--- trunk/SophyaLib/TArray/fioarr.cc	(revision 784)
+++ trunk/SophyaLib/TArray/fioarr.cc	(revision 785)
@@ -74,4 +74,5 @@
   is.Get(dobj->step_, TARRAY_MAXNDIMS);
   is.Get(dobj->minstep_);
+  is.Get(dobj->moystep_);
   is.Get(dobj->offset_);
   is.Get(dobj->marowi_);
@@ -101,4 +102,5 @@
   os.Put(dobj->step_, TARRAY_MAXNDIMS);
   os.Put(dobj->minstep_);
+  os.Put(dobj->moystep_);
   os.Put(dobj->offset_);
   os.Put(dobj->marowi_);
Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 784)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 785)
@@ -9,10 +9,4 @@
 
 
-// Classe utilitaires 
-Sequence::Sequence(double start, double step)
-{
-  start_ = start;
-  step_ = step;
-}
 
 // Variables statiques globales
@@ -53,4 +47,5 @@
   totsize_ = 0;
   minstep_ = 0;
+  moystep_ = 0;
   offset_ = 0;
   // Default for matrices : Lines = Y , Columns = X
@@ -123,27 +118,4 @@
 
 template <class T>
-TArray<T> TArray<T>::SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos)
-{
-  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
-  int k;
-  for(k=0; k<ndim; k++) 
-    if ( (siz[k]+pos[k]) > size_[k] )  
-      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (1)") );
-  for(k=ndim; k<ndim_; k++)   
-    if ( (pos[k]) >= size_[k] )  
-      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (2)") );
-  TArray<T> ra(*this);
-  ra.ndim_ = ndim;
-  for(k=ndim; k<TARRAY_MAXNDIMS; k++) {
-    ra.size_[k] = 1;
-    ra.step_[k] = 0;
-  }
-  ra.offset_ = offset_;
-  for(k=0; k<ndim_; k++) ra.offset_ += pos[k]*step_[k]; 
-  ra.SetTemp(true);
-  return(ra);   
-}
-
-template <class T>
 TArray<T>& TArray<T>::operator = (const TArray<T>& a)
 {
@@ -192,4 +164,95 @@
 }
 
+template <class T>
+TArray<T>& TArray<T>::CompactDim()
+{
+  if (ndim_ < 2)  return(*this);
+  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];
+  size[0] = rx.Size();
+  size[1] = ry.Size();
+  size[2] = rz.Size();
+  size[3] = rt.Size();
+  size[4] = ru.Size();
+
+  step[0] = rx.Step();
+  step[1] = ry.Step();
+  step[2] = rz.Step();
+  step[3] = rt.Step();
+  step[4] = ru.Step();
+
+  pos[0] = rx.Start();
+  pos[1] = ry.Start();
+  pos[2] = rz.Start();
+  pos[3] = rt.Start();
+  pos[4] = ru.Start();
+
+  ndim = ndim_;
+  TArray<T> ra;
+  SubArray(ra, ndim, size, pos, step); 
+  ra.SetTemp(true);
+  return(ra);
+}
+
 //   ...... Operation de calcul sur les tableaux ......
 //              ------- Attention --------
@@ -197,9 +260,22 @@
 //  Possibilite de // , vectorisation 
 template <class T>
-TArray<T>& TArray<T>::operator = (Sequence & seq)
-{
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
+TArray<T>& TArray<T>::operator = (Sequence seq)
+{
+  T * pe;
+  uint_8 j,k;
+  if (AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    pe = Data();
+    for(k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = MaxSizeKA();
+    uint_8 step = Step(ka);
+    uint_8 gpas = Size(ka);
+    for(j=0; j<totsize_; j += Size(ka))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      for(k=0; k<gpas; k++)  pe[k*step] = seq(j+k);
+    }
+  }
   return(*this);
 }
@@ -210,7 +286,21 @@
 TArray<T>& TArray<T>::operator = (T x)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = x;
+  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] = x;
+  }
+  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] = x;
+    }
+  }
   return(*this);
 }
@@ -219,7 +309,21 @@
 TArray<T>& TArray<T>::operator += (T x)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] += x;
+  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] += x;
+  }
+  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] += x;
+    }
+  }
   return(*this);
 }
@@ -228,7 +332,21 @@
 TArray<T>& TArray<T>::operator -= (T x)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] -= x;
+  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] -= x;
+  }
+  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] -= x;
+    }
+  }
   return(*this);
 }
@@ -237,7 +355,21 @@
 TArray<T>& TArray<T>::operator *= (T x)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_; k += step)  pe[k] *= x;
+  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] *= x;
+  }
+  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] *= x;
+    }
+  }
   return(*this);
 }
@@ -246,7 +378,21 @@
 TArray<T>& TArray<T>::operator /= (T x)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] /= x;
+  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] /= x;
+  }
+  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] /= x;
+    }
+  }
   return(*this);
 }
@@ -258,13 +404,26 @@
 {
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator += (const TArray<T>&)")) ;
-  T * pe = Data();
-  const T * pea = a.Data();
-  uint_8 step = Step();
-  uint_8 stepa = a.Step();
-  uint_8 k, ka;
-  k = ka = 0;
-  for(k=0; k<totsize_; k += step)  {
-    pe[k] += pea[ka];
-    ka += stepa;
+
+  T * pe;
+  const T * pea;
+  uint_8 j,k,ka;
+  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 stepa = a.AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    pea = a.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] += pea[ka] ;
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ax = MaxSizeKA();
+    uint_8 step = Step(ax);
+    uint_8 stepa = a.Step(ax);
+    uint_8 gpas = Size(ax)*step;
+    for(j=0; j<totsize_; j += Size(ax))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      pea = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka]; 
+    }
   }
   return(*this);
@@ -275,13 +434,26 @@
 {
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
-  T * pe = Data();
-  const T * pea = a.Data();
-  uint_8 step = Step();
-  uint_8 stepa = a.Step();
-  uint_8 k, ka;
-  k = ka = 0;
-  for(k=0; k<totsize_; k += step)  {
-    pe[k] -= pea[ka];
-    ka += stepa;
+
+  T * pe;
+  const T * pea;
+  uint_8 j,k,ka;
+  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 stepa = a.AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    pea = a.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] -= pea[ka] ;
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ax = MaxSizeKA();
+    uint_8 step = Step(ax);
+    uint_8 stepa = a.Step(ax);
+    uint_8 gpas = Size(ax)*step;
+    for(j=0; j<totsize_; j += Size(ax))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      pea = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka]; 
+    }
   }
   return(*this);
@@ -293,13 +465,26 @@
 {
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
-  T * pe = Data();
-  const T * pea = a.Data();
-  uint_8 step = Step();
-  uint_8 stepa = a.Step();
-  uint_8 k, ka;
-  k = ka = 0;
-  for(k=0; k<totsize_; k += step)  {
-    pe[k] *= pea[ka];
-    ka += stepa;
+
+  T * pe;
+  const T * pea;
+  uint_8 j,k,ka;
+  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 stepa = a.AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    pea = a.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] *= pea[ka] ;
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ax = MaxSizeKA();
+    uint_8 step = Step(ax);
+    uint_8 stepa = a.Step(ax);
+    uint_8 gpas = Size(ax)*step;
+    for(j=0; j<totsize_; j += Size(ax))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      pea = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka]; 
+    }
   }
   return(*this);
@@ -310,13 +495,26 @@
 {
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
-  T * pe = Data();
-  const T * pea = a.Data();
-  uint_8 step = Step();
-  uint_8 stepa = a.Step();
-  uint_8 k, ka;
-  k = ka = 0;
-  for(k=0; k<totsize_; k += step)  {
-    pe[k] /= pea[ka];
-    ka += stepa;
+
+  T * pe;
+  const T * pea;
+  uint_8 j,k,ka;
+  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 stepa = a.AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    pea = a.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] /= pea[ka] ;
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ax = MaxSizeKA();
+    uint_8 step = Step(ax);
+    uint_8 stepa = a.Step(ax);
+    uint_8 gpas = Size(ax)*step;
+    for(j=0; j<totsize_; j += Size(ax))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      pea = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka]; 
+    }
   }
   return(*this);
@@ -330,7 +528,21 @@
 TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((double)pe[k]));
+  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);
 }
@@ -339,7 +551,21 @@
 TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
 {
-  T * pe = Data();
-  uint_8 step = Step();
-  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((float)pe[k]));
+  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);
 }
@@ -382,5 +608,5 @@
   os <<  " Step(X Y ...)= " ;
   for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
-  os << "\n " << endl;
+  os << " Offset= " << offset_  << "\n " << endl;
   if (si && (mInfo != NULL)) os << (*mInfo) << endl;
  
@@ -395,13 +621,14 @@
   int k0,k1,k2,k3,k4;
   for(k4=0; k4<size_[4]; k4++) {
-    if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4;
+    if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
     for(k3=0; k3<size_[3]; k3++) {
-      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3;
+      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3 << endl;
       for(k2=0; k2<size_[2]; k2++) {
-	if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2;
+	if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2 << endl;
 	for(k1=0; k1<size_[1]; k1++) {
-	  if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1;
+	  if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1 << endl;
 	  for(k0=0; k0<size_[0]; k0++) {
-            os << Elem(k0, k1, k2, k3, k4) << ", ";    npr++;
+            if(k0 > 0) os << ", ";  
+            os << Elem(k0, k1, k2, k3, k4);     npr++;
 	    if (npr >= maxprt) {
 	      if (npr < totsize_)  os << "\n     .... " << endl; return;
@@ -434,5 +661,5 @@
   }
 
-  minstep_ = 1;
+  minstep_ = moystep_ = step;
 
   // Flagging bad updates ...
@@ -478,12 +705,12 @@
     step_[k] = 0;
   }
-  minstep_ = step[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) {
+    if (step_[k] < minstep)  minstep = step_[k];
+  }
+  if (minstep < 1) {
     exmsg += " Step(=0) Error";  return false;
   }
@@ -491,4 +718,9 @@
     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
@@ -499,4 +731,5 @@
   return true;
 }
+
 template <class T>
 bool TArray<T>::UpdateSizes(const TArray<T>& a, string & exmsg)
@@ -515,12 +748,12 @@
     step_[k] = 0;
   }
-  minstep_ = a.step_[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) {
+    if (step_[k] < minstep)  minstep = step_[k];
+  }
+  if (minstep < 1) {
     exmsg += " Step(=0) Error";  return false;
   }
@@ -529,4 +762,6 @@
   }
 
+  minstep_ = a.minstep_;
+  moystep_ = a.moystep_;
   offset_ = a.offset_;
   macoli_ = a.macoli_;
@@ -553,4 +788,25 @@
 }
 
+
+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;
+}
 
 ///////////////////////////////////////////////////////////////
Index: trunk/SophyaLib/TArray/tarray.h
===================================================================
--- trunk/SophyaLib/TArray/tarray.h	(revision 784)
+++ trunk/SophyaLib/TArray/tarray.h	(revision 785)
@@ -14,4 +14,5 @@
 #include "ndatablock.h"
 #include "dvlist.h"
+#include "utilarr.h"
 
 
@@ -20,18 +21,4 @@
 
 namespace SOPHYA {
-
-/* Quelques utilitaires */
-typedef double (* Arr_DoubleFunctionOfX) (double x);
-typedef float  (* Arr_FloatFunctionOfX)  (float x);
-
-class Sequence { 
-public:
-  Sequence (double start=0., double step=1.);
-  inline double & Start() { return start_; }
-  inline double & Step() { return step_; }
-  inline double operator () (uint_4 k) { return(start_+(double)k*step_); }
-protected:
-  double start_, step_;
-};
 
 // Forward declaration
@@ -56,6 +43,4 @@
   virtual TArray<T>& operator = (const TArray<T>& a);
 
-  // Extraction de sous-tableau
-  virtual TArray<T> SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos);
 
   // Gestion taille/Remplissage
@@ -65,5 +50,7 @@
   // Returns true if ndim and sizes are equal
   virtual bool CompareSizes(const TArray<T>& a);
- 
+  // Compacts size=1 array dimensions
+  virtual TArray<T>& CompactDim();
+
   // Array dimensions
   inline uint_4 NbDimensions() const { return( ndim_ ); }
@@ -75,6 +62,8 @@
   inline uint_4 Size(int ka) const { return(size_[CheckDI(ka,1)]); }
 
+  uint_4 MaxSizeKA() const ;
+
   // memory organisation - packing information
-  inline bool   IsPacked() const { return(minstep_ == 1); }
+  inline bool   IsPacked() const { return(moystep_ == 1); }
   inline bool   IsPackedX() const { return(step_[0] == 1); }
   inline bool   IsPackedY() const { return(step_[1] == 1); }
@@ -82,5 +71,6 @@
   inline bool   IsPacked(int ka) const { return(step_[CheckDI(ka,2)] == 1); }
 
-  inline uint_4 Step() const  { return(minstep_); }
+  inline uint_4 MinStep() const  { return(minstep_); }
+  inline uint_4 AvgStep() const  { return(moystep_); }
   inline uint_4 StepX() const { return(step_[0]); }
   inline uint_4 StepZ() const { return(step_[1]); }
@@ -88,5 +78,7 @@
   inline uint_4 Step(int ka) const { return(step_[CheckDI(ka,3)]); }
 
-  inline uint_8 Offset() const { return(offset_); }
+  uint_4 MinStepKA() const ;
+
+  uint_8 Offset(uint_8 ip=0) const ;
 
 
@@ -94,4 +86,7 @@
   inline bool   IsTemp(void) const {return mNDBlock.IsTemp();}
   inline void   SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
+
+  // SubArrays
+  virtual TArray<T> operator () (Range rx, Range ry, Range rz=0, Range rt=0, Range ru=0); 
 
   // Acces to data
@@ -100,6 +95,6 @@
   inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0) const ;
   inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0);
-  inline T const& operator[](uint_4 ip) const ;
-  inline T&       operator[](uint_4 ip);
+  inline T const& operator[](uint_8 ip) const ;
+  inline T&       operator[](uint_8 ip);
 
   inline T const& Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
@@ -117,5 +112,5 @@
   inline operator T();
 // Met les elements a une suite de valeurs
-  virtual TArray<T>&  operator = (Sequence & seq);
+  virtual TArray<T>&  operator = (Sequence seq);
 // A = x (tous les elements a x)
   virtual TArray<T>&  operator = (T x);
@@ -162,4 +157,6 @@
   // Share: partage les donnees de "a"
   void Share(const TArray<T>& a);
+  // Extraction de sous-tableau
+  virtual void SubArray(TArray<T> & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step);
 
   uint_4 ndim_;                   // nb of dimensions
@@ -168,4 +165,5 @@
   uint_4 step_[TARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension 
   uint_4 minstep_;                   // minimal step (in any axes)
+  uint_4 moystep_;                   // mean step 0 non regular steps
   uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
   NDataBlock<T> mNDBlock;      // Le bloc des donnees
@@ -254,4 +252,5 @@
 }
 
+
 // --------------------------------------------------
 //        inline element acces methods
@@ -329,21 +328,23 @@
 }
 
-template <class T>
-inline T const& TArray<T>::operator[](uint_4 ip) const
+
+template <class T>
+inline T const& TArray<T>::operator[](uint_8 ip) const
 {
 #ifdef SO_BOUNDCHECKING
   if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
 #endif
-return *(mNDBlock.Begin()+offset_+ip*minstep_);
-}
-
-template <class T>
-inline T & TArray<T>::operator[](uint_4 ip) 
+return *(mNDBlock.Begin()+Offset(ip));
+}
+
+template <class T>
+inline T & TArray<T>::operator[](uint_8 ip) 
 {
 #ifdef SO_BOUNDCHECKING
   if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
 #endif
-return *(mNDBlock.Begin()+offset_+ip*minstep_);
-}
+return *(mNDBlock.Begin()+Offset(ip));
+}
+
 
 template <class T>
Index: trunk/SophyaLib/TArray/utilarr.cc
===================================================================
--- trunk/SophyaLib/TArray/utilarr.cc	(revision 785)
+++ trunk/SophyaLib/TArray/utilarr.cc	(revision 785)
@@ -0,0 +1,35 @@
+//     Utility classes for template numerical arrays
+//                     R. Ansari, C.Magneville   03/2000
+
+#include "machdefs.h"
+#include "utilarr.h"
+
+// Classe utilitaires
+Sequence::Sequence(double start, double step, Arr_DoubleFunctionOfX f)
+{
+  start_ = start;
+  step_ = step;
+  myf_ = f;
+}
+
+double Sequence::operator () (uint_4 k)
+{
+  double x = start_+(double)k*step_;
+  if (myf_)  return(myf_(x));
+  return x;
+}
+
+Range::Range(uint_4 start, uint_4 size, uint_4 step)
+{
+  start_ = start;
+  size_ = (size > 0) ? size : 1;
+  step_ = (step > 0) ? step : 1;
+}
+
+Range & Range::operator = (uint_4 start)
+{
+  start_ = start;
+  size_ = 1;
+  step_ = 1;
+  return (*this);
+}
Index: trunk/SophyaLib/TArray/utilarr.h
===================================================================
--- trunk/SophyaLib/TArray/utilarr.h	(revision 785)
+++ trunk/SophyaLib/TArray/utilarr.h	(revision 785)
@@ -0,0 +1,43 @@
+// This may look like C code, but it is really -*- C++ -*-
+//     Utility classes for template numerical arrays 
+//                     R. Ansari, C.Magneville   03/2000
+
+#ifndef UtilArray_SEEN
+#define UtilArray_SEEN
+
+#include "machdefs.h"
+#include <stdlib.h>
+
+namespace SOPHYA {   
+
+/* Quelques utilitaires pour les tableaux (Array) */             
+
+typedef double (* Arr_DoubleFunctionOfX) (double x);
+typedef float  (* Arr_FloatFunctionOfX)  (float x);
+
+class Sequence {
+public:
+  Sequence (double start=0., double step=1., Arr_DoubleFunctionOfX f=NULL);
+  inline double & Start() { return start_; }
+  inline double & Step() { return step_; }
+  double operator () (uint_4 k);
+protected:
+  double start_, step_;
+  Arr_DoubleFunctionOfX myf_;
+};
+
+class Range {
+public:
+  Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
+  inline uint_4 & Start()  {  return start_; }
+  inline uint_4 & Size()  {  return size_; }
+  inline uint_4 & Step()  {  return step_; }
+  Range & operator = (uint_4 start);
+protected:
+  uint_4 start_, size_, step_;
+};
+
+} // Fin du namespace
+
+#endif
+
