30d29
< #include <complex>
422,727d420
<     // **********************************************************************
< 		/* complex vector classes to follow */
< 
<     template <int N>
<     CVectorBase<N>::CVectorBase(std::complex<double> value)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] = value;
<     }
< 
<     template <int N>
<     void CVectorBase<N>::set(std::complex<double> value)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] = value;
<     }
< 
<     template <int N>
<     std::complex<double> & CVectorBase<N>::operator[](int n)
<     {
<       return _val[n];
<     }
< 
<     template <int N>
<     std::complex<double> CVectorBase<N>::operator[](int n) const
<     {
<       return _val[n];
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator+(const CVectorBase<N> &v) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = _val[i] + v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator-(const CVectorBase<N> &v) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = _val[i] - v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator-() const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = -_val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> & CVectorBase<N>::neg()
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] = -_val[i];
< 
<       return *this;
<     }
< 
<     template <int N>
<     const CVectorBase<N> & CVectorBase<N>::operator+=(const CVectorBase<N> &v)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] += v._val[i];
< 
<       return *this;
<     }
< 
<     template <int N>
<     const CVectorBase<N> & CVectorBase<N>::operator-=(const CVectorBase<N> &v)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] -= v._val[i];
< 
<       return *this;
<     }
< 
<     template <int N>
<     std::complex<double> CVectorBase<N>::operator*(const CVectorBase<N> &v) const
<     {
<       std::complex<double> r = 0;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r += _val[i] * v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::mul(const CVectorBase &v) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = _val[i] * v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator/(const CVectorBase &v) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = _val[i] / v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator*(std::complex<double> scale) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = scale * _val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator/(std::complex<double> scale) const
<     {
<       Math::CVectorBase<N> r;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = _val[i] / scale;
< 
<       return r;
<     }
< 
<     template <int N>
<     const CVectorBase<N> & CVectorBase<N>::operator*=(std::complex<double> scale)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] *= scale;
< 
<       return *this;
<     }
< 
<     template <int N>
<     const CVectorBase<N> & CVectorBase<N>::operator/=(std::complex<double> scale)
<     {
<       for (unsigned int i = 0; i < N; i++)
< 	_val[i] /= scale;
< 
<       return *this;
<     }
< 
<     template <int N>
<     std::complex<double> CVectorBase<N>::len() const
<     {
<       std::complex<double> r = 0;
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r += square(_val[i]);
< 
<       return sqrt(r);
<     }
< 
<     template <int N>
<     const CVectorBase<N> & CVectorBase<N>::normalize()
<     {
<       return *this /= len();
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::normalized() const
<     {
<       CVectorBase<N> r = *this / len();
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::magnitude(std::complex<double> newlen) const
<     {
<       return *this * (newlen / len());
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::select(unsigned int index, const CVectorBase<N> v) const
<     {
<       Math::CVectorBase<N> r;
< 
<       assert(index < N);
< 
<       for (unsigned int i = 0; i < N; i++)
< 	r._val[i] = (i == index) ? _val[i] : v._val[i];
< 
<       return r;
<     }
< 
<     template <int N>
<     CVectorBase<N> CVectorBase<N>::operator*(const Matrix<N> &m)
<     {
<       CVectorBase<N> r;
< 
<       for (int j = 0; j < N; j++)
< 	{
< 	  std::complex<double> s = 0;
< 
< 	  for (int k = 0; k < N; k++)
< 	    s += _val[k][j] * m._val[k];
< 
< 	  r[j] = s;
< 	}
< 
<       return r;
<     }
< 
<     template <int N>
<     inline bool CVectorBase<N>::operator==(const CVectorBase &m) const
<     {
<       for (int j = 0; j < N; j++)
< 	if (_val[j] != m._val[j])
< 	  return false;
< 
<       return true;
<     }
< 
<     template <int N>
<     std::ostream & operator<<(std::ostream &o, const CVectorBase<N> &v)
<     {
<       o << "[";
< 
<       for (unsigned int i = 0; i < N; i++)
< 	{
< 	  o << v[i];
< 	  if (i + 1 < N)
< 	    o << ", ";
< 	}
< 
<       o << "]";
< 
<       return o;
<     }
< 
<     // **********************************************************************
< 
<     template <int N>
<     CVector<N>::CVector()
<     {
<     }
< 
<     template <int N>
<     CVector<N>::CVector(const CVectorBase<N> &v)
<       : CVectorBase<N>(v)
<     {
<     }
< 
< 
<     // **********************************************************************
< 
<     CVector<2>::CVector()
<     {
<     }
< 
<     CVector<2>::CVector(const CVectorBase<2> &v)
<       : CVectorBase<2>(v)
<     {
<     }
< 
<     CVector<2>::CVector(std::complex<double> x, std::complex<double> y)
<     {
<       _val[0] = x;
<       _val[1] = y;
<     }
< 
<     std::complex<double> CVector<2>::perp_product(const CVector<2> &v) const
<     {
<       return _val[0] * v._val[1] - _val[1] * v._val[0];
<     }
< 
<     std::complex<double> & CVector<2>::x()
<     {
<       return _val[0];
<     }
< 
<     std::complex<double> & CVector<2>::y()
<     {
<       return _val[1];
<     }
< 
<     std::complex<double> CVector<2>::x() const
<     {
<       return _val[0];
<     }
< 
<     std::complex<double> CVector<2>::y() const
<     {
<       return _val[1];
<     }
<