#ifndef VECTORFIELD3D_H
#define VECTORFIELD3D_H

#include <iostream.h>
#include <assert.h>

class VectorField3D;
class Mask3D;


class VectorField3D
{
    public:
	VectorField3D(int,int,int);
	VectorField3D();
	~VectorField3D();
	void reset(int,int,int);

	int nx() const;
	int ny() const;
	int nz() const;
	int numX() const;
	int numY() const;
	int numZ() const;

	float valX(int,int,int) const;
	float valY(int,int,int) const;
	float valZ(int,int,int) const;
	float& valXRef(int,int,int) const;
	float& valYRef(int,int,int) const;
	float& valZRef(int,int,int) const;

	float uncenteredValX(int,int,int) const;
	float uncenteredValY(int,int,int) const;
	float uncenteredValZ(int,int,int) const;
	float& uncenteredValXRef(int,int,int);
	float& uncenteredValYRef(int,int,int);
	float& uncenteredValZRef(int,int,int);

	float valX(int) const;
	float valY(int) const;
	float valZ(int) const;
	float& valXRef(int);
	float& valYRef(int);
	float& valZRef(int);

	void trimField(int,int,int,VectorField3D&);
        float giveField(int);

    protected:
	float* f;
	int sizeX, sizeY, sizeZ;
	int centerX, centerY, centerZ;
	int index(int,int,int) const;
};

//CenterX is defined as sizeX/2, likewise for Y and Z
inline int VectorField3D::nx() const
{
    return centerX;
}

inline int VectorField3D::ny() const
{
    return centerY;
}

inline int VectorField3D::nz() const
{
    return centerZ;
}

//Sizes of the dimensions
inline int VectorField3D::numX() const
{
    return sizeX;
}

inline int VectorField3D::numY() const
{
    return sizeY;
}

inline int VectorField3D::numZ() const
{
    return sizeZ;
}

//For a point (i,j,k), index() returns the integer valued index of that point
//in the one dimensional array
inline int VectorField3D::index(int i, int j, int k) const
{
    assert(0 <= i && i <= sizeX &&
	   0 <= j && j <= sizeY &&
	   0 <= k && k <= sizeZ);
    return 3 * ( k + sizeZ * (j + sizeY  * i));
}

//Returns X value at 1-D index i
inline float VectorField3D::valX(int i) const
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i];
}

//Returns Y value at 1-D index i
inline float VectorField3D::valY(int i) const
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i + 1];
}

//Returns Z value at 1-D index i
inline float VectorField3D::valZ(int i) const
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i + 2];
}

//Returns references to the X,Y,Z values at i
inline float& VectorField3D::valXRef(int i)
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i];
}

inline float& VectorField3D::valYRef(int i)
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i + 1];
}

inline float& VectorField3D::valZRef(int i)
{
    assert(0 <= i && i < sizeX * sizeY * sizeZ);
    return f[3*i + 2];
}

//Returns the X value at point (i,j,k)
inline float VectorField3D::valX(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)];
}

//As above, Y value
inline float VectorField3D::valY(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)+1];
}

//As above, Z value
inline float VectorField3D::valZ(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)+2];
}

//As above, X reference
inline float& VectorField3D::valXRef(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)];
}

//As above, Y reference
inline float& VectorField3D::valYRef(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)+1];
}

//As above, Z reference
inline float& VectorField3D::valZRef(int i, int j, int k) const
{
    assert(-centerX <= i && i <= centerX &&
	   -centerY <= j && j <= centerY &&
	   -centerZ <= k && k <= centerZ);
   return f[index(i+centerX, j+centerY, k+centerZ)+2];
}

//Returns the value at i
inline float VectorField3D::giveField(int i)
{
    return *(f + i);
}
#endif