// This file is part of the brownmove simulation package
// 
// (C) 2003 - 2009 Tihamer Geyer
// tihamer.geyer@bioinformatik.uni-saarland.de
// 
// Check for Updates at the Project Homepage 
// http://service.bioinformatik.uni-saarland.de/brownmove

#include "interactingBody.h"
#include "molecule.h"

#include "hishape.h"


namespace brown{
	
	void (*HiShape::tensor) (diffusionMatrix& D, const doublevec &myPos, const double myR, const doublevec &otherPos, const double otherR, bool doesRotate) = NULL;
	
	HiShape::HiShape(const HiShape &other){
		doublevec_set(myPosition, other.myPosition);
		radius = other.radius;
		HICutOff = other.HICutOff;
	}
	
	HiShape::~HiShape(){
	}
	
	
	std::string HiShape::describe(){
		ostringstream ausgabe("", ios::out);
		int i=0;
		ausgabe << "radius = "     << radius;
		return(ausgabe.str());
	}
	
	
	
	void HiShape::recalculate(const fullPosition &p) {
		doublevec_set(myPosition, p.pos);
	}
	
	void HiShape::initForces(const fullForce &rF, const fullForce &eF)
	{
// 		std::cerr << "HiShape::initForces" << std::endl;
		
		doublevec_set(rawRandomForce.f,   rF.f);
		doublevec_set(rawRandomForce.t,   rF.t);
		doublevec_set(rawExternalForce.f, eF.f);
		doublevec_set(rawExternalForce.t, eF.t);
		
		// reset effective forces
		doublevec_set(effRandomForce.f,   0.0);
		doublevec_set(effRandomForce.t,   0.0);
		doublevec_set(effExternalForce.f, 0.0);
		doublevec_set(effExternalForce.t, 0.0);
		
		// also reset HI corrections factors
		for(int i = 0; i<6; i++) {
			Ci_sum[i] = 0.0;
		}
	}
	

