#ifndef _FORTRANARRAY_H_
#define _FORTRANARRAY_H_

/**
	@file	FORTRANarray.hpp
	@brief   A templated class to simplify multi-dimensional array interfaces with FORTRAN methods.

	@date $Date: 2007/07/10 03:46:55 $

	@version $Rev$

	@author
		- Lance Deaver

	@copyright
		(©) Copyright 2006 by GATS Inc.
		11864 Canon Blvd., Suite 101, Newport News, VA 23606

		All Rights Reserved. No part of this software or publication may be
		reproduced, stored in a retrieval system, or transmitted, in any form
		or by any means, electronic, mechanical, photocopying, recording, or
		otherwise without the prior written permission of GATS Inc.

	@note	Only allows up to 4 dimensional arrays.
	@bug		None known
	@todo  Uses asserts to check for array overflow..  would exceptions be better?
	
	<em> $Id: FORTRANarray.hpp,v 1.1 2007/07/10 03:46:55 deaver Exp deaver $ </em>

*/
#include <vector>
#include <cassert>

template <class T>
class FORTRANarray
{
private:
	union {
		unsigned N_[4];
		struct {
			unsigned N1_;
			unsigned N2_;
			unsigned N3_;
			unsigned N4_;
		};
	};
	unsigned _numberOfDimensions; 
	std::vector<T> rep_;
public:
	FORTRANarray():
		N1_(0),N2_(0),N3_(0),N4_(0),_numberOfDimensions(0) { rep_.clear(); }
	FORTRANarray(unsigned n1):
		N1_(n1),N2_(0),N3_(0),N4_(0),_numberOfDimensions(1)
	    { assert(n1>0); rep_.resize(n1); }
	FORTRANarray(unsigned n1, unsigned n2):
		N1_(n1),N2_(n2),N3_(0),N4_(0),_numberOfDimensions(2)
		{assert( n1>0 && n2>0); rep_.resize(n1*n2); }
	FORTRANarray(unsigned n1, unsigned n2, unsigned n3):
		N1_(n1),N2_(n2),N3_(n3),N4_(0),_numberOfDimensions(3)
		{ assert(n1>0 && n2>0 && n3> 0) ; rep_.resize(n1*n2*n3); }
	FORTRANarray(unsigned n1, unsigned n2, unsigned n3, unsigned n4):
		N1_(n1),N2_(n2),N3_(n3),N4_(n4),_numberOfDimensions(4)
		{ assert(n1>0 && n2>0 && n3>0 && n4>0); rep_.resize(n1*n2*n3*n4);  }
		

	void setVal(const T& val,  const unsigned i)
	{ 
		assert(i <= N1_ && i > 0); 
		rep_[i-1] = val; 
	}
	void setVal(const T& val, const unsigned i, const unsigned j) 
	{ 
		assert(i <= N1_ && j <= N2_ && i > 0 && j > 0); 
		rep_[i-1 + N1_*(j-1) ] = val;
	}
	void setVal(const T& val, const unsigned i, const unsigned j, const unsigned k) 
	{
		assert(i <= N1_ && j <= N2_ && k <= N3_ && i > 0 && j > 0 && k > 0); 
		rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1) ] = val; 
	} 
	void setVal(const T& val, const unsigned i, const unsigned j, const  unsigned k, const unsigned l) 
	{
		assert(i <= N1_ && j <= N2_ && k <= N3_ && l <= N4_ && i > 0 && j > 0 && k > 0 && l > 0); 
		rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1)+ N1_*N2_*N3_*(l-1) ] = val; 
	}  
	

	T getVal(const unsigned i=1) const
	{
		assert(i <= N1_ && i > 0);
		return rep_[i-1];
	} 
	T getVal(const unsigned i, const unsigned j) const 
	{
		assert(i <= N1_ && j <= N2_ && i > 0 && j > 0); 
		return rep_[i-1+N1_*(j-1)]; 
	} 
	T getVal(const unsigned i, const unsigned j, const unsigned k) const 
	{ 
		assert(i <= N1_ && j <= N2_ && k <= N3_ && i > 0 && j >0 && k >0); 
		return rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1) ]; 
	}
	T getVal(const unsigned i, const unsigned j,const unsigned k,const unsigned l) const
	{
		assert(i <= N1_ && j <= N2_ && k <= N3_ && l <= N4_ && i > 0 && j >0 && k>0 && l>0);
		return rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1)+ N1_*N2_*N3_*(l-1) ] ; 
	}

	T operator() (const unsigned i) const 
		{ return getVal(i); }
	T operator() (const unsigned i, const unsigned j) const 
		{ return getVal(i,j); }
	T operator() (const unsigned i, const unsigned j, const unsigned k) const 
		{ return getVal(i,j,k); }
	T operator() (const unsigned i, const unsigned j,const unsigned k,const unsigned l) const 
		{ return getVal(i,j,k,l); }

	T* getArray(const unsigned i=1)
	{
		assert(i <= N1_ && i > 0);
		return &rep_[i-1]; 
	} 
	T* getArray(const unsigned i, const unsigned j) 
	{
		assert(i <= N1_ && j <= N2_ && i > 0 && j > 0); 
		return &rep_[i-1+N1_*(j-1)]; 
	} 
	T* getArray(const unsigned i, const unsigned j, const unsigned k) 
	{ 
		assert(i <= N1_ && j <= N2_ && k <= N3_ && i > 0 && j >0 && k >0); 
		return &rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1) ]; 
	}
	T* getArray(const unsigned i, const unsigned j,const unsigned k,const unsigned l)
	{
		assert(i <= N1_ && j <= N2_ && k <= N3_ && l <= N4_ && i > 0 && j >0 && k>0 && l>0);
		return &rep_[i-1+N1_*(j-1) + N1_*N2_*(k-1)+ N1_*N2_*N3_*(l-1) ] ; 
	}
	unsigned getNdims() const { return _numberOfDimensions; }
    unsigned getDims(const unsigned i) const {
    		assert( i>=1 && i <= _numberOfDimensions);
    		return N_[i-1];
    	}
    	unsigned size() const { return rep_.size() ; }

	~FORTRANarray() {  }
};


#endif //_FORTRANARRAY_H_