//  ComPhys file: Vector.h

#ifndef _VECTOR_H
#define _VECTOR_H 1

#include <iostream.h>
#include <stdlib.h>

template<class T> 
class Vector {
public:
    // constructors
    
    Vector() {
        _nullify();
    }

    Vector(int dimension) {
        _nullify();
        _resize(dimension);
    }

    Vector(const Vector& vector) {
        _dimension = vector._dimension;
        _components = _allocate(_dimension);
        for (int dim = 0; dim < _dimension; ++dim)
            _components[dim] = vector._components[dim];
    }

    // set or change dimension

    void setDimension(int dimension) {
        _resize(dimension);
    }

    // get dimension

    int dimension() const {
        return _dimension;
    }

    // destructor

    ~Vector() {
        delete [] _components;
    }

    // allow access to components by dimension

    T& operator()(int component) {
        if (component < 1 || component > _dimension) {
            cerr << "Vector: non-existent component " << component << endl;
            exit(1);
        }
        return _components[component-1];
    }

    // allow access to components by index

    T& operator[](int index) {
        if (index < 0 || index >= _dimension) {
            cerr << "Vector: bad index " << index << endl;
            exit(1);
        }
        return _components[index];
    }

    // assignment operators

    Vector& operator=(const Vector& vector) {
        _checkDimension(vector._dimension);
        for (int dim = 0; dim < _dimension; ++dim)
            _components[dim] = vector._components[dim];
        return *this;
    }

    Vector& operator+=(const Vector& vector) {
        _checkDimension(vector._dimension);
        for (int dim = 0; dim < _dimension; ++dim)
            _components[dim] += vector._components[dim];
        return *this;
    }

    Vector& operator-=(const Vector& vector) {
        _checkDimension(vector._dimension);
        for (int dim = 0; dim < _dimension; ++dim)
            _components[dim] -= vector._components[dim];
        return *this;
    }

private:
    int _dimension;
    T* _components;

    void _nullify() {
        _dimension = 0;
        _components = NULL;
    }

    T* _allocate(int dimension) {
        T* components = new T[dimension];
        if (components == NULL) {
            cerr << "Vector: memory allocation failure" << endl;
            exit(1);
        }
        return components;
    }
    
    void _resize(int dimension) {
        if (dimension == _dimension)
            return;
        if (dimension < 1) {
            cerr << "Vector: dimension " << dimension << " < 1" << endl;
            exit(1);
        }
        delete [] _components;
        _dimension = dimension;
        _components = _allocate(dimension);
    }

    void _checkDimension(int dimension) {
        if (dimension != _dimension) {
            cerr << "Vector: assignment: dimension " << _dimension
                 << " not = " << dimension << endl;
            exit(1);
        }
    }
    
    // friends
    
	friend istream& operator>>(istream& is, Vector& vector) {
        for (int dim = 0; dim < vector._dimension; ++dim)
            is >> vector._components[dim];
        return is;
    }
    
	friend ostream& operator<<(ostream& os, Vector& vector) {
        for (int dim = 0; dim < vector._dimension-1; ++dim)
            os << vector._components[dim] << '\t';
        os << vector._components[vector._dimension-1];
        return os;
    }
    
};

//	overload operators for vector algebra

template<class T>
Vector<T> operator+(const Vector<T>& vector) {
    Vector<T> v = vector;
    return v;
}

template<class T>
Vector<T> operator-(const Vector<T>& vector) {
    Vector<T> v = vector;
    for (int dim = 0; dim < v.dimension(); ++dim)
        v[dim] = -v[dim];
    return v;
}

template<class T, class S>
Vector<T>& operator*=(Vector<T>& vector, const S scalar) {
    for (int dim = 0; dim < vector.dimension(); ++dim) 
        vector[dim] *= scalar;
    return vector;
}

template<class T, class S>
Vector<T>& operator/=(Vector<T>& vector, const S scalar) {
    for (int dim = 0; dim < vector.dimension(); ++dim) 
        vector[dim] /= scalar;
    return vector;
}

template<class T, class S>
Vector<T> operator*(S scalar, Vector<T> vector) {
    int dimension = vector.dimension();
    for (int dim = 0; dim < dimension; ++dim)
        vector[dim] *= scalar;
    return vector;
}

template<class T, class S>
Vector<T> operator/(Vector<T> vector, S scalar) {
    int dimension = vector.dimension();
    for (int dim = 0; dim < dimension; ++dim)
        vector[dim] /= scalar;
    return vector;
}

template<class T, class S>
Vector<T> operator*(Vector<T> vector, S scalar) {
    int dimension = vector.dimension();
    for (int dim = 0; dim < dimension; ++dim)
        vector[dim] *= scalar;
    return vector;
}

inline void _compareDimensions(int dimension1, int dimension2) {
    if (dimension1 != dimension2) {
        cerr << "Vector algebra: dimension " << dimension1
             << " not = " << dimension2 << endl;
        exit(1);
    }
}

template<class T>
Vector<T> operator+(const Vector<T>& vector1, const Vector<T>& vector2) {
    int dimension = vector1.dimension();
    _compareDimensions(dimension, vector2.dimension());
    Vector<T> vector = vector1;
    vector += vector2;
    return vector;
}

template<class T>
Vector<T> operator-(const Vector<T>& vector1, const Vector<T>& vector2) {
    int dimension = vector1.dimension();
    _compareDimensions(dimension, vector2.dimension());
    Vector<T> vector = vector1;
    vector -= vector2;
    return vector;
}

#endif _VECTOR_H