	double HiShape::getHiCorrections(HiShape *otherHiS, bool includeRotation) 
	{
		diffusionMatrix coupling;
		int i, k;
		fullForce otherForce;
		fullForce myContrib, otherContrib;
		
// 		if(includeRotation) {
// 			std::cerr << "HiShape::getHiCorrections with Rotation" << std::endl;
// 		} else {
// 			std::cerr << "HiShape::getHiCorrections - no Rotation" << std::endl;
// 		}
		
		// determine distance and compare to cut-off
		doublevec abstand;
		doublevec_set(abstand, myPosition);
		doublevec_sub(abstand, otherHiS->getPosition());
		double wieWeit = doublevec_length(abstand);
		
		// std::cerr << "HiShape::getHiCorrections: distance = " << wieWeit << std::endl;
		
		if(wieWeit <= HICutOff) 
		{
			// std::cerr << "HiShape::getHiCorrections:inside cutoff radius "<< HICutOff << " - determining HI" << std::endl;
			
			// set up 6x6 HI-Tensor with includeRotation
			double otherRadius = otherHiS->getRadius();
			tensor(coupling, myPosition, radius, otherHiS->getPosition(), otherRadius, includeRotation);
			
			
			// calculate HI corrections for external forces
			doublevec_set(otherForce.f, otherHiS->getRawExternalForce().f);
			doublevec_set(otherForce.t, otherHiS->getRawExternalForce().t);
			
			transform(myContrib.f, myContrib.t, coupling, otherForce.f, otherForce.t);						// other acting on me
			transform_T(otherContrib.f, otherContrib.t, coupling, rawExternalForce.f, rawExternalForce.t);	// me onto other
			
			// 		std::cerr << "HiShape::getHiCorrections: otherForce.f       = " <<otherForce.f[0]<<", "<<otherForce.f[1]<<", "<<otherForce.f[2]<< std::endl;
			// 		std::cerr << "HiShape::getHiCorrections: rawExternalForce.f = " <<rawExternalForce.f[0]<<", "<<rawExternalForce.f[1]<<", "<<rawExternalForce.f[2]<< std::endl;
			// 		std::cerr << "HiShape::getHiCorrections: myContrib.f        = " <<myContrib.f[0]<<", "<<myContrib.f[1]<<", "<<myContrib.f[2]<< std::endl;
			// 		std::cerr << "HiShape::getHiCorrections: otherContrib.f     = " <<otherContrib.f[0]<<", "<<otherContrib.f[1]<<", "<<otherContrib.f[2]<< std::endl;
			
			// correct for radii, which are not considered in the Rotne-Prager Tensor
			doublevec_mult(myContrib.f, radius);
			doublevec_mult(myContrib.t, radius);
			doublevec_mult(otherContrib.f, otherRadius);
			doublevec_mult(otherContrib.t, otherRadius);

			addEffExternal(myContrib);
			otherHiS->addEffExternal(otherContrib);
			
			// 		std::cerr << "HiShape::getHiCorrections: effExternalForce.f = " <<effExternalForce.f[0]<<", "<<effExternalForce.f[1]<<", "<<effExternalForce.f[2]<< std::endl;
			
			// same thing for random contributions
			doublevec_set(otherForce.f, otherHiS->getRawRandomForce().f);
			doublevec_set(otherForce.t, otherHiS->getRawRandomForce().t);
			
			transform(myContrib.f, myContrib.t, coupling, otherForce.f, otherForce.t);					// other acting on me
			transform_T(otherContrib.f, otherContrib.t, coupling, rawRandomForce.f, rawRandomForce.t);	// me onto other
			
			doublevec_mult(myContrib.f, radius);
			doublevec_mult(myContrib.t, radius);
			doublevec_mult(otherContrib.f, otherRadius);
			doublevec_mult(otherContrib.t, otherRadius);
			
			addEffRandom(myContrib);
			otherHiS->addEffRandom(otherContrib);
			
			
			// correction factors: epsilon
			double eps_sum = 0.0;
			double myCi[6], otherCi[6];
			double Dik, Dii;

			
			for(i=0; i<6; i++) {
				myCi[i] = 0.0;
				otherCi[i] = 0.0;
			}
			
			#ifdef HI_WITH_ROTATION
			for(i=0; i<6; i++) {
				for(k=i+1; k<6; k++) 
			#else
			for(i=0; i<3; i++) {			// ### Hi-Change ###
				for(k=0; k<3; k++)			// ### Hi-Change ###
			#endif
				{
					Dik = coupling[i][k];
					eps_sum += Dik;
					myCi[i] += Dik*Dik;
					otherCi[k] += Dik*Dik;
					
					// std::cerr << "HiShape::getHiCorrections:  myCi    = " << myCi[0]<<", "<<myCi[1]<<", "<<myCi[2]<<" - "<<myCi[3]<<", "<<myCi[4]<<", "<<myCi[5] << std::endl;
					// std::cerr << "HiShape::getHiCorrections:  otherCi = " << otherCi[0]<<", "<<otherCi[1]<<", "<<otherCi[2]<<" - "<<otherCi[3]<<", "<<otherCi[4]<<", "<<otherCi[5] << std::endl;
				}
			}

			eps_sum *= (0.5*(radius + otherRadius));
			double radiusProduct = radius*otherRadius;
			for(i=0; i<6; i++) {
				myCi[i] *= radiusProduct;
				otherCi[i] *= radiusProduct;
			}
			
			// std::cerr << "HiShape::getHiCorrections: eps_sum = " << eps_sum << std::endl;
			
			addCorrection(myCi);
			otherHiS->addCorrection(otherCi);
			
			return(eps_sum);
		}	
		else // distance > HI cutoff
		{
			return(0.0);
		}
	}
	
	
	void HiShape::addCorrection(const double cs[]) {
		for(int i=0; i<6; i++) {
			Ci_sum[i] += cs[i];
		}
	}
	


	void HiShape::doEffectiveForces(fullForce &totalForce, const double &beta_prime)
	{
		// Ci
		double localCi[6];
		double help;
		double beta2 = beta_prime*beta_prime;
		int i;
		
		for(i=0; i<6; i++) {
			help = 1.0 + beta2 * Ci_sum[i];
			localCi[i] = 1.0 / sqrt(help);
		}
		
// 		std::cerr << "HiShape::doEffectiveForces: beta_prime = " << beta_prime << "   Ci = " << localCi[0]<<", "<<localCi[1]<<", "<<localCi[2]<<" - "<<localCi[3]<<", "<<localCi[4]<<", "<<localCi[5] << std::endl;
		
// 		std::cerr << "HiShape::doEffectiveForces: rawExternalForce.f = " <<rawExternalForce.f[0]<<", "<<rawExternalForce.f[1]<<", "<<rawExternalForce.f[2]<< std::endl;
// 		std::cerr << "HiShape::doEffectiveForces: effExternalForce.f = " <<effExternalForce.f[0]<<", "<<effExternalForce.f[1]<<", "<<effExternalForce.f[2]<< std::endl;
// 		std::cerr << "HiShape::doEffectiveForces: rawRandomForce.f   = " <<rawRandomForce.f[0]<<", "<<rawRandomForce.f[1]<<", "<<rawRandomForce.f[2]<< std::endl;
// 		std::cerr << "HiShape::doEffectiveForces: effRandomForce.f   = " <<effRandomForce.f[0]<<", "<<effRandomForce.f[1]<<", "<<effRandomForce.f[2]<< std::endl;

		// start with random forces
		doublevec_set(totalForce.f, effRandomForce.f);
		doublevec_set(totalForce.t, effRandomForce.t);
		doublevec_mult(totalForce.f, beta_prime);
		doublevec_mult(totalForce.t, beta_prime);
		doublevec_add(totalForce.f, rawRandomForce.f);
		doublevec_add(totalForce.t, rawRandomForce.t);
		
		for(i=0; i<3; i++) {
			totalForce.f[i] *= localCi[i];
			totalForce.t[i] *= localCi[i+3];
		}
		
// 		std::cerr << "HiShape::doEffectiveForces: totalForce.f V1    = " <<totalForce.f[0]<<", "<<totalForce.f[1]<<", "<<totalForce.f[2]<< std::endl;

		// add external forces
		doublevec_add(totalForce.f, effExternalForce.f);
		doublevec_add(totalForce.t, effExternalForce.t);
		doublevec_add(totalForce.f, rawExternalForce.f);
		doublevec_add(totalForce.t, rawExternalForce.t);

// 		std::cerr << "HiShape::doEffectiveForces: totalForce.f V2    = " <<totalForce.f[0]<<", "<<totalForce.f[1]<<", "<<totalForce.f[2]<< std::endl;
}
	 

	

	
	void HiShape::setupDyadMat(const doublevec& dist, tripleMatrix& dyad )
	{
		for (int i=0;i<3;++i) {
			for (int j=0; j<3; ++j) {
				dyad[i][j] = dist[i] * dist[j];
			}
		}
	}
	

	
	//! This RotnePrager implementation does not take different sized molecules into account. 
	//! A Way to include these can be found in "Hydrodynamic Properties of Macromolecular Complexes", 
	//! Jose Garcia de la Torre, Biopolymeres 16, 1747 - 1763 (1977)
	void HiShape::RotnePrager(diffusionMatrix& mat, const doublevec &p_i, const double radius, const doublevec &p_j, const double otherRadius, bool doesRotate)
	{
		doublevec distvec;
		doublevec tmprow;
		tripleMatrix dyadmat,tmp1,tmp2,tmp3;
		double distance;
		// In VARIABLES _ indicates an over
		// one is o
		// numbers denote the power
		double o_rij;	// 1/r
		double o_rij2;  // 1/r^2
		double o_rij3;  // 1/r^3
		double ai2aj2;	// ai^2 + aj^2

		//setup the normalized distance vector and scalar and distance dependent constants
		doublevec_set(distvec, p_j);
		doublevec_sub(distvec, p_i);
		distance = doublevec_length(distvec);
		doublevec_mult(distvec, 1.0/distance);	// normalized distance
		
		o_rij   = 1.0 / distance;
		o_rij2  = o_rij / distance;
		o_rij3  = o_rij2 / distance;
		// a_rij3  = a_rij2 / distance;

		ai2aj2 = radius*radius + otherRadius*otherRadius;
		
		//setup the dyadic product
		setupDyadMat(distvec,dyadmat); 
		
		// Fill the diffusion matrix
		
		// distance > 2 radii
		// if(distance >= (radius + otherRadius))
		// {
			// Trans-Trans-Part: 3/4rij (1 + rij*rij) + 1/4 (ai^2+aj^2)/rij^3 ( 1 - 3rij*rij)
			//                         tmp1                        tmp2
			//                         
			trMat_set(tmp1,1.0); 			// set identity
			trMat_add(tmp1,dyadmat);		// add dyadic product
			trMat_mult(tmp1,0.75/distance);	// prefactor
			
			trMat_set(tmp2, 1.0); 			// set identity
			trMat_set(tmp3, dyadmat);		// tmp3 = 3 (r_ij * r_ij)
			trMat_mult(tmp3, 3.0);
			trMat_sub(tmp2, tmp3);			// (1 - 3 rij*rij)
			trMat_mult(tmp2, 0.25*ai2aj2*o_rij3);	// tmp2 = 1/4 (ai^2+aj^2)/rij^3 (I - 3 (r_ij*r_ij))
			
			trMat_add(tmp1,tmp2);			// tmp1 = 3/4r_ij * ( I + r_ij*r_ij ) + 1/4 (ai^2+aj^2)/rij^3 (I - 3 (r_ij*r_ij))
		// }
		//done, set the translation part of the diffusion Matrix
		mat.set(tmp1,0);
		
		if(doesRotate)
		{
			// Rotation-Rotation: 3/8rij^3 (3 rij*rij - 1)
			trMat_set(tmp1,1.0); 			// set identity
			trMat_set(tmp2, dyadmat);		// 3 rij*rij
			trMat_mult(tmp2, 3.0);
			trMat_sub(tmp2, tmp1);			// 3 rij*rij - 1
			trMat_mult(tmp2, 0.375*o_rij3);	// prefactor 3/8rij^3
			
			mat.set(tmp2, 3);
			
			// Trans-Rot = Rot-Trans^T: -3/4r^2 eijk rij
			// tmp1 = epsilon r_ij  (don't forget the minus signs!!!)
			tmp1[0][0] = 0.0        ; tmp1[0][1] = -distvec[2]; tmp1[0][2] = distvec[1] ;
			tmp1[1][0] = distvec[2] ; tmp1[1][1] = 0.0        ; tmp1[1][2] = -distvec[0];
			tmp1[2][0] = -distvec[1]; tmp1[2][1] = distvec[0] ; tmp1[2][2] = 0.0        ;
			
			// trMat_mult(tmp1, (-0.75 * a_rij2));		// -0.75 * a/rij^2 * epsilon r_ij
			trMat_mult(tmp1, (0.75 * o_rij2));		// -0.75 * o/rij^2 * epsilon r_ij
			
			// set the Trans-Rot part of the diffusionMatrix
			mat.set(tmp1,1);
			
			// RT = TR^T -> multiply by (-1) and set the Rot-Trans part of the diffusionMatrix
			// trMat_mult(tmp1, -1.0);
			mat.set(tmp1,2);
		}
		else
		{
			trMat_set(tmp1, 0.0);	// set 3x3 Null-matrix
			mat.set(tmp1, 1);		// and zero the rest of the diffusion matrix
			mat.set(tmp1, 2);
			mat.set(tmp1, 3);
		}
		
		
		// debug output
// 		std::cerr << "RotnePrager: "<<mat[0][0]<<"\t"<<mat[0][1]<<"\t"<<mat[0][2]<<"\t"<<mat[0][3]<<"\t"<<mat[0][4]<<"\t"<<mat[0][5]<<std::endl;
// 		std::cerr << "             "<<mat[1][0]<<"\t"<<mat[1][1]<<"\t"<<mat[1][2]<<"\t"<<mat[1][3]<<"\t"<<mat[1][4]<<"\t"<<mat[1][5]<<std::endl;
// 		std::cerr << "             "<<mat[2][0]<<"\t"<<mat[2][1]<<"\t"<<mat[2][2]<<"\t"<<mat[2][3]<<"\t"<<mat[2][4]<<"\t"<<mat[2][5]<<std::endl;
// 		std::cerr << "             "<<mat[3][0]<<"\t"<<mat[3][1]<<"\t"<<mat[3][2]<<"\t"<<mat[3][3]<<"\t"<<mat[3][4]<<"\t"<<mat[3][5]<<std::endl;
// 		std::cerr << "             "<<mat[4][0]<<"\t"<<mat[4][1]<<"\t"<<mat[4][2]<<"\t"<<mat[4][3]<<"\t"<<mat[4][4]<<"\t"<<mat[4][5]<<std::endl;
// 		std::cerr << "             "<<mat[5][0]<<"\t"<<mat[5][1]<<"\t"<<mat[5][2]<<"\t"<<mat[5][3]<<"\t"<<mat[5][4]<<"\t"<<mat[5][5]<<std::endl;
	}
	
	

//	// *** Below is old code left for reference ***
//	// 	
// 	//! This RotnePrager implementation does not take different sized molecules into account. 
// 	//! A Way to include these can be found in "Hydrodynamic Properties of Macromolecular Complexes", 
// 	//! Jose Garcia de la Torre, Biopolymeres 16, 1747 - 1763 (1977)
// 	void HiShape::RotnePrager(diffusionMatrix& mat, const doublevec &p_i, const double radius, const doublevec &p_j, const double otherRadius, bool doesRotate)
// 	{
// 		doublevec distvec;
// 		doublevec tmprow;
// 		tripleMatrix dyadmat,tmp1,tmp2,tmp3;
// 		double distance;
// 		// In VARIABLES _ indicates an over
// 		// one is o
// 		// numbers denote the power
// 		double a_rij;	// a/r
// 		double a3_rij3; // a^3/r^3
// 		double a_rij2;  // a/r^2
// 		double a_rij3;  // a/r^3
// 		
// 		//setup the normalized distance vector and scalar and distance dependent constants
// 		doublevec_set(distvec, p_j);
// 		doublevec_sub(distvec, p_i);
// 		distance = doublevec_length(distvec);
// 		doublevec_mult(distvec,1.0/distance);
// 		
// 		a_rij   = radius / distance;
// 		a3_rij3 = a_rij * a_rij * a_rij;
// 		a_rij2  = a_rij / distance;
// 		a_rij3  = a_rij2 / distance;
// 		
// 		//setup the dyadic product
// 		setupDyadMat(distvec,dyadmat); 
// 		
// 		// Fill the diffusion matrix
// 		
// 		// distance > 2 radii
// 		if(distance >= (radius + otherRadius))
// 		{
// 			// Trans-Trans-Part: 3a/4rij (1 + rij*rij) + 1/2 (a/rij)^3 ( 1 - 3rij*rij)
// 			//                         tmp1                        tmp2
// 			//                         
// 			trMat_set(tmp1,1.0); 			// set identity
// 			trMat_add(tmp1,dyadmat);		// add dyadic product
// 			trMat_mult(tmp1,0.75*a_rij);	// prefactor
// 			
// 			trMat_set(tmp2, 1.0); 			// set identity
// 			trMat_set(tmp3, dyadmat);		// tmp3 = 3 (r_ij * r_ij)
// 			trMat_mult(tmp3, 3.0);
// 			trMat_sub(tmp2, tmp3);			// (1 - 3 rij*rij)
// 			trMat_mult(tmp2, 0.5*a3_rij3);	// tmp2 = 1/2 (a/r_ij)^3 (I - 3 (r_ij*r_ij))
// 			
// 			trMat_add(tmp1,tmp2);			// tmp1 = 3a/4r_ij * ( I + r_ij*r_ij ) + 1/2 (a/r_ij)^3 (I - 3 (r_ij*r_ij))
// 		}
// 		else
// 		{
// 			// Trans-Trans: (1 - 9rij/32a) * 1  +  3rij/32a * (rij*rij)
// 			//                      tmp1                 tmp2
// 			//                      
// 			trMat_set(tmp1, 1.0);					// set identity
// 			trMat_mult(tmp1, 1.0 - 9/(32*a_rij));	// prefactor
// 			
// 			trMat_set(tmp2, dyadmat);				// dyadic product for second term
// 			trMat_mult(tmp2, 3.0/(32*a_rij));		// prefactor
// 			
// 			trMat_add(tmp1, tmp2);					// sum up
// 		}
// 		//done, set the translation part of the diffusion Matrix
// 		mat.set(tmp1,0);
// 		
// 		if(doesRotate)
// 		{
// 			// Rotation-Rotation: 3a/8rij^3 (3 rij*rij - 1)
// 			trMat_set(tmp1,1.0); 			// set identity
// 			trMat_set(tmp2, dyadmat);		// 3 rij*rij
// 			trMat_mult(tmp2, 3.0);
// 			trMat_sub(tmp2, tmp1);			// 3 rij*rij - 1
// 			trMat_mult(tmp2, 0.375*a_rij3);	// prefactor 3a/8rij^3
// 			
// 			mat.set(tmp2, 3);
// 			
// 			// Trans-Rot = Rot-Trans^T: -3a/4r^2 eijk rij
// 			// tmp1 = epsilon r_ij  (don't forget the minus signs!!!)
// 			tmp1[0][0] = 0.0        ; tmp1[0][1] = -distvec[2]; tmp1[0][2] = distvec[1] ;
// 			tmp1[1][0] = distvec[2] ; tmp1[1][1] = 0.0        ; tmp1[1][2] = -distvec[0];
// 			tmp1[2][0] = -distvec[1]; tmp1[2][1] = distvec[0] ; tmp1[2][2] = 0.0        ;
// 			
// 			// trMat_mult(tmp1, (-0.75 * a_rij2));		// -0.75 * a/rij^2 * epsilon r_ij
// 			trMat_mult(tmp1, (0.75 * a_rij2));		// -0.75 * a/rij^2 * epsilon r_ij
// 			
// 			// set the Trans-Rot part of the diffusionMatrix
// 			mat.set(tmp1,1);
// 			
// 			// RT = TR^T -> multiply by (-1) and set the Rot-Trans part of the diffusionMatrix
// 			// trMat_mult(tmp1, -1.0);
// 			mat.set(tmp1,2);
// 		}
// 		else
// 		{
// 			trMat_set(tmp1, 0.0);	// set 3x3 Null-matrix
// 			mat.set(tmp1, 1);		// and zero the rest of the diffusion matrix
// 			mat.set(tmp1, 2);
// 			mat.set(tmp1, 3);
// 		}
// 		
// 		
// 		// debug output
// 		std::cerr << "RotnePrager: "<<mat[0][0]<<"\t"<<mat[0][1]<<"\t"<<mat[0][2]<<"\t"<<mat[0][3]<<"\t"<<mat[0][4]<<"\t"<<mat[0][5]<<std::endl;
// 		std::cerr << "             "<<mat[1][0]<<"\t"<<mat[1][1]<<"\t"<<mat[1][2]<<"\t"<<mat[1][3]<<"\t"<<mat[1][4]<<"\t"<<mat[1][5]<<std::endl;
// 		std::cerr << "             "<<mat[2][0]<<"\t"<<mat[2][1]<<"\t"<<mat[2][2]<<"\t"<<mat[2][3]<<"\t"<<mat[2][4]<<"\t"<<mat[2][5]<<std::endl;
// 		std::cerr << "             "<<mat[3][0]<<"\t"<<mat[3][1]<<"\t"<<mat[3][2]<<"\t"<<mat[3][3]<<"\t"<<mat[3][4]<<"\t"<<mat[3][5]<<std::endl;
// 		std::cerr << "             "<<mat[4][0]<<"\t"<<mat[4][1]<<"\t"<<mat[4][2]<<"\t"<<mat[4][3]<<"\t"<<mat[4][4]<<"\t"<<mat[4][5]<<std::endl;
// 		std::cerr << "             "<<mat[5][0]<<"\t"<<mat[5][1]<<"\t"<<mat[5][2]<<"\t"<<mat[5][3]<<"\t"<<mat[5][4]<<"\t"<<mat[5][5]<<std::endl;
// 	}
	
	
	void HiShape::NoHI(diffusionMatrix& mat, const doublevec &myPos, const double myR, const doublevec &otherPos, const double otherR, bool doesRotate)
	{
		makeDiagonal(mat, 0.0);
		
// 		std::cerr << "NoHI       : "<<mat[0][0]<<"\t"<<mat[0][1]<<"\t"<<mat[0][2]<<"\t"<<mat[0][3]<<"\t"<<mat[0][4]<<"\t"<<mat[0][5]<<std::endl;
// 		std::cerr << "             "<<mat[1][0]<<"\t"<<mat[1][1]<<"\t"<<mat[1][2]<<"\t"<<mat[1][3]<<"\t"<<mat[1][4]<<"\t"<<mat[1][5]<<std::endl;
// 		std::cerr << "             "<<mat[2][0]<<"\t"<<mat[2][1]<<"\t"<<mat[2][2]<<"\t"<<mat[2][3]<<"\t"<<mat[2][4]<<"\t"<<mat[2][5]<<std::endl;
// 		std::cerr << "             "<<mat[3][0]<<"\t"<<mat[3][1]<<"\t"<<mat[3][2]<<"\t"<<mat[3][3]<<"\t"<<mat[3][4]<<"\t"<<mat[3][5]<<std::endl;
// 		std::cerr << "             "<<mat[4][0]<<"\t"<<mat[4][1]<<"\t"<<mat[4][2]<<"\t"<<mat[4][3]<<"\t"<<mat[4][4]<<"\t"<<mat[4][5]<<std::endl;
// 		std::cerr << "             "<<mat[5][0]<<"\t"<<mat[5][1]<<"\t"<<mat[5][2]<<"\t"<<mat[5][3]<<"\t"<<mat[5][4]<<"\t"<<mat[5][5]<<std::endl;
	}

	
	void HiShape::dumpStructure() {
		std::cout << "\n\t\tHiShape myOneAndOnlySphere" << std::endl;
		std::cout << "\t\t\t"<< radius << std::endl;
		std::cout << "\t\tEnd" << std::endl;
	}
	
}