#ifndef __VEC_H
#define __VEC_H
/*!
  \file    vec.hh
  \brief   static vector
  \date    30 April 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/

#include <cassert>
#include <cmath>
#include <iostream>

using namespace std;

template <int Size, typename T=double>
class Vec
{
public:
  /// default constructor: set vector to zero
  Vec()
  {
    for (int i=0; i<Size; ++i) data_[i] = 0;
  }

  /// alternative constructor, works only for Vec<2> and Vec<3>
  Vec(T v0, T v1, T v2=0)
  {
    data_[0] = v0;
    data_[1] = v1;
    if (Size == 3) data_[2] = v2;
  }
  
  /// destructor
  ~Vec()
  {
    ; 
  }

  /// copy constructor
  Vec (const Vec& v)
  {
    for (int i=0; i<Size; ++i) data_[i] = v.data_[i];
  }
  
  /// assignment
  Vec& operator= (const Vec& v)
  {
    for (int i=0; i<Size; ++i) data_[i] = v.data_[i];
    return *this;
  }

  /// assign all entries to val
  Vec& operator= (const T val)
  {
    for (int i=0; i<Size; ++i) data_[i] = val;
    return *this;
  }
  
  /// element access
  const T operator[] (int i) const
  {
    assert (0 <= i && i < Size && "range check");
    return data_[i]; 
  }
  
  /// element access by reference
  T& operator[] (int i)
  {
    assert (0 <= i && i < Size && "range check");
    return data_[i]; 
  }
  
  /// set first entries, works only for Size >= 2
  void set (T v0, T v1, T v2=0)
  {
    data_[0] = v0;
    data_[1] = v1;
    if (Size > 2) data_[2] = v2;
  }

  /// increment
  Vec& operator+= (const Vec& v)
  {
    for (int i=0; i<Size; ++i) data_[i] += v.data_[i];
    return *this;
  }
  
  /// decrement
  Vec& operator-= (const Vec& v)
  {
    for (int i=0; i<Size; ++i) data_[i] -= v.data_[i];
    return *this;
  }
  
  /// scale
  Vec& operator*= (const T factor)
  {
    for (int i=0; i<Size; ++i) data_[i] *= factor;
    return *this;
  }

  /// output
  void print (std::ostream& os) const
  {
    os << "[";
    for (int i=0; i<Size; ++i) os << ' ' << data_[i];
    os << " ]";
  }
  
private:
  
  T data_[Size];
  
}; // class Vec


  
/// vector addition
template <int Size>
Vec<Size> operator+ (const Vec<Size>& v, const Vec<Size>& w)
{
  Vec<Size> res(v); // copy
  res += w;
  return res;
}

 
/// vector subtraction
template <int Size>
Vec<Size> operator- (const Vec<Size>& v, const Vec<Size>& w)
{
  Vec<Size> res(v); // copy
  res -= w;
  return res;
}


/// scale vector
template <int Size, typename T>
Vec<Size, T> operator* (const T factor, const Vec<Size>& v)
{
  Vec<Size, T> res(v);
  res *= factor;
  return res;
}


/// Euclidean inner product
template <int Size, typename T>
T inner_product (const Vec<Size, T>& v, const Vec<Size, T>& w)
{
  T sum = 0;
  for (int i=0; i<Size; ++i) sum += v[i] * w[i];
  return sum;
}

/// cross product
Vec<3> cross_product (const Vec<3>& a, const Vec<3>& b)
{
  Vec<3> c;
  c[0] = a[1] * b[2] - a[2] * b[1];
  c[1] = a[2] * b[0] - a[0] * b[2];
  c[2] = a[0] * b[1] - a[1] * b[0];
  return c;
}

/// Euclidean norm
template <int Size>
double l2norm (const Vec<Size, double>& v)
{
  double norm=0;
  for (int i=0; i<Size; ++i) norm += v[i]*v[i];
  return sqrt(norm);
}


/// output
template <int Size, typename T>
std::ostream& operator<< (std::ostream& os, const Vec<Size, T>& v)
{
  v.print (os);
  return os;
}


typedef Vec<2, double> Point2D;


#endif // __VEC_H