*** ./liboctave/CNDArray.cc.orig	2004-03-26 18:17:12.000000000 +0100
--- ./liboctave/CNDArray.cc	2004-03-26 20:42:05.000000000 +0100
***************
*** 683,688 ****
--- 683,874 ----
    return ::cat_ra(*this, ra_arg, dim, iidx, move);
  }
  
+ static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
+ 
+ ComplexNDArray
+ ComplexNDArray::max (int dim) const
+ {
+   ArrayN<int> dummy_idx;
+   return max (dummy_idx, dim);
+ }
+ 
+ ComplexNDArray
+ ComplexNDArray::max (ArrayN<int>& idx_arg, int dim) const
+ {
+   dim_vector dv = dims ();
+   dim_vector dr = dims ();
+ 
+   if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+     return ComplexNDArray ();
+   
+   dr(dim) = 1;
+ 
+   ComplexNDArray result (dr);
+   idx_arg.resize (dr);
+ 
+   int x_stride = 1;
+   int x_len = dv(dim);
+   for (int i = 0; i < dim; i++)
+     x_stride *= dv(i);
+ 
+   for (int i = 0; i < dr.numel (); i++)
+     {
+       int x_offset;
+       if (x_stride == 1)
+ 	x_offset = i * x_len;
+       else
+ 	{
+ 	  int x_offset2 = 0;
+ 	  x_offset = i;
+ 	  while (x_offset >= x_stride)
+ 	    {
+ 	      x_offset -= x_stride;
+ 	      x_offset2++;
+ 	    }
+ 	  x_offset += x_offset2 * x_stride * x_len;
+ 	}
+ 
+       int idx_j;
+ 
+       Complex tmp_max;
+ 
+       double abs_max = octave_NaN;
+ 
+       for (idx_j = 0; idx_j < x_len; idx_j++)
+ 	{
+ 	  tmp_max = elem (idx_j * x_stride + x_offset);
+ 	  
+ 	  if (! octave_is_NaN_or_NA (tmp_max))
+ 	    {
+ 	      abs_max = ::abs(tmp_max);
+ 	      break;
+ 	    }
+ 	}
+ 
+       for (int j = idx_j+1; j < x_len; j++)
+ 	{
+ 	  Complex tmp = elem (j * x_stride + x_offset);
+ 
+ 	  if (octave_is_NaN_or_NA (tmp))
+ 	    continue;
+ 
+ 	  double abs_tmp = ::abs (tmp);
+ 
+ 	  if (abs_tmp > abs_max)
+ 	    {
+ 	      idx_j = j;
+ 	      tmp_max = tmp;
+ 	      abs_max = abs_tmp;
+ 	    }
+ 	}
+ 
+       if (octave_is_NaN_or_NA (tmp_max))
+ 	{
+ 	  result.elem (i) = Complex_NaN_result;
+ 	  idx_arg.elem (i) = 0;
+ 	}
+       else
+ 	{
+ 	  result.elem (i) = tmp_max;
+ 	  idx_arg.elem (i) = idx_j;
+ 	}
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ ComplexNDArray::min (int dim) const
+ {
+   ArrayN<int> dummy_idx;
+   return min (dummy_idx, dim);
+ }
+ 
+ ComplexNDArray
+ ComplexNDArray::min (ArrayN<int>& idx_arg, int dim) const
+ {
+   dim_vector dv = dims ();
+   dim_vector dr = dims ();
+ 
+   if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+     return ComplexNDArray ();
+   
+   dr(dim) = 1;
+ 
+   ComplexNDArray result (dr);
+   idx_arg.resize (dr);
+ 
+   int x_stride = 1;
+   int x_len = dv(dim);
+   for (int i = 0; i < dim; i++)
+     x_stride *= dv(i);
+ 
+   for (int i = 0; i < dr.numel (); i++)
+     {
+       int x_offset;
+       if (x_stride == 1)
+ 	x_offset = i * x_len;
+       else
+ 	{
+ 	  int x_offset2 = 0;
+ 	  x_offset = i;
+ 	  while (x_offset >= x_stride)
+ 	    {
+ 	      x_offset -= x_stride;
+ 	      x_offset2++;
+ 	    }
+ 	  x_offset += x_offset2 * x_stride * x_len;
+ 	}
+ 
+       int idx_j;
+ 
+       Complex tmp_min;
+ 
+       double abs_min = octave_NaN;
+ 
+       for (idx_j = 0; idx_j < x_len; idx_j++)
+ 	{
+ 	  tmp_min = elem (idx_j * x_stride + x_offset);
+ 	  
+ 	  if (! octave_is_NaN_or_NA (tmp_min))
+ 	    {
+ 	      abs_min = ::abs(tmp_min);
+ 	      break;
+ 	    }
+ 	}
+ 
+       for (int j = idx_j+1; j < x_len; j++)
+ 	{
+ 	  Complex tmp = elem (j * x_stride + x_offset);
+ 
+ 	  if (octave_is_NaN_or_NA (tmp))
+ 	    continue;
+ 
+ 	  double abs_tmp = ::abs (tmp);
+ 
+ 	  if (abs_tmp < abs_min)
+ 	    {
+ 	      idx_j = j;
+ 	      tmp_min = tmp;
+ 	      abs_min = abs_tmp;
+ 	    }
+ 	}
+ 
+       if (octave_is_NaN_or_NA (tmp_min))
+ 	{
+ 	  result.elem (i) = Complex_NaN_result;
+ 	  idx_arg.elem (i) = 0;
+ 	}
+       else
+ 	{
+ 	  result.elem (i) = tmp_min;
+ 	  idx_arg.elem (i) = idx_j;
+ 	}
+     }
+ 
+   return result;
+ }
+ 
  NDArray
  ComplexNDArray::abs (void) const
  {
***************
*** 836,841 ****
--- 1022,1162 ----
    return is;
  }
  
+ // XXX FIXME XXX -- it would be nice to share code among the min/max
+ // functions below.
+ 
+ #define EMPTY_RETURN_CHECK(T) \
+   if (nel == 0)	\
+     return T (dv);
+ 
+ ComplexNDArray
+ min (Complex c, const ComplexNDArray& m)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (c, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ min (const ComplexNDArray& m, Complex c)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (c, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ min (const ComplexNDArray& a, const ComplexNDArray& b)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   if (dv != b.dims ())
+     {
+       (*current_liboctave_error_handler)
+ 	("two-arg min expecting args of same size");
+       return ComplexNDArray ();
+     }
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (a (i), b (i));
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ max (Complex c, const ComplexNDArray& m)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (c, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ max (const ComplexNDArray& m, Complex c)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (c, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ ComplexNDArray
+ max (const ComplexNDArray& a, const ComplexNDArray& b)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   if (dv != b.dims ())
+     {
+       (*current_liboctave_error_handler)
+ 	("two-arg max expecting args of same size");
+       return ComplexNDArray ();
+     }
+ 
+   EMPTY_RETURN_CHECK (ComplexNDArray);
+ 
+   ComplexNDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (a (i), b (i));
+     }
+ 
+   return result;
+ }
+ 
  NDS_CMP_OPS(ComplexNDArray, real, Complex, real)
  NDS_BOOL_OPS(ComplexNDArray, Complex, 0.0)
  
*** ./liboctave/CNDArray.h.orig	2004-03-26 18:17:23.000000000 +0100
--- ./liboctave/CNDArray.h	2004-03-26 18:26:50.000000000 +0100
***************
*** 89,94 ****
--- 89,99 ----
    ComplexNDArray sumsq (int dim = -1) const;
    int cat (const ComplexNDArray& ra_arg, int dim, int iidx, int move);
  
+   ComplexNDArray max (int dim = 0) const;
+   ComplexNDArray max (ArrayN<int>& index, int dim = 0) const;
+   ComplexNDArray min (int dim = 0) const;
+   ComplexNDArray min (ArrayN<int>& index, int dim = 0) const;
+ 
    ComplexNDArray& insert (const NDArray& a, int r, int c);
    ComplexNDArray& insert (const ComplexNDArray& a, int r, int c);
    
***************
*** 130,135 ****
--- 135,148 ----
      : MArrayN<Complex> (d, dv) { }
  };
  
+ extern ComplexNDArray min (const Complex& c, const ComplexNDArray& m);
+ extern ComplexNDArray min (const ComplexNDArray& m, const Complex& c);
+ extern ComplexNDArray min (const ComplexNDArray& a, const ComplexNDArray& b);
+ 
+ extern ComplexNDArray max (const Complex& c, const ComplexNDArray& m);
+ extern ComplexNDArray max (const ComplexNDArray& m, const Complex& c);
+ extern ComplexNDArray max (const ComplexNDArray& a, const ComplexNDArray& b);
+ 
  NDS_CMP_OP_DECLS (ComplexNDArray, Complex)
  NDS_BOOL_OP_DECLS (ComplexNDArray, Complex)
  
*** ./liboctave/dNDArray.cc.orig	2004-03-29 10:27:13.000000000 +0200
--- ./liboctave/dNDArray.cc	2004-03-26 20:39:23.000000000 +0100
***************
*** 656,661 ****
--- 656,811 ----
    MX_ND_REAL_OP_REDUCTION (+= elem (iter_idx), 0);
  }
  
+ NDArray
+ NDArray::max (int dim) const
+ {
+   ArrayN<int> dummy_idx;
+   return max (dummy_idx, dim);
+ }
+ 
+ NDArray
+ NDArray::max (ArrayN<int>& idx_arg, int dim) const
+ {
+   dim_vector dv = dims ();
+   dim_vector dr = dims ();
+ 
+   if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+     return NDArray ();
+   
+   dr(dim) = 1;
+ 
+   NDArray result (dr);
+   idx_arg.resize (dr);
+ 
+   int x_stride = 1;
+   int x_len = dv(dim);
+   for (int i = 0; i < dim; i++)
+     x_stride *= dv(i);
+ 
+   for (int i = 0; i < dr.numel (); i++)
+     {
+       int x_offset;
+       if (x_stride == 1)
+ 	x_offset = i * x_len;
+       else
+ 	{
+ 	  int x_offset2 = 0;
+ 	  x_offset = i;
+ 	  while (x_offset >= x_stride)
+ 	    {
+ 	      x_offset -= x_stride;
+ 	      x_offset2++;
+ 	    }
+ 	  x_offset += x_offset2 * x_stride * x_len;
+ 	}
+ 
+       int idx_j;
+ 
+       double tmp_max = octave_NaN;
+ 
+       for (idx_j = 0; idx_j < x_len; idx_j++)
+ 	{
+ 	  tmp_max = elem (idx_j * x_stride + x_offset);
+ 	  
+ 	  if (! octave_is_NaN_or_NA (tmp_max))
+ 	    break;
+ 	}
+ 
+       for (int j = idx_j+1; j < x_len; j++)
+ 	{
+ 	  double tmp = elem (j * x_stride + x_offset);
+ 
+ 	  if (octave_is_NaN_or_NA (tmp))
+ 	    continue;
+ 	  else if (tmp > tmp_max)
+ 	    {
+ 	      idx_j = j;
+ 	      tmp_max = tmp;
+ 	    }
+ 	}
+ 
+       result.elem (i) = tmp_max;
+       idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j;
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ NDArray::min (int dim) const
+ {
+   ArrayN<int> dummy_idx;
+   return min (dummy_idx, dim);
+ }
+ 
+ NDArray
+ NDArray::min (ArrayN<int>& idx_arg, int dim) const
+ {
+   dim_vector dv = dims ();
+   dim_vector dr = dims ();
+ 
+   if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+     return NDArray ();
+   
+   dr(dim) = 1;
+ 
+   NDArray result (dr);
+   idx_arg.resize (dr);
+ 
+   int x_stride = 1;
+   int x_len = dv(dim);
+   for (int i = 0; i < dim; i++)
+     x_stride *= dv(i);
+ 
+   for (int i = 0; i < dr.numel (); i++)
+     {
+       int x_offset;
+       if (x_stride == 1)
+ 	x_offset = i * x_len;
+       else
+ 	{
+ 	  int x_offset2 = 0;
+ 	  x_offset = i;
+ 	  while (x_offset >= x_stride)
+ 	    {
+ 	      x_offset -= x_stride;
+ 	      x_offset2++;
+ 	    }
+ 	  x_offset += x_offset2 * x_stride * x_len;
+ 	}
+ 
+       int idx_j;
+ 
+       double tmp_min = octave_NaN;
+ 
+       for (idx_j = 0; idx_j < x_len; idx_j++)
+ 	{
+ 	  tmp_min = elem (idx_j * x_stride + x_offset);
+ 	  
+ 	  if (! octave_is_NaN_or_NA (tmp_min))
+ 	    break;
+ 	}
+ 
+       for (int j = idx_j+1; j < x_len; j++)
+ 	{
+ 	  double tmp = elem (j * x_stride + x_offset);
+ 
+ 	  if (octave_is_NaN_or_NA (tmp))
+ 	    continue;
+ 	  else if (tmp < tmp_min)
+ 	    {
+ 	      idx_j = j;
+ 	      tmp_min = tmp;
+ 	    }
+ 	}
+ 
+       result.elem (i) = tmp_min;
+       idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j;
+     }
+ 
+   return result;
+ }
+ 
  int
  NDArray::cat (const NDArray& ra_arg, int dim, int iidx, int move)
  {
***************
*** 776,781 ****
--- 926,1066 ----
    return is;
  }
  
+ // XXX FIXME XXX -- it would be nice to share code among the min/max
+ // functions below.
+ 
+ #define EMPTY_RETURN_CHECK(T) \
+   if (nel == 0)	\
+     return T (dv);
+ 
+ NDArray
+ min (double d, const NDArray& m)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (d, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ min (const NDArray& m, double d)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (d, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ min (const NDArray& a, const NDArray& b)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   if (dv != b.dims ())
+     {
+       (*current_liboctave_error_handler)
+ 	("two-arg min expecting args of same size");
+       return NDArray ();
+     }
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmin (a (i), b (i));
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ max (double d, const NDArray& m)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (d, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ max (const NDArray& m, double d)
+ {
+   dim_vector dv = m.dims ();
+   int nel = dv.numel ();
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (d, m (i));
+     }
+ 
+   return result;
+ }
+ 
+ NDArray
+ max (const NDArray& a, const NDArray& b)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   if (dv != b.dims ())
+     {
+       (*current_liboctave_error_handler)
+ 	("two-arg max expecting args of same size");
+       return NDArray ();
+     }
+ 
+   EMPTY_RETURN_CHECK (NDArray);
+ 
+   NDArray result (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       result (i) = xmax (a (i), b (i));
+     }
+ 
+   return result;
+ }
+ 
  NDS_CMP_OPS(NDArray, , double, )
  NDS_BOOL_OPS(NDArray, double, 0.0)
  
*** ./liboctave/dNDArray.h.orig	2004-03-26 18:17:38.000000000 +0100
--- ./liboctave/dNDArray.h	2004-03-26 18:26:17.000000000 +0100
***************
*** 87,93 ****
    NDArray sum (int dim = -1) const;  
    NDArray sumsq (int dim = -1) const;
    int cat (const NDArray& ra_arg, int dim, int iidx, int move);
!      
    NDArray abs (void) const;
  
    ComplexNDArray fourier (int dim = 1) const;
--- 87,98 ----
    NDArray sum (int dim = -1) const;  
    NDArray sumsq (int dim = -1) const;
    int cat (const NDArray& ra_arg, int dim, int iidx, int move);
! 
!   NDArray max (int dim = 0) const;
!   NDArray max (ArrayN<int>& index, int dim = 0) const;
!   NDArray min (int dim = 0) const;
!   NDArray min (ArrayN<int>& index, int dim = 0) const;
!   
    NDArray abs (void) const;
  
    ComplexNDArray fourier (int dim = 1) const;
***************
*** 125,130 ****
--- 130,143 ----
    NDArray (double *d, const dim_vector& dv) : MArrayN<double> (d, dv) { }
  };
  
+ extern NDArray min (double d, const NDArray& m);
+ extern NDArray min (const NDArray& m, double d);
+ extern NDArray min (const NDArray& a, const NDArray& b);
+ 
+ extern NDArray max (double d, const NDArray& m);
+ extern NDArray max (const NDArray& m, double d);
+ extern NDArray max (const NDArray& a, const NDArray& b);
+ 
  NDS_CMP_OP_DECLS (NDArray, double)
  NDS_BOOL_OP_DECLS (NDArray, double)
  
*** ./liboctave/ChangeLog.orig	2004-03-27 01:44:39.000000000 +0100
--- ./liboctave/ChangeLog	2004-03-29 15:35:22.000000000 +0200
***************
*** 1,3 ****
--- 1,23 ----
+ 2004-03-29  David Bateman <address@hidden>
+ 
+ 	* lo-specfun.cc (besselj, bessely, besseli, besselk, besselh1, 
+ 	besselh2, airy, biry, betainc, gammainc): Convert for NDArrays.
+ 	* lo-specfun.h (besselj, bessely, besseli, besselk, besselh1, 
+ 	besselh2, airy, biry, betainc, gammainc): Declarations
+ 	* lo-specfun.cc (SN_BESSEL, NS_BESSEL, NN_BESSEL): New macros.
+ 	* lo-specfun.cc (do_bessel): Convert for NDArrays.
+ 	
+ 2004-03-26  David Bateman <address@hidden>
+ 
+ 	* dNDArray.cc (NDArray::min, NDArray::max, min, max):
+ 	Min and Max functions for NDArrays.
+ 	* dNDArray.h (NDArray::min, NDArray::max, min, max): Decl.
+ 
+ 	* CNDArray.cc (ComplexNDArray::min, ComplexNDArray::max, min, max):
+ 	Min and Max functions for ComplexNDArrays.
+ 	* CNDArray.h (ComplexNDArray::min, ComplexNDArray::max, min, max): 
+ 	Decl.
+ 	
  2004-03-11  John W. Eaton  <address@hidden>
  
  	* so-array.cc (SND_CMP_OP, NDS_CMP_OP, NDND_CMP_OP):
*** ./liboctave/lo-specfun.h.orig	2004-03-29 13:13:38.000000000 +0200
--- ./liboctave/lo-specfun.h	2004-03-29 14:07:14.000000000 +0200
***************
*** 24,33 ****
--- 24,36 ----
  #define octave_liboctave_specfun_h 1
  
  #include "oct-cmplx.h"
+ #include "ArrayN.h"
  
  template <class T> class Array2;
  class Matrix;
  class ComplexMatrix;
+ class NDArray;
+ class ComplexNDArray;
  class RowVector;
  class ComplexColumnVector;
  class Range;
***************
*** 145,150 ****
--- 148,225 ----
  besselh2 (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
  	  Array2<int>& ierr);
  
+ extern ComplexNDArray
+ besselj (double alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ bessely (double alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besseli (double alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselk (double alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh1 (double alpha, const ComplexNDArray& x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh2 (double alpha, const ComplexNDArray& x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselj (const NDArray& alpha, Complex x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ bessely (const NDArray& alpha, Complex x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besseli (const NDArray& alpha, Complex x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselk (const NDArray& alpha, Complex x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh1 (const NDArray& alpha, Complex x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh2 (const NDArray& alpha, Complex x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselj (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ bessely (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besseli (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselk (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	 ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh1 (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ besselh2 (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+ 	  ArrayN<int>& ierr);
+ 
  extern ComplexMatrix
  besselj (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
  	 Array2<int>& ierr);
***************
*** 178,198 ****
--- 253,292 ----
  extern ComplexMatrix
  biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr);
  
+ extern ComplexNDArray
+ airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr);
+ 
+ extern ComplexNDArray
+ biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr);
+ 
  extern double betainc (double x, double a, double b);
  extern Matrix betainc (double x, double a, const Matrix& b);
  extern Matrix betainc (double x, const Matrix& a, double b);
  extern Matrix betainc (double x, const Matrix& a, const Matrix& b);
  
+ extern NDArray betainc (double x, double a, const NDArray& b);
+ extern NDArray betainc (double x, const NDArray& a, double b);
+ extern NDArray betainc (double x, const NDArray& a, const NDArray& b);
+ 
  extern Matrix betainc (const Matrix& x, double a, double b);
  extern Matrix betainc (const Matrix& x, double a, const Matrix& b);
  extern Matrix betainc (const Matrix& x, const Matrix& a, double b);
  extern Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b);
  
+ extern NDArray betainc (const NDArray& x, double a, double b);
+ extern NDArray betainc (const NDArray& x, double a, const NDArray& b);
+ extern NDArray betainc (const NDArray& x, const NDArray& a, double b);
+ extern NDArray betainc (const NDArray& x, const NDArray& a, const NDArray& b);
+ 
  extern double gammainc (double x, double a, bool& err);
  extern Matrix gammainc (double x, const Matrix& a);
  extern Matrix gammainc (const Matrix& x, double a);
  extern Matrix gammainc (const Matrix& x, const Matrix& a);
  
+ extern NDArray gammainc (double x, const NDArray& a);
+ extern NDArray gammainc (const NDArray& x, double a);
+ extern NDArray gammainc (const NDArray& x, const NDArray& a);
+ 
  inline double gammainc (double x, double a)
  {
    bool err;
*** ./liboctave/lo-specfun.cc.orig	2004-03-29 13:13:42.000000000 +0200
--- ./liboctave/lo-specfun.cc	2004-03-29 14:17:46.000000000 +0200
***************
*** 29,34 ****
--- 29,36 ----
  #include "CMatrix.h"
  #include "dRowVector.h"
  #include "dMatrix.h"
+ #include "dNDArray.h"
+ #include "CNDArray.h"
  #include "f77-fcn.h"
  #include "lo-error.h"
  #include "lo-ieee.h"
***************
*** 607,612 ****
--- 609,670 ----
    return retval;
  }
  
+ static inline ComplexNDArray
+ do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x,
+ 	   bool scaled, ArrayN<int>& ierr)
+ {
+   dim_vector dv = x.dims ();
+   int nel = dv.numel ();
+   ComplexNDArray retval (dv);
+ 
+   ierr.resize (dv);
+ 
+   for (int i = 0; i < nel; i++)
+       retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
+ 
+   return retval;
+ }
+ 
+ static inline ComplexNDArray
+ do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x,
+ 	   bool scaled, ArrayN<int>& ierr)
+ {
+   dim_vector dv = alpha.dims ();
+   int nel = dv.numel ();
+   ComplexNDArray retval (dv);
+ 
+   ierr.resize (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
+ 
+   return retval;
+ }
+ 
+ static inline ComplexNDArray
+ do_bessel (fptr f, const char *fn, const NDArray& alpha,
+ 	   const ComplexNDArray& x, bool scaled, ArrayN<int>& ierr)
+ {
+   dim_vector dv = x.dims ();
+   ComplexNDArray retval;
+ 
+   if (dv == alpha.dims ())
+     {
+       int nel = dv.numel ();
+ 
+       retval.resize (dv);
+       ierr.resize (dv);
+ 
+       for (int i = 0; i < nel; i++)
+ 	retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
+     }
+   else
+     (*current_liboctave_error_handler)
+       ("%s: the sizes of alpha and x must conform", fn);
+ 
+   return retval;
+ }
+ 
  static inline ComplexMatrix
  do_bessel (fptr f, const char *, const RowVector& alpha,
  	   const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
***************
*** 656,661 ****
--- 714,743 ----
      return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
    }
  
+ #define SN_BESSEL(name, fcn) \
+   ComplexNDArray \
+   name (double alpha, const ComplexNDArray& x, bool scaled, \
+ 	ArrayN<int>& ierr) \
+   { \
+     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+   }
+ 
+ #define NS_BESSEL(name, fcn) \
+   ComplexNDArray \
+   name (const NDArray& alpha, const Complex& x, bool scaled, \
+ 	ArrayN<int>& ierr) \
+   { \
+     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+   }
+ 
+ #define NN_BESSEL(name, fcn) \
+   ComplexNDArray \
+   name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
+ 	ArrayN<int>& ierr) \
+   { \
+     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+   }
+ 
  #define RC_BESSEL(name, fcn) \
    ComplexMatrix \
    name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
***************
*** 669,674 ****
--- 751,759 ----
    SM_BESSEL (name, fcn) \
    MS_BESSEL (name, fcn) \
    MM_BESSEL (name, fcn) \
+   SN_BESSEL (name, fcn) \
+   NS_BESSEL (name, fcn) \
+   NN_BESSEL (name, fcn) \
    RC_BESSEL (name, fcn)
  
  ALL_BESSEL (besselj, zbesj)
***************
*** 778,783 ****
--- 863,898 ----
    return retval;
  }
  
+ ComplexNDArray
+ airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
+ {
+   dim_vector dv = z.dims ();
+   int nel = dv.numel ();
+   ComplexNDArray retval (dv);
+ 
+   ierr.resize (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval (i) = airy (z(i), deriv, scaled, ierr(i));
+ 
+   return retval;
+ }
+ 
+ ComplexNDArray
+ biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
+ {
+   dim_vector dv = z.dims ();
+   int nel = dv.numel ();
+   ComplexNDArray retval (dv);
+ 
+   ierr.resize (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval (i) = biry (z(i), deriv, scaled, ierr(i));
+ 
+   return retval;
+ }
+ 
  static void
  gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
  			     int c3)
***************
*** 787,792 ****
--- 902,922 ----
       r1, c1, r2, c2, r3, c3);
  }
  
+ static dim_vector null_dims (0);
+ 
+ static void
+ gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2,
+ 			     const dim_vector& d3)
+ {
+   std::string d1_str = d1.str ();
+   std::string d2_str = d2.str ();
+   std::string d3_str = d3.str ();
+ 
+   (*current_liboctave_error_handler)
+   ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
+    d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
+ }
+ 
  double
  betainc (double x, double a, double b)
  {
***************
*** 850,855 ****
--- 980,1035 ----
    return retval;
  }
  
+ NDArray
+ betainc (double x, double a, const NDArray& b)
+ {
+   dim_vector dv = b.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval (i) = betainc (x, a, b(i));
+ 
+   return retval;
+ }
+ 
+ NDArray
+ betainc (double x, const NDArray& a, double b)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval (i) = betainc (x, a(i), b);
+ 
+   return retval;
+ }
+ 
+ NDArray
+ betainc (double x, const NDArray& a, const NDArray& b)
+ {
+   NDArray retval;
+   dim_vector dv = a.dims ();
+ 
+   if (dv == b.dims ())
+     {
+       int nel = dv.numel ();
+ 
+       retval.resize (dv);
+ 
+       for (int i = 0; i < nel; i++)
+ 	retval (i) = betainc (x, a(i), b(i));
+     }
+   else
+     gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ());
+   
+   return retval;
+ }
+ 
+ 
  Matrix
  betainc (const Matrix& x, double a, double b)
  {
***************
*** 943,948 ****
--- 1123,1205 ----
    return retval;
  }
  
+ NDArray
+ betainc (const NDArray& x, double a, double b)
+ {
+   dim_vector dv = x.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     retval (i) = betainc (x(i), a, b);
+ 
+   return retval;
+ }
+ 
+ NDArray
+ betainc (const NDArray& x, double a, const NDArray& b)
+ {
+   NDArray retval;
+   dim_vector dv = x.dims ();
+ 
+   if (dv == b.dims ())
+     {
+       int nel = dv.numel ();
+ 
+       retval.resize (dv);
+ 
+       for (int i = 0; i < nel; i++)
+ 	retval (i) = betainc (x(i), a, b(i));
+     }
+   else
+     gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ());
+   
+   return retval;
+ }
+ 
+ NDArray
+ betainc (const NDArray& x, const NDArray& a, double b)
+ {
+   NDArray retval;
+   dim_vector dv = x.dims ();
+ 
+   if (dv == a.dims ())
+     {
+       int nel = dv.numel ();
+ 
+       retval.resize (dv);
+ 
+       for (int i = 0; i < nel; i++)
+ 	retval (i) = betainc (x(i), a(i), b);
+     }
+   else
+     gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0));
+   
+   return retval;
+ }
+ 
+ NDArray
+ betainc (const NDArray& x, const NDArray& a, const NDArray& b)
+ {
+   NDArray retval;
+   dim_vector dv = x.dims ();
+ 
+   if (dv == a.dims () && dv == b.dims ())
+     {
+       int nel = dv.numel ();
+ 
+       retval.resize (dv);
+ 
+       for (int i = 0; i < nel; i++)
+ 	retval (i) = betainc (x(i), a(i), b(i));
+     }
+   else
+     gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
+ 
+   return retval;
+ }
+ 
  // XXX FIXME XXX -- there is still room for improvement here...
  
  double
***************
*** 1058,1063 ****
--- 1315,1412 ----
    return retval;
  }
  
+ NDArray
+ gammainc (double x, const NDArray& a)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval;
+   NDArray result (dv);
+ 
+   bool err;
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       result (i) = gammainc (x, a(i), err);
+ 
+       if (err)
+ 	goto done;
+     }
+ 
+   retval = result;
+ 
+  done:
+ 
+   return retval;
+ }
+ 
+ NDArray
+ gammainc (const NDArray& x, double a)
+ {
+   dim_vector dv = x.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval;
+   NDArray result (dv);
+ 
+   bool err;
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       result (i) = gammainc (x(i), a, err);
+ 
+       if (err)
+ 	goto done;
+     }
+ 
+   retval = result;
+ 
+  done:
+ 
+   return retval;
+ }
+ 
+ NDArray
+ gammainc (const NDArray& x, const NDArray& a)
+ {
+   dim_vector dv = x.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval;
+   NDArray result;
+ 
+   if (dv == a.dims ())
+     {
+       result.resize (dv);
+ 
+       bool err;
+ 
+       for (int i = 0; i < nel; i++)
+ 	{
+ 	  result (i) = gammainc (x(i), a(i), err);
+ 	  
+ 	  if (err)
+ 	    goto done;
+ 	}
+ 
+       retval = result;
+     }
+   else
+     {
+       std::string x_str = dv.str ();
+       std::string a_str = a.dims ().str ();
+ 
+       (*current_liboctave_error_handler)
+ 	("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
+ 	 x_str.c_str (), a_str. c_str ());
+     }
+ 
+  done:
+ 
+   return retval;
+ }
+ 
  /*
  ;;; Local Variables: ***
  ;;; mode: C++ ***
*** ./scripts/control/base/bode.m.orig	2004-02-16 20:57:22.000000000 +0100
--- ./scripts/control/base/bode.m	2004-03-25 21:44:58.000000000 +0100
***************
*** 67,73 ****
  ##
  ## @strong{Example}
  ## @example
! ## bode(sys,[],"y_3",list("u_1","u_4");
  ## @end example
  ## @end table
  ## @strong{Outputs}
--- 67,73 ----
  ##
  ## @strong{Example}
  ## @example
! ## bode(sys,[],"y_3", @{"u_1","u_4"@});
  ## @end example
  ## @end table
  ## @strong{Outputs}
*** ./scripts/control/base/lqg.m.orig	2004-02-16 20:57:22.000000000 +0100
--- ./scripts/control/base/lqg.m	2004-03-25 18:43:57.000000000 +0100
***************
*** 40,46 ****
  ## @itemx  r
  ## state, control weighting respectively.  Control ARE is
  ## @item  in_idx
! ## names or indices of controlled inputs (see @code{sysidx}, @code{listidx})
  ##
  ## default: last dim(R) inputs are assumed to be controlled inputs, all
  ## others are assumed to be noise inputs.
--- 40,46 ----
  ## @itemx  r
  ## state, control weighting respectively.  Control ARE is
  ## @item  in_idx
! ## names or indices of controlled inputs (see @code{sysidx}, @code{cellidx})
  ##
  ## default: last dim(R) inputs are assumed to be controlled inputs, all
  ## others are assumed to be noise inputs.
*** ./scripts/control/system/ss2sys.m.orig	2004-02-16 20:57:22.000000000 +0100
--- ./scripts/control/system/ss2sys.m	2004-03-25 21:46:25.000000000 +0100
***************
*** 46,63 ****
  ## see below for system partitioning
  ##
  ## @item  stname
! ## list of strings of state signal names
  ##
  ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
  ##                     @code{xd_n} for discrete states
  ##
  ## @item inname
! ## list of strings of input signal names
  ##
  ## default (@var{inname} = [] on input): @code{u_n}
  ##
  ## @item outname
! ## list of strings of input signal names
  ##
  ## default (@var{outname} = [] on input): @code{y_n}
  ##
--- 46,63 ----
  ## see below for system partitioning
  ##
  ## @item  stname
! ## cell array of strings of state signal names
  ##
  ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
  ##                     @code{xd_n} for discrete states
  ##
  ## @item inname
! ## cell array of strings of input signal names
  ##
  ## default (@var{inname} = [] on input): @code{u_n}
  ##
  ## @item outname
! ## cell array of strings of input signal names
  ##
  ## default (@var{outname} = [] on input): @code{y_n}
  ##
***************
*** 140,146 ****
  ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
  ## octave:2> b = [0 0 ; 0 1 ; 1 0];
  ## octave:3> c = eye(3);
! ## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules"));
  ## octave:5> sysout(sys);
  ## Input(s)
  ##         1: u_1
--- 140,146 ----
  ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
  ## octave:2> b = [0 0 ; 0 1 ; 1 0];
  ## octave:3> c = eye(3);
! ## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@});
  ## octave:5> sysout(sys);
  ## Input(s)
  ##         1: u_1
*** ./scripts/control/system/cellidx.m.orig	2004-02-17 03:34:33.000000000 +0100
--- ./scripts/control/system/cellidx.m	2004-03-25 18:43:01.000000000 +0100
***************
*** 17,23 ****
  ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
  
  ## -*- texinfo -*-
! ## @deftypefn {Function File} address@hidden, @var{errmsg}] =} listidx (@var{listvar}, @var{strlist})
  ## Return indices of string entries in @var{listvar} that match strings
  ## in @var{strlist}.
  ##
--- 17,23 ----
  ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
  
  ## -*- texinfo -*-
! ## @deftypefn {Function File} address@hidden, @var{errmsg}] =} cellidx (@var{listvar}, @var{strlist})
  ## Return indices of string entries in @var{listvar} that match strings
  ## in @var{strlist}.
  ##
***************
*** 29,35 ****
  ##
  ## If @var{strlist} contains a string not in @var{listvar}, then
  ## an error message is returned in @var{errmsg}.  If only one output
! ## argument is requested, then @var{listidx} prints @var{errmsg} to the
  ## screen and exits with an error.
  ## @end deftypefn
  
--- 29,35 ----
  ##
  ## If @var{strlist} contains a string not in @var{listvar}, then
  ## an error message is returned in @var{errmsg}.  If only one output
! ## argument is requested, then @var{cellidx} prints @var{errmsg} to the
  ## screen and exits with an error.
  ## @end deftypefn
  
*** ./scripts/control/system/dmr2d.m.orig	2004-02-16 20:57:22.000000000 +0100
--- ./scripts/control/system/dmr2d.m	2004-03-25 18:43:20.000000000 +0100
***************
*** 29,35 ****
  ## @code{dmr2d} exits with an error if @var{sys} is not discrete
  ## @item   idx
  ## indices or names of states with sampling time 
! ## @code{sysgettsam(@var{sys})} (may be empty); see @code{listidx}
  ## @item   sprefix
  ## list of string prefixes of states with sampling time
  ## @code{sysgettsam(@var{sys})} (may be empty)
--- 29,35 ----
  ## @code{dmr2d} exits with an error if @var{sys} is not discrete
  ## @item   idx
  ## indices or names of states with sampling time 
! ## @code{sysgettsam(@var{sys})} (may be empty); see @code{cellidx}
  ## @item   sprefix
  ## list of string prefixes of states with sampling time
  ## @code{sysgettsam(@var{sys})} (may be empty)
*** ./scripts/control/system/ss.m.orig	2004-02-17 03:34:33.000000000 +0100
--- ./scripts/control/system/ss.m	2004-03-25 21:46:39.000000000 +0100
***************
*** 46,63 ****
  ## see below for system partitioning
  ##
  ## @item  stname
! ## list of strings of state signal names
  ##
  ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
  ##                     @code{xd_n} for discrete states
  ##
  ## @item inname
! ## list of strings of input signal names
  ##
  ## default (@var{inname} = [] on input): @code{u_n}
  ##
  ## @item outname
! ## list of strings of input signal names
  ##
  ## default (@var{outname} = [] on input): @code{y_n}
  ##
--- 46,63 ----
  ## see below for system partitioning
  ##
  ## @item  stname
! ## cell array of strings of state signal names
  ##
  ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
  ##                     @code{xd_n} for discrete states
  ##
  ## @item inname
! ## cell array of strings of input signal names
  ##
  ## default (@var{inname} = [] on input): @code{u_n}
  ##
  ## @item outname
! ## cell array of strings of input signal names
  ##
  ## default (@var{outname} = [] on input): @code{y_n}
  ##
***************
*** 140,146 ****
  ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
  ## octave:2> b = [0 0 ; 0 1 ; 1 0];
  ## octave:3> c = eye(3);
! ## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules"));
  ## octave:5> sysout(sys);
  ## Input(s)
  ##         1: u_1
--- 140,146 ----
  ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
  ## octave:2> b = [0 0 ; 0 1 ; 1 0];
  ## octave:3> c = eye(3);
! ## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@});
  ## octave:5> sysout(sys);
  ## Input(s)
  ##         1: u_1
*** ./scripts/control/system/sysidx.m.orig	2002-08-09 20:58:14.000000000 +0200
--- ./scripts/control/system/sysidx.m	2004-03-25 18:44:41.000000000 +0100
***************
*** 38,50 ****
    endif
  
    ## extract correct set of signal names values
!   [idxvec, msg] = listidx (list ("in", "out", "st", "yd"), sigtype);
    if (msg)
      error ("invalid sigtype = %s", sigtype);
    endif
  
    syssiglist = sysgetsignals (sys, sigtype);
!   [idxvec, msg] = listidx (syssiglist, signamelist);
    if (length (msg))
      error ("sysidx (sigtype = %s): %s", sigtype,
  	   strrep (msg, "strlist", "signamelist"));
--- 38,50 ----
    endif
  
    ## extract correct set of signal names values
!   [idxvec, msg] = cellidx ({"in", "out", "st", "yd"}, sigtype);
    if (msg)
      error ("invalid sigtype = %s", sigtype);
    endif
  
    syssiglist = sysgetsignals (sys, sigtype);
!   [idxvec, msg] = cellidx (syssiglist, signamelist);
    if (length (msg))
      error ("sysidx (sigtype = %s): %s", sigtype,
  	   strrep (msg, "strlist", "signamelist"));
*** ./scripts/control/system/sysprune.m.orig	2004-02-16 20:57:22.000000000 +0100
--- ./scripts/control/system/sysprune.m	2004-03-25 21:58:22.000000000 +0100
***************
*** 33,39 ****
  ##
  ## @example
  ## retsys = sysprune(Asys,[1:3,4],"u_1");
! ## retsys = sysprune(Asys,list("tx","ty","tz"), 4);
  ## @end example
  ##
  ## @end table
--- 33,39 ----
  ##
  ## @example
  ## retsys = sysprune(Asys,[1:3,4],"u_1");
! ## retsys = sysprune(Asys,@{"tx","ty","tz"@}, 4);
  ## @end example
  ##
  ## @end table
*** ./scripts/ChangeLog.orig	2004-03-12 20:26:56.000000000 +0100
--- ./scripts/ChangeLog	2004-03-25 22:01:06.000000000 +0100
***************
*** 1,3 ****
--- 1,10 ----
+ 2004-03-25  David Bateman <address@hidden>
+ 	* control/base/bode.m, control/base/lqg.m, control/system/ss2sys.m,
+ 	control/system/cellidx.m, control/system/dmr2d.m control/system/ss.m,
+ 	control/system/sysprune.m: Doc update for usage of cell arrays.
+ 
+ 	* control/system/sysidx.m: Use cellidx and not listidx.
+ 	
  2004-03-12  Stefan van der Walt  <address@hidden>
  
  	* image/imshow.m: Accept "truesize" argument.
*** ./src/DLD-FUNCTIONS/betainc.cc.orig	2004-03-29 15:27:02.000000000 +0200
--- ./src/DLD-FUNCTIONS/betainc.cc	2004-03-29 15:26:52.000000000 +0200
***************
*** 88,94 ****
  		    }
  		  else
  		    {
! 		      Matrix b = b_arg.matrix_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
--- 88,94 ----
  		    }
  		  else
  		    {
! 		      NDArray b = b_arg.array_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
***************
*** 97,103 ****
  	    }
  	  else
  	    {
! 	      Matrix a = a_arg.matrix_value ();
  
  	      if (! error_state)
  		{
--- 97,103 ----
  	    }
  	  else
  	    {
! 	      NDArray a = a_arg.array_value ();
  
  	      if (! error_state)
  		{
***************
*** 110,116 ****
  		    }
  		  else
  		    {
! 		      Matrix b = b_arg.matrix_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
--- 110,116 ----
  		    }
  		  else
  		    {
! 		      NDArray b = b_arg.array_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
***************
*** 120,126 ****
  	}
        else
  	{
! 	  Matrix x = x_arg.matrix_value ();
  
  	  if (a_arg.is_scalar_type ())
  	    {
--- 120,126 ----
  	}
        else
  	{
! 	  NDArray x = x_arg.array_value ();
  
  	  if (a_arg.is_scalar_type ())
  	    {
***************
*** 137,143 ****
  		    }
  		  else
  		    {
! 		      Matrix b = b_arg.matrix_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
--- 137,143 ----
  		    }
  		  else
  		    {
! 		      NDArray b = b_arg.array_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
***************
*** 146,152 ****
  	    }
  	  else
  	    {
! 	      Matrix a = a_arg.matrix_value ();
  
  	      if (! error_state)
  		{
--- 146,152 ----
  	    }
  	  else
  	    {
! 	      NDArray a = a_arg.array_value ();
  
  	      if (! error_state)
  		{
***************
*** 159,165 ****
  		    }
  		  else
  		    {
! 		      Matrix b = b_arg.matrix_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
--- 159,165 ----
  		    }
  		  else
  		    {
! 		      NDArray b = b_arg.array_value ();
  
  		      if (! error_state)
  			retval = betainc (x, a, b);
*** ./src/DLD-FUNCTIONS/filter.cc.orig	2004-03-26 12:21:58.000000000 +0100
--- ./src/DLD-FUNCTIONS/filter.cc	2004-03-26 15:49:32.000000000 +0100
***************
*** 39,170 ****
  #include "oct-obj.h"
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&);
  
! extern MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&);
  #endif
  
  template <class T>
! MArray<T>
! filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si)
  {
!   MArray<T> y;
  
    int a_len  = a.length ();
    int b_len  = b.length ();
-   int x_len  = x.length ();
- 
-   int si_len = si.length ();
  
    int ab_len = a_len > b_len ? a_len : b_len;
  
    b.resize (ab_len, 0.0);
  
!   if (si.length () != ab_len - 1)
      {
!       error ("filter: si must be a vector of length max (length (a), length (b)) - 1");
        return y;
      }
  
!   T norm = a (0);
  
!   if (norm == 0.0)
      {
!       error ("filter: the first element of a must be non-zero");
        return y;
      }
  
!   y.resize (x_len, 0.0);
  
!   if (norm != 1.0)
!     b = b / norm;
  
!   if (a_len > 1)
      {
!       a.resize (ab_len, 0.0);
  
!       if (norm != 1.0)
! 	a = a / norm;
  
!       for (int i = 0; i < x_len; i++)
  	{
! 	  y (i) = si (0) + b (0) * x (i);
  
! 	  if (si_len > 1)
  	    {
! 	      for (int j = 0; j < si_len - 1; j++)
  		{
! 		  OCTAVE_QUIT;
  
! 		  si (j) = si (j+1) - a (j+1) * y (i)
! 		    + b (j+1) * x (i);
  		}
! 
! 	      si (si_len-1) = b (si_len) * x (i)
! 		- a (si_len) * y (i);
  	    }
- 	  else
- 	    si (0) = b (si_len) * x (i)
- 	      - a (si_len) * y (i);
  	}
!     }
!   else if (si_len > 0)
!     {
!       for (int i = 0; i < x_len; i++)
  	{
! 	  y (i) = si (0) + b (0) * x (i);
! 
! 	  if (si_len > 1)
  	    {
! 	      for (int j = 0; j < si_len - 1; j++)
  		{
! 		  OCTAVE_QUIT;
  
! 		  si (j) = si (j+1) + b (j+1) * x (i);
  		}
! 
! 	      si (si_len-1) = b (si_len) * x (i);
  	    }
- 	  else
- 	    si (0) = b (1) * x (i);
  	}
      }
!   else
!     y = b (0) * x;
! 
    return y;
  }
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&,
! 	MArray<double>&);
! 
! extern MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
! 	MArray<Complex>&);
  #endif
  
  template <class T>
! MArray<T>
! filter (MArray<T>& b, MArray<T>& a, MArray<T>& x)
  {
    int a_len = a.length ();
    int b_len = b.length ();
  
    int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
!   MArray<T> si (si_len, T (0.0));
! 
!   return filter (b, a, x, si);
  }
  
  DEFUN_DLD (filter, args, nargout,
    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
  @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
  Return the solution to the following linear, time-invariant difference\n\
  equation:\n\
  @iftex\n\
--- 39,252 ----
  #include "oct-obj.h"
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
  
! extern MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
  #endif
  
  template <class T>
! MArrayN<T>
! filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, MArrayN<T>& si, 
! 	int dim = 0)
  {
!   MArrayN<T> y;
  
    int a_len  = a.length ();
    int b_len  = b.length ();
  
    int ab_len = a_len > b_len ? a_len : b_len;
  
    b.resize (ab_len, 0.0);
+   if (a_len > 1)
+     a.resize (ab_len, 0.0);
  
!   T norm = a (0);
! 
!   if (norm == 0.0)
      {
!       error ("filter: the first element of a must be non-zero");
        return y;
      }
  
!   dim_vector x_dims = x.dims ();
!   if ((dim < 0) || (dim > x_dims.length ()))
!     {
!       error ("filter: filtering over invalid dimension");
!       return y;
!     }
  
!   int x_len = x_dims (dim);
! 
!   dim_vector si_dims = si.dims ();
!   int si_len = si_dims (0);
! 
!   if (si_len != ab_len - 1)
      {
!       error ("filter: first dimension of si must be of length max (length (a), length (b)) - 1");
        return y;
      }
  
!   if (si_dims.length () != x_dims.length ())
!     {
!       error ("filter: dimensionality of si and x must agree");
!       return y;
!     }
  
!   int si_dim = 0;
!   for (int i = 0; i < x_dims.length (); i++)
!     {
!       if (i == dim)
! 	continue;
!       
!       if (si_dims (++si_dim) != x_dims (i))
! 	{
! 	  error ("filter: dimensionality of si and x must agree");
! 	  return y;
! 	}
!     }
  
!   if (norm != 1.0)
      {
!       a = a / norm;
!       b = b / norm;
!     }
  
!   if ((a_len <= 1) && (si_len <= 1))
!     return b(0) * x;
  
!   y.resize (x_dims, 0.0);
! 
!   int x_stride = 1;
!   for (int i = 0; i < dim; i++)
!     x_stride *= x_dims(i);
! 
!   int x_num = x_dims.numel () / x_len;
!   for (int num = 0; num < x_num; num++)
!     {
!       int x_offset;
!       if (x_stride == 1)
! 	x_offset = num * x_len;
!       else
  	{
! 	  int x_offset2 = 0;
! 	  x_offset = num;
! 	  while (x_offset >= x_stride)
! 	    {
! 	      x_offset -= x_stride;
! 	      x_offset2++;
! 	    }
! 	  x_offset += x_offset2 * x_stride * x_len;
! 	}
!       int si_offset = num * si_len;
  
!       if (a_len > 1)
! 	{
! 	  for (int i = 0; i < x_len; i++)
  	    {
! 	      int idx = i * x_stride + x_offset; 
! 	      y (idx) = si (si_offset) + b (0) * x (idx);
! 
! 	      if (si_len > 1)
  		{
! 		  for (int j = 0; j < si_len - 1; j++)
! 		    {
! 		      OCTAVE_QUIT;
! 
! 		      si (j + si_offset) =  si (j + 1 + si_offset) - 
! 			a (j+1) * y (idx) + b (j+1) * x (idx);
! 		    }
  
! 		  si (si_len - 1 + si_offset) = b (si_len) * x (idx)
! 		    - a (si_len) * y (idx);
  		}
! 	      else
! 		si (si_offset) = b (si_len) * x (idx)
! 		  - a (si_len) * y (idx);
  	    }
  	}
!       else if (si_len > 0)
  	{
! 	  for (int i = 0; i < x_len; i++)
  	    {
! 	      int idx = i * x_stride + x_offset; 
! 	      y (idx) = si (si_offset) + b (0) * x (idx);
! 
! 	      if (si_len > 1)
  		{
! 		  for (int j = 0; j < si_len - 1; j++)
! 		    {
! 		      OCTAVE_QUIT;
! 
! 		      si (j + si_offset) = si (j + 1 + si_offset) + 
! 			b (j+1) * x (idx);
! 		    }
  
! 		  si (si_len - 1 + si_offset) = b (si_len) * x (idx);
  		}
! 	      else
! 		si (si_offset) = b (1) * x (idx);
  	    }
  	}
      }
!   
    return y;
  }
  
  #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
! extern MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
! 	MArrayN<double>&, int dim);
! 
! extern MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
! 	MArrayN<Complex>&, int dim);
  #endif
  
  template <class T>
! MArrayN<T>
! filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, int dim = -1)
  {
+   dim_vector x_dims = x.dims ();
+ 
+   if (dim < 0)
+     {
+       // Find first non-singleton dimension
+       while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
+ 	dim++;
+   
+       // All dimensions singleton, pick first dimension
+       if (dim == x_dims.length ())
+ 	dim = 0;
+     }
+   else
+     if (dim < 0 || dim > x_dims.length ())
+       {
+ 	error ("filter: filtering over invalid dimension");
+ 	return MArrayN<T> ();
+       }
+ 
    int a_len = a.length ();
    int b_len = b.length ();
  
    int si_len = (a_len > b_len ? a_len : b_len) - 1;
+   dim_vector si_dims = x.dims ();
+   for (int i = dim; i > 0; i--)
+     si_dims (i) = si_dims (i-1);
+   si_dims (0) = si_len;
+   
+   MArrayN<T> si (si_dims, T (0.0));
  
!   return filter (b, a, x, si, dim);
  }
  
  DEFUN_DLD (filter, args, nargout,
    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
  @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
+ @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\
+ @deftypefnx {Loadable Function} address@hidden, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\
  Return the solution to the following linear, time-invariant difference\n\
  equation:\n\
  @iftex\n\
***************
*** 194,200 ****
   $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
  @end tex\n\
  @end iftex\n\
! An equivalent form of this equation is:\n\
  @iftex\n\
  @tex\n\
  $$\n\
--- 276,283 ----
   $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
  @end tex\n\
  @end iftex\n\
! over the first non-singleton dimension of @var{x} or over @var{dim} if\n\
! supplied. An equivalent form of this equation is:\n\
  @iftex\n\
  @tex\n\
  $$\n\
***************
*** 259,317 ****
  
    int nargin  = args.length ();
  
!   if (nargin < 3 || nargin > 4)
      {
        print_usage ("filter");
        return retval;
      }
  
!   const char *errmsg = "filter: arguments must be vectors";
  
!   bool x_is_row_vector = (args(2).rows () == 1);
  
!   bool si_is_row_vector = (nargin == 4 && args(3).rows () == 1);
  
    if (args(0).is_complex_type ()
        || args(1).is_complex_type ()
        || args(2).is_complex_type ()
!       || (nargin == 4 && args(3).is_complex_type ()))
      {
        ComplexColumnVector b (args(0).complex_vector_value ());
        ComplexColumnVector a (args(1).complex_vector_value ());
!       ComplexColumnVector x (args(2).complex_vector_value ());
  
        if (! error_state)
  	{
! 	  ComplexColumnVector si;
  
! 	  if (nargin == 3)
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      si.resize (si_len, 0.0);
  	    }
  	  else
! 	    si = ComplexColumnVector (args(3).complex_vector_value ());
  
  	  if (! error_state)
  	    {
! 	      ComplexColumnVector y (filter (b, a, x, si));
  
  	      if (nargout == 2)
! 		{
! 		  if (si_is_row_vector)
! 		    retval(1) = si.transpose ();
! 		  else
! 		    retval(1) = si;
! 		}
  
! 	      if (x_is_row_vector)
! 		retval(0) = y.transpose ();
! 	      else
! 		retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
--- 342,432 ----
  
    int nargin  = args.length ();
  
!   if (nargin < 3 || nargin > 5)
      {
        print_usage ("filter");
        return retval;
      }
  
!   const char *errmsg = "filter: arguments a and b must be vectors";
  
!   int dim;
!   dim_vector x_dims = args(2).dims ();
  
!   if (nargin == 5)
!     {
!       dim = args(4).nint_value() - 1;
!       if (dim < 0 || dim >= x_dims.length ())
! 	{
! 	  error ("filter: filtering over invalid dimension");
! 	  return retval;
! 	}
!     }
!   else
!     {
!       // Find first non-singleton dimension
!       dim = 0;
!       while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
! 	dim++;
!   
!       // All dimensions singleton, pick first dimension
!       if (dim == x_dims.length ())
! 	dim = 0;
!     }
  
    if (args(0).is_complex_type ()
        || args(1).is_complex_type ()
        || args(2).is_complex_type ()
!       || (nargin >= 4 && args(3).is_complex_type ()))
      {
        ComplexColumnVector b (args(0).complex_vector_value ());
        ComplexColumnVector a (args(1).complex_vector_value ());
! 
!       ComplexNDArray x (args(2).complex_array_value ());
  
        if (! error_state)
  	{
! 	  ComplexNDArray si;
  
! 	  if (nargin == 3 || args(3).is_empty ())
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      dim_vector si_dims = x.dims ();
! 	      for (int i = dim; i > 0; i--)
! 		si_dims (i) = si_dims (i-1);
! 	      si_dims (0) = si_len;
! 
! 	      si.resize (si_dims, 0.0);
  	    }
  	  else
! 	    {
! 	      dim_vector si_dims = args (3).dims ();
! 	      bool si_is_vector = true;
! 	      for (int i=0; i < si_dims.length (); i++)
! 		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
! 		  {
! 		    si_is_vector = false;
! 		    break;
! 		  }
! 
! 	      if (si_is_vector)
! 		si = ComplexNDArray (args(3).complex_vector_value ());
! 	      else
! 		si = args(3).complex_array_value ();
! 	    }
  
  	  if (! error_state)
  	    {
! 	      ComplexNDArray y (filter (b, a, x, si, dim));
  
  	      if (nargout == 2)
! 		retval(1) = si;
  
! 	      retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
***************
*** 323,362 ****
      {
        ColumnVector b (args(0).vector_value ());
        ColumnVector a (args(1).vector_value ());
!       ColumnVector x (args(2).vector_value ());
  
        if (! error_state)
  	{
! 	  ColumnVector si;
  
! 	  if (nargin == 3)
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      si.resize (si_len, 0.0);
  	    }
  	  else
! 	    si = ColumnVector (args(3).vector_value ());
  
  	  if (! error_state)
  	    {
! 	      ColumnVector y (filter (b, a, x, si));
  
  	      if (nargout == 2)
! 		{
! 		  if (si_is_row_vector)
! 		    retval(1) = si.transpose ();
! 		  else
! 		    retval(1) = si;
! 		}
  
! 	      if (x_is_row_vector)
! 		retval(0) = y.transpose ();
! 	      else
! 		retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
--- 438,489 ----
      {
        ColumnVector b (args(0).vector_value ());
        ColumnVector a (args(1).vector_value ());
! 
!       NDArray x (args(2).array_value ());
  
        if (! error_state)
  	{
! 	  NDArray si;
  
! 	  if (nargin == 3 || args(3).is_empty ())
  	    {
  	      int a_len = a.length ();
  	      int b_len = b.length ();
  
  	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
  
! 	      dim_vector si_dims = x.dims ();
! 	      for (int i = dim; i > 0; i--)
! 		si_dims (i) = si_dims (i-1);
! 	      si_dims (0) = si_len;
! 
! 	      si.resize (si_dims, 0.0);
  	    }
  	  else
! 	    {
! 	      dim_vector si_dims = args (3).dims ();
! 	      bool si_is_vector = true;
! 	      for (int i=0; i < si_dims.length (); i++)
! 		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
! 		  {
! 		    si_is_vector = false;
! 		    break;
! 		  }
! 
! 	      if (si_is_vector)
! 		si = NDArray (args(3).vector_value ());
! 	      else
! 		si = args(3).array_value ();
! 	    }
  
  	  if (! error_state)
  	    {
! 	      NDArray y (filter (b, a, x, si, dim));
  
  	      if (nargout == 2)
! 		retval(1) = si;
  
! 	      retval(0) = y;
  	    }
  	  else
  	    error (errmsg);
***************
*** 368,386 ****
    return retval;
  }
  
! template MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&,
! 	MArray<double>&);
! 
! template MArray<double>
! filter (MArray<double>&, MArray<double>&, MArray<double>&);
! 
! template MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
! 	MArray<Complex>&);
  
! template MArray<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&);
  
  /*
  ;;; Local Variables: ***
--- 495,513 ----
    return retval;
  }
  
! template MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
! 	MArrayN<double>&, int dim);
! 
! template MArrayN<double>
! filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
! 
! template MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
! 	MArrayN<Complex>&, int dim);
  
! template MArrayN<Complex>
! filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
  
  /*
  ;;; Local Variables: ***
*** ./src/DLD-FUNCTIONS/minmax.cc.orig	2004-03-26 18:17:53.000000000 +0100
--- ./src/DLD-FUNCTIONS/minmax.cc	2004-03-27 01:37:06.000000000 +0100
***************
*** 28,35 ****
  
  #include "lo-ieee.h"
  #include "lo-mappers.h"
! #include "dMatrix.h"
! #include "CMatrix.h"
  #include "quit.h"
  
  #include "defun-dld.h"
--- 28,35 ----
  
  #include "lo-ieee.h"
  #include "lo-mappers.h"
! #include "dNDArray.h"
! #include "CNDArray.h"
  #include "quit.h"
  
  #include "defun-dld.h"
***************
*** 37,49 ****
  #include "gripes.h"
  #include "oct-obj.h"
  
  #define MINMAX_BODY(FCN) \
   \
    octave_value_list retval;  \
   \
    int nargin = args.length (); \
   \
!   if (nargin < 1 || nargin > 2 || nargout > 2) \
      { \
        print_usage (#FCN); \
        return retval; \
--- 37,51 ----
  #include "gripes.h"
  #include "oct-obj.h"
  
+ #include "ov-cx-mat.h"
+ 
  #define MINMAX_BODY(FCN) \
   \
    octave_value_list retval;  \
   \
    int nargin = args.length (); \
   \
!   if (nargin < 1 || nargin > 3 || nargout > 2) \
      { \
        print_usage (#FCN); \
        return retval; \
***************
*** 51,59 ****
--- 53,65 ----
   \
    octave_value arg1; \
    octave_value arg2; \
+   octave_value arg3; \
   \
    switch (nargin) \
      { \
+     case 3: \
+       arg3 = args(2); \
+  \
      case 2: \
        arg2 = args(1); \
   \
***************
*** 66,162 ****
        break; \
      } \
   \
!   if (nargin == 1 && (nargout == 1 || nargout == 0)) \
      { \
        if (arg1.is_real_type ()) \
  	{ \
! 	  Matrix m = arg1.matrix_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      if (m.rows () == 1) \
! 		retval(0) = m.row_ ## FCN (); \
! 	      else \
! 		{ \
! 		  if (m.rows () == 0 || m.columns () == 0) \
! 		    retval(0) = Matrix (); \
! 		  else \
! 		    retval(0) = m.column_ ## FCN (); \
! 		} \
  	    } \
  	} \
        else if (arg1.is_complex_type ()) \
  	{ \
! 	  ComplexMatrix m = arg1.complex_matrix_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      if (m.rows () == 1) \
! 		retval(0) = m.row_ ## FCN (); \
! 	      else \
! 		{ \
! 		  if (m.rows () == 0 || m.columns () == 0) \
! 		    retval(0) = Matrix (); \
! 		  else \
! 		    retval(0) = m.column_ ## FCN (); \
! 		} \
  	    } \
  	} \
        else \
  	gripe_wrong_type_arg (#FCN, arg1); \
      } \
!   else if (nargin == 1 && nargout == 2) \
      { \
!       Array<int> index; \
   \
        if (arg1.is_real_type ()) \
  	{ \
! 	  Matrix m = arg1.matrix_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      retval.resize (2); \
!  \
! 	      if (m.rows () == 1) \
! 		retval(0) = m.row_ ## FCN (index); \
! 	      else \
! 		{ \
! 		  if (m.rows () == 0 || m.columns () == 0) \
! 		    retval(0) = Matrix (); \
! 		  else \
! 		    retval(0) = m.column_ ## FCN (index); \
! 		} \
  	    } \
  	} \
        else if (arg1.is_complex_type ()) \
  	{ \
! 	  ComplexMatrix m = arg1.complex_matrix_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      retval.resize (2); \
!  \
! 	      if (m.rows () == 1) \
! 		retval(0) = m.row_ ## FCN (index); \
! 	      else \
! 		{ \
! 		  if (m.rows () == 0 || m.columns () == 0) \
! 		    retval(0) = Matrix (); \
! 		  else \
! 		    retval(0) = m.column_ ## FCN (index); \
! 		} \
  	    } \
  	} \
        else \
  	gripe_wrong_type_arg (#FCN, arg1); \
   \
!       int len = index.length (); \
   \
        if (len > 0) \
  	{ \
  	  double nan_val = lo_ieee_nan_value (); \
   \
! 	  RowVector idx (len); \
   \
  	  for (int i = 0; i < len; i++) \
  	    { \
--- 72,182 ----
        break; \
      } \
   \
!   int dim; \
!   dim_vector dv = ((const octave_complex_matrix&) arg1) .dims (); \
!   if (error_state) \
!     { \
!       gripe_wrong_type_arg (#FCN, arg1);  \
!       return retval; \
!     } \
!  \
!   if (nargin == 3) \
!     { \
!       dim = arg3.nint_value () - 1;  \
!       if (dim < 0 || dim >= dv.length ()) \
!         { \
! 	  error ("%s: invalid dimension", #FCN); \
! 	  return retval; \
! 	} \
!     } \
!   else \
!     { \
!       dim = 0; \
!       while ((dim < dv.length ()) && (dv (dim) <= 1)) \
! 	dim++; \
!       if (dim == dv.length ()) \
! 	dim = 0; \
!     } \
!  \
!   bool single_arg = (nargin == 1) || arg2.is_empty();	\
!  \
!   if (single_arg) \
!     { \
!       dv(dim) = 1; \
!       int n_dims = dv.length (); \
!       for (int i = n_dims; i > 1; i--) \
! 	{ \
! 	  if (dv(i-1) == 1) \
! 	    n_dims--; \
! 	  else \
! 	    break; \
! 	} \
!       dv.resize (n_dims); \
!     } \
!  \
!   if (single_arg && (nargout == 1 || nargout == 0)) \
      { \
        if (arg1.is_real_type ()) \
  	{ \
! 	  NDArray m = arg1.array_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      NDArray n = m. FCN (dim); \
! 	      n.resize (dv); \
! 	      retval(0) = n; \
  	    } \
  	} \
        else if (arg1.is_complex_type ()) \
  	{ \
! 	  ComplexNDArray m = arg1.complex_array_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      ComplexNDArray n = m. FCN (dim); \
! 	      n.resize (dv); \
! 	      retval(0) = n; \
  	    } \
  	} \
        else \
  	gripe_wrong_type_arg (#FCN, arg1); \
      } \
!   else if (single_arg && nargout == 2) \
      { \
!       ArrayN<int> index; \
   \
        if (arg1.is_real_type ()) \
  	{ \
! 	  NDArray m = arg1.array_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      NDArray n = m. FCN (index, dim);	\
! 	      n.resize (dv); \
! 	      retval(0) = n; \
  	    } \
  	} \
        else if (arg1.is_complex_type ()) \
  	{ \
! 	  ComplexNDArray m = arg1.complex_array_value (); \
   \
  	  if (! error_state) \
  	    { \
! 	      ComplexNDArray n = m. FCN (index, dim);	\
! 	      n.resize (dv); \
! 	      retval(0) = n; \
  	    } \
  	} \
        else \
  	gripe_wrong_type_arg (#FCN, arg1); \
   \
!       int len = index.numel (); \
   \
        if (len > 0) \
  	{ \
  	  double nan_val = lo_ieee_nan_value (); \
   \
! 	  NDArray idx (index.dims ()); \
   \
  	  for (int i = 0; i < len; i++) \
  	    { \
***************
*** 169,177 ****
  	  retval(1) = idx; \
  	} \
        else \
! 	retval(1) = Matrix (); \
      } \
!   else if (nargin == 2) \
      { \
        int arg1_is_scalar = arg1.is_scalar_type (); \
        int arg2_is_scalar = arg2.is_scalar_type (); \
--- 189,197 ----
  	  retval(1) = idx; \
  	} \
        else \
! 	retval(1) = NDArray (); \
      } \
!   else \
      { \
        int arg1_is_scalar = arg1.is_scalar_type (); \
        int arg2_is_scalar = arg2.is_scalar_type (); \
***************
*** 184,193 ****
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
  	      Complex c1 = arg1.complex_value (); \
! 	      ComplexMatrix m2 = arg2.complex_matrix_value (); \
  	      if (! error_state) \
  		{ \
! 		  ComplexMatrix result = FCN (c1, m2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
--- 204,213 ----
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
  	      Complex c1 = arg1.complex_value (); \
! 	      ComplexNDArray m2 = arg2.complex_array_value (); \
  	      if (! error_state) \
  		{ \
! 		  ComplexNDArray result = FCN (c1, m2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
***************
*** 195,205 ****
  	  else \
  	    { \
  	      double d1 = arg1.double_value (); \
! 	      Matrix m2 = arg2.matrix_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  Matrix result = FCN (d1, m2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
--- 215,225 ----
  	  else \
  	    { \
  	      double d1 = arg1.double_value (); \
! 	      NDArray m2 = arg2.array_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  NDArray result = FCN (d1, m2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
***************
*** 209,232 ****
  	{ \
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
! 	      ComplexMatrix m1 = arg1.complex_matrix_value (); \
   \
  	      if (! error_state) \
  		{ \
  		  Complex c2 = arg2.complex_value (); \
! 		  ComplexMatrix result = FCN (m1, c2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
  	    } \
  	  else \
  	    { \
! 	      Matrix m1 = arg1.matrix_value (); \
   \
  	      if (! error_state) \
  		{ \
  		  double d2 = arg2.double_value (); \
! 		  Matrix result = FCN (m1, d2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
--- 229,252 ----
  	{ \
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
! 	      ComplexNDArray m1 = arg1.complex_array_value (); \
   \
  	      if (! error_state) \
  		{ \
  		  Complex c2 = arg2.complex_value (); \
! 		  ComplexNDArray result = FCN (m1, c2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
  	    } \
  	  else \
  	    { \
! 	      NDArray m1 = arg1.array_value (); \
   \
  	      if (! error_state) \
  		{ \
  		  double d2 = arg2.double_value (); \
! 		  NDArray result = FCN (m1, d2); \
  		  if (! error_state) \
  		    retval(0) = result; \
  		} \
***************
*** 236,250 ****
  	{ \
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
! 	      ComplexMatrix m1 = arg1.complex_matrix_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  ComplexMatrix m2 = arg2.complex_matrix_value (); \
   \
  		  if (! error_state) \
  		    { \
! 		      ComplexMatrix result = FCN (m1, m2); \
  		      if (! error_state) \
  			retval(0) = result; \
  		    } \
--- 256,270 ----
  	{ \
  	  if (arg1_is_complex || arg2_is_complex) \
  	    { \
! 	      ComplexNDArray m1 = arg1.complex_array_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  ComplexNDArray m2 = arg2.complex_array_value (); \
   \
  		  if (! error_state) \
  		    { \
! 		      ComplexNDArray result = FCN (m1, m2); \
  		      if (! error_state) \
  			retval(0) = result; \
  		    } \
***************
*** 252,266 ****
  	    } \
  	  else \
  	    { \
! 	      Matrix m1 = arg1.matrix_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  Matrix m2 = arg2.matrix_value (); \
   \
  		  if (! error_state) \
  		    { \
! 		      Matrix result = FCN (m1, m2); \
  		      if (! error_state) \
  			retval(0) = result; \
  		    } \
--- 272,286 ----
  	    } \
  	  else \
  	    { \
! 	      NDArray m1 = arg1.array_value (); \
   \
  	      if (! error_state) \
  		{ \
! 		  NDArray m2 = arg2.array_value (); \
   \
  		  if (! error_state) \
  		    { \
! 		      NDArray result = FCN (m1, m2); \
  		      if (! error_state) \
  			retval(0) = result; \
  		    } \
***************
*** 268,288 ****
  	    } \
  	} \
      } \
-   else \
-     panic_impossible (); \
   \
    return retval
  
  DEFUN_DLD (min, args, nargout,
    "-*- texinfo -*-\n\
! @deftypefn {Mapping Function} {} min (@var{x}, @var{y})\n\
  @deftypefnx {Mapping Function} address@hidden, @var{iw}] =} min (@var{x})\n\
  @cindex Utility Functions\n\
  For a vector argument, return the minimum value.  For a matrix\n\
  argument, return the minimum value from each column, as a row\n\
! vector.\n\
! For two matrices (or a matrix and scalar),\n\
! return the pair-wise minimum.\n\
  Thus,\n\
  \n\
  @example\n\
--- 288,305 ----
  	    } \
  	} \
      } \
   \
    return retval
  
  DEFUN_DLD (min, args, nargout,
    "-*- texinfo -*-\n\
! @deftypefn {Mapping Function} {} min (@var{x}, @var{y}, @var{dim})\n\
  @deftypefnx {Mapping Function} address@hidden, @var{iw}] =} min (@var{x})\n\
  @cindex Utility Functions\n\
  For a vector argument, return the minimum value.  For a matrix\n\
  argument, return the minimum value from each column, as a row\n\
! vector, or over the dimension @var{dim} if defined. For two matrices\n\
! (or a matrix and scalar), return the pair-wise minimum.\n\
  Thus,\n\
  \n\
  @example\n\
***************
*** 323,336 ****
  
  DEFUN_DLD (max, args, nargout,
    "-*- texinfo -*-\n\
! @deftypefn {Mapping Function} {} max (@var{x}, @var{y})\n\
  @deftypefnx {Mapping Function} address@hidden, @var{iw}] =} max (@var{x})\n\
  @cindex Utility Functions\n\
  For a vector argument, return the maximum value.  For a matrix\n\
  argument, return the maximum value from each column, as a row\n\
! vector.\n\
! For two matrices (or a matrix and scalar),\n\
! return the pair-wise maximum.\n\
  Thus,\n\
  \n\
  @example\n\
--- 340,352 ----
  
  DEFUN_DLD (max, args, nargout,
    "-*- texinfo -*-\n\
! @deftypefn {Mapping Function} {} max (@var{x}, @var{y}, @var{dim})\n\
  @deftypefnx {Mapping Function} address@hidden, @var{iw}] =} max (@var{x})\n\
  @cindex Utility Functions\n\
  For a vector argument, return the maximum value.  For a matrix\n\
  argument, return the maximum value from each column, as a row\n\
! vector, or over the dimension @var{dim} if defined. For two matrices\n\
! (or a matrix and scalar), return the pair-wise maximum.\n\
  Thus,\n\
  \n\
  @example\n\
*** ./src/DLD-FUNCTIONS/besselj.cc.orig	2004-03-29 15:08:42.000000000 +0200
--- ./src/DLD-FUNCTIONS/besselj.cc	2004-03-29 15:21:21.000000000 +0200
***************
*** 97,102 ****
--- 97,120 ----
    return retval;
  }
  
+ static inline NDArray
+ int_arrayN_to_array (const ArrayN<int>& a)
+ {
+   dim_vector dv = a.dims ();
+   int nel = dv.numel ();
+ 
+   NDArray retval (dv);
+ 
+   for (int i = 0; i < nel; i++)
+     {
+       OCTAVE_QUIT;
+       
+       retval(i) = (double) (a(i));
+     }
+ 
+   return retval;
+ }
+ 
  static void
  gripe_bessel_arg (const char *fn, const char *arg)
  {
***************
*** 145,161 ****
  		}
  	      else
  		{
! 		  ComplexMatrix x = x_arg.complex_matrix_value ();
  
  		  if (! error_state)
  		    {
! 		      Array2<int> ierr;
  		      octave_value result;
  
  		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
  
  		      if (nargout > 1)
! 			retval(1) = int_array2_to_matrix (ierr);
  
  		      retval(0) = result;
  		    }
--- 163,179 ----
  		}
  	      else
  		{
! 		  ComplexNDArray x = x_arg.complex_array_value ();
  
  		  if (! error_state)
  		    {
! 		      ArrayN<int> ierr;
  		      octave_value result;
  
  		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
  
  		      if (nargout > 1)
! 			retval(1) = int_arrayN_to_array (ierr);
  
  		      retval(0) = result;
  		    }
***************
*** 168,188 ****
  	}
        else
  	{
! 	  Matrix alpha = args(0).matrix_value ();
  
! 	  if (! error_state)
  	    {
! 	      if (x_arg.is_scalar_type ())
  		{
! 		  Complex x = x_arg.complex_value ();
  
  		  if (! error_state)
  		    {
  		      Array2<int> ierr;
  		      octave_value result;
  
! 		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
! 
  		      if (nargout > 1)
  			retval(1) = int_array2_to_matrix (ierr);
  
--- 186,213 ----
  	}
        else
  	{
! 	  dim_vector dv0 = args(0).dims ();
! 	  dim_vector dv1 = args(1).dims ();
! 	  
! 	  bool args0_is_row_vector = (dv0 (1) == dv0.numel ());
! 	  bool args1_is_col_vector = (dv1 (0) == dv1.numel ());
  
! 	  if (args0_is_row_vector && args1_is_col_vector)
  	    {
! 	      RowVector ralpha = args(0).row_vector_value ();
! 
! 	      if (! error_state)
  		{
! 		  ComplexColumnVector cx = 
! 		    x_arg.complex_column_vector_value ();
  
  		  if (! error_state)
  		    {
  		      Array2<int> ierr;
  		      octave_value result;
  
! 		      DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
! 		      
  		      if (nargout > 1)
  			retval(1) = int_array2_to_matrix (ierr);
  
***************
*** 192,236 ****
  		    gripe_bessel_arg (fn, "second");
  		}
  	      else
! 		{
! 		  ComplexMatrix x = x_arg.complex_matrix_value ();
  
! 		  if (! error_state)
  		    {
! 		      if (alpha.rows () == 1 && x.columns () == 1)
! 			{
! 			  RowVector ralpha = alpha.row (0);
! 			  ComplexColumnVector cx = x.column (0);
  
! 			  Array2<int> ierr;
  			  octave_value result;
  
! 			  DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
  
  			  if (nargout > 1)
! 			    retval(1) = int_array2_to_matrix (ierr);
  
  			  retval(0) = result;
  			}
  		      else
  			{
! 			  Array2<int> ierr;
  			  octave_value result;
! 
  			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
! 
  			  if (nargout > 1)
! 			    retval(1) = int_array2_to_matrix (ierr);
! 
  			  retval(0) = result;
  			}
  		    }
- 		  else
- 		    gripe_bessel_arg (fn, "second");
  		}
  	    }
- 	  else
- 	    gripe_bessel_arg (fn, "first");
  	}
      }
    else
--- 217,272 ----
  		    gripe_bessel_arg (fn, "second");
  		}
  	      else
! 		gripe_bessel_arg (fn, "first");
! 	    }
! 	  else
! 	    {
! 	      NDArray alpha = args(0).array_value ();
  
! 	      if (! error_state)
! 		{
! 		  if (x_arg.is_scalar_type ())
  		    {
! 		      Complex x = x_arg.complex_value ();
  
! 		      if (! error_state)
! 			{
! 			  ArrayN<int> ierr;
  			  octave_value result;
  
! 			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
  
  			  if (nargout > 1)
! 			    retval(1) = int_arrayN_to_array (ierr);
  
  			  retval(0) = result;
  			}
  		      else
+ 			gripe_bessel_arg (fn, "second");
+ 		    }
+ 		  else
+ 		    {
+ 		      ComplexNDArray x = x_arg.complex_array_value ();
+ 
+ 		      if (! error_state)
  			{
! 			  ArrayN<int> ierr;
  			  octave_value result;
! 			  
  			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
! 			  
  			  if (nargout > 1)
! 			    retval(1) = int_arrayN_to_array (ierr);
! 			  
  			  retval(0) = result;
  			}
+ 		      else
+ 			gripe_bessel_arg (fn, "second");
  		    }
  		}
+ 	      else
+ 		gripe_bessel_arg (fn, "first");
  	    }
  	}
      }
    else
***************
*** 422,428 ****
  
        int kind = 0;
  
!       ComplexMatrix z;
  
        if (nargin > 1)
  	{
--- 458,464 ----
  
        int kind = 0;
  
!       ComplexNDArray z;
  
        if (nargin > 1)
  	{
***************
*** 441,451 ****
  
        if (! error_state)
  	{
! 	  z = args(nargin == 1 ? 0 : 1).complex_matrix_value ();
  
  	  if (! error_state)
  	    {
! 	      Array2<int> ierr;
  	      octave_value result;
  
  	      if (kind > 1)
--- 477,487 ----
  
        if (! error_state)
  	{
! 	  z = args(nargin == 1 ? 0 : 1).complex_array_value ();
  
  	  if (! error_state)
  	    {
! 	      ArrayN<int> ierr;
  	      octave_value result;
  
  	      if (kind > 1)
***************
*** 454,460 ****
  		result = airy (z, kind == 1, scale, ierr);
  
  	      if (nargout > 1)
! 		retval(1) = int_array2_to_matrix (ierr);
  
  	      retval(0) = result;
  	    }
--- 490,496 ----
  		result = airy (z, kind == 1, scale, ierr);
  
  	      if (nargout > 1)
! 		retval(1) = int_arrayN_to_array (ierr);
  
  	      retval(0) = result;
  	    }
*** ./src/DLD-FUNCTIONS/gammainc.cc.orig	2004-03-29 15:25:53.000000000 +0200
--- ./src/DLD-FUNCTIONS/gammainc.cc	2004-03-29 15:25:55.000000000 +0200
***************
*** 86,92 ****
  		}
  	      else
  		{
! 		  Matrix a = a_arg.matrix_value ();
  
  		  if (! error_state)
  		    retval = gammainc (x, a);
--- 86,92 ----
  		}
  	      else
  		{
! 		  NDArray a = a_arg.array_value ();
  
  		  if (! error_state)
  		    retval = gammainc (x, a);
***************
*** 95,101 ****
  	}
        else
  	{
! 	  Matrix x = x_arg.matrix_value ();
  
  	  if (! error_state)
  	    {
--- 95,101 ----
  	}
        else
  	{
! 	  NDArray x = x_arg.array_value ();
  
  	  if (! error_state)
  	    {
***************
*** 108,114 ****
  		}
  	      else
  		{
! 		  Matrix a = a_arg.matrix_value ();
  
  		  if (! error_state)
  		    retval = gammainc (x, a);
--- 108,114 ----
  		}
  	      else
  		{
! 		  NDArray a = a_arg.array_value ();
  
  		  if (! error_state)
  		    retval = gammainc (x, a);
*** ./src/ChangeLog.orig	2004-03-26 15:55:00.000000000 +0100
--- ./src/ChangeLog	2004-03-29 15:35:26.000000000 +0200
***************
*** 1,3 ****
--- 1,22 ----
+ 2004-03-29 David Bateman <address@hidden>
+ 
+ 	* DLD-FUNCTIONS/besselj.cc: Convert for NDArray, better Matlab
+ 	compatibility.
+ 	* DLD-FUNCTIONS/betainc.cc: ditto
+ 	* DLD-FUNCTIONS/gammainc.cc: ditto
+ 	
+ 2004-03-26 David Bateman <address@hidden>
+ 
+ 	* DLD-FUNCTIONS/minmax.cc: Convert for NDArray, better Matlab
+ 	compatibility.
+ 
+ 	* DLD-FUNCTIONS/filter.cc: Convert for NDArray, better Matlab
+ 	compatibility.
+ 	
+ 2004-03-25 David Bateman <address@hidden>
+ 
+ 	* load-save.cc (Fload): Better handling of non existent files.
+ 	
  2004-03-12  John W. Eaton  <address@hidden>
  
  	* ov-cell.cc (octave_cell::save_hdf5): Handle empty cells.
*** ./src/load-save.cc.orig	2004-02-20 22:16:54.000000000 +0100
--- ./src/load-save.cc	2004-03-25 13:10:12.000000000 +0100
***************
*** 306,311 ****
--- 306,318 ----
  {
    load_save_format retval = LS_UNKNOWN;
  
+   // If the file doesn't exist do nothing
+   std::ifstream file_exist (fname.c_str ());
+   if (file_exist)
+     file_exist.close ();
+   else
+     return LS_UNKNOWN;
+ 
  #ifdef HAVE_HDF5
    // check this before we open the file
    if (H5Fis_hdf5 (fname.c_str ()) > 0)
***************
*** 722,740 ****
  	{
  	  i++;
  
! 	  hdf5_ifstream hdf5_file (fname.c_str ());
! 
! 	  if (hdf5_file.file_id >= 0)
  	    {
! 	      retval = do_load (hdf5_file, orig_fname, force, format,
! 				flt_fmt, list_only, swap, verbose,
! 				argv, i, argc, nargout);
  
! 	      hdf5_file.close ();
  	    }
- 	  else
- 	    error ("load: couldn't open input file `%s'",
- 		   orig_fname.c_str ());
  	}
        else
  #endif /* HAVE_HDF5 */
--- 729,754 ----
  	{
  	  i++;
  
! 	  // If the file doesn't exist do nothing
! 	  std::ifstream file (fname.c_str (), std::ios::in);
! 	  if (file)
  	    {
! 	      file.close ();
! 	      
! 	      hdf5_ifstream hdf5_file (fname.c_str ());
  
! 	      if (hdf5_file.file_id >= 0)
! 		{
! 		  retval = do_load (hdf5_file, orig_fname, force, format,
! 				    flt_fmt, list_only, swap, verbose,
! 				    argv, i, argc, nargout);
! 
! 		  hdf5_file.close ();
! 		}
! 	      else
! 		error ("load: couldn't open input file `%s'",
! 		       orig_fname.c_str ());
  	    }
  	}
        else
  #endif /* HAVE_HDF5 */