/***************************************************************************
                          Mat4d.h  -  description
                             -------------------
    begin                : Thu Sept 17 2003
    copyright            : (C) 2003 by Brian J. Keay
    email                : ulisse2k@yahoo.com
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef _4D_TYPES_H_
#define _4D_TYPES_H_

//#include <complex>
//#include <fstream>
//#include "nrutil_nr.h"
using namespace std;


template <class T>
class NRMat4d {
private:
	int nn;
	int mm;
	int kk;
	int qq;
	T ****v;
public:
	NRMat4d();
	NRMat4d(int n, int m, int k, int q);
	inline T*** operator[](const int i);	//subscripting: pointer to row i
	inline const T* const * const * operator[](const int i) const;		// added const *
	inline int dim1() const;
	inline int dim2() const;
	inline int dim3() const;
	inline int dim4() const;
	~NRMat4d();
};

template <class T>
NRMat4d<T>::NRMat4d(): nn(0), mm(0), kk(0), qq(0), v(0) {}

template <class T>
NRMat4d<T>::NRMat4d(int n, int m, int k, int q) : nn(n), mm(m), kk(k), qq(q), v(new T***[n])
{
	int i,j,r;
	v[0] = new T**[n*m];
	v[0][0] = new T*[n*m*k];
	v[0][0][0] = new T[n*m*k*q];

	for(r=1; r<k; r++)
				v[0][0][r] = v[0][0][r-1] + q;
	for(j=1; j<m; j++){
				v[0][j] = v[0][j-1] + k;
				v[0][j][0] = v[0][j-1][0] + k*q;
				}
	for(i=1; i<n; i++) {
		v[i] = v[i-1] + m;
		v[i][0] = v[i-1][0] + m*k;
		v[i][0][0] = v[i-1][0][0] + m*k*q;
		for(j=1; j<m; j++){
			v[i][j] = v[i][j-1] + k;			
			v[i][j][0] = v[i][j-1][0] + k*q;
			for(r=1; r<k; r++)
				v[i][j][r] = v[i][j][r-1] + q;
	}}
for(j=1; j<m; j++)
for(r=1; r<k; r++) v[0][j][r] = v[0][j][r-1] + q;

for(i=1; i<n; i++)
for(r=1; r<k; r++) v[i][0][r] = v[i][0][r-1] + q;

}	

template <class T>
inline T*** NRMat4d<T>::operator[](const int i) //subscripting: pointer to row i
{
	return v[i];
}


template <class T>
inline const T* const * const * NRMat4d<T>::operator[](const int i) const           // added const *
{
	return v[i];
}

/**/

template <class T>
inline int NRMat4d<T>::dim1() const
{
	return nn;
}

template <class T>
inline int NRMat4d<T>::dim2() const
{
	return mm;
}

template <class T>
inline int NRMat4d<T>::dim3() const
{
	return kk;
}

template <class T>
inline int NRMat4d<T>::dim4() const
{
	return qq;
}


template <class T>
NRMat4d<T>::~NRMat4d()
{
	if (v != 0) {
		delete[] (v[0][0][0]);
		delete[] (v[0][0]);
		delete[] (v[0]);
		delete[] (v);
	}
}


// 4D Matrix Types

typedef const NRMat4d<DP> Mat4D_I_DP;
typedef NRMat4d<DP> Mat4D_DP, Mat4D_O_DP, Mat4D_IO_DP;

//////////////////////////////////////////////

////////////////////////////////////////////

#endif /* _NR_TYPES_H_ */