// 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 "vdwshape.h"

namespace brown{
	
	vdWShape::vdWShape(){
		paramsArray = NULL;
		myWallColor = -1;
		myWallOffset = 0.0;
		doublevec_set(myActualWallPos, 0.0);
		doublevec_set(myActualWallNormal, 1.0, 0.0, 0.0);
	}
	
	vdWShape::vdWShape(const vdWShape &other){
		spheres = other.spheres;
		actual_spheres = other.actual_spheres;
		blocks = other.blocks;
		actual_blocks = other.actual_blocks;
		paramsArray = other.paramsArray;
		myWallColor = other.myWallColor;
		myWallOffset = other.myWallOffset;
		doublevec_set(myActualWallPos, other.myActualWallPos);
		doublevec_set(myActualWallNormal, other.myActualWallNormal);
	}
	
	vdWShape::~vdWShape(){
		spheres.clear();
		actual_spheres.clear();
		blocks.clear();
		actual_blocks.clear();
	}
	
	
	std::string vdWShape::describe(){
		ostringstream ausgabe("", ios::out);
		int i=0;
		
		ausgabe << "actual sphere positions";
		vector<vdWPos>::iterator mlc;
		for(mlc = actual_spheres.begin(); mlc != actual_spheres.end(); ++mlc, i++){
			ausgabe << " : "<<(*mlc).pos[0]<<" "<<(*mlc).pos[1]<<" "<<(*mlc).pos[2] ; 
		}
		// ausgabe <<" :";
		// for(mlc = spheres.begin(); mlc != spheres.end(); ++mlc, i++){
		// 	ausgabe << " : "<<(*mlc).pos[0]<<" "<<(*mlc).pos[1]<<" "<<(*mlc).pos[2] ; 
		// }
		
		if(actual_blocks.size() > 0) {
			ausgabe << "  blocks";
			vector<vdWBlock>::iterator ibl;
			for(ibl = actual_blocks.begin(); ibl != actual_blocks.end(); ibl++) {
				ausgabe << " : " << (*ibl).pos[0]<<" "<<(*ibl).pos[1]<<" "<<(*ibl).pos[2];
			}
		}
				
		if(myWallColor >= 0) {
			ausgabe << "  wall: " <<myActualWallPos[0]<<" "<<myActualWallPos[1]<<" "<<myActualWallPos[2]<<" / "<<myActualWallNormal[0]<<" "<<myActualWallNormal[1]<<" "<<myActualWallNormal[2];
		}
		
		return(ausgabe.str());
	}
	
	
	
	void vdWShape::recalculate(const fullPosition &p) 
	{
		vector<vdWPos>::iterator act = actual_spheres.begin();
		for(vector<vdWPos>::iterator mlc = spheres.begin(); mlc != spheres.end(); ++mlc, ++act)
		{
			doublevec_set((*act).pos, (*mlc).pos);
			doublevec_rotate((*act).pos, p.rot);
			doublevec_set((*act).dr, (*act).pos);
			positionMove((*act).pos, p.pos);
		}

		if(blocks.size() > 0) 
		{
			vector<vdWBlock>::iterator abl = actual_blocks.begin();
			for(vector<vdWBlock>::iterator ibl = blocks.begin(); ibl != blocks.end(); ++ibl, ++abl) 
			{
				// position
				doublevec_set((*abl).pos, (*ibl).pos);
				doublevec_rotate((*abl).pos, p.rot);
				positionMove((*abl).pos, p.pos);
				
				// coordinate system - only rotation
				doublevec_set((*abl).ex, (*ibl).ex);
				doublevec_rotate((*abl).ex, p.rot);
				doublevec_set((*abl).ey, (*ibl).ey);
				doublevec_rotate((*abl).ey, p.rot);
				doublevec_set((*abl).ey, (*ibl).ey);
				doublevec_rotate((*abl).ey, p.rot);
			}
		}
		
		if(myWallColor >= 0) 
		{
			doublevec shift;
			doublevec_set(shift, 1.0, 0.0, 0.0);		// this was the original normal of the wall in the un-teleported gestalt
			doublevec_rotate(shift, p.rot);				// new direction of the normal vector
			doublevec_set(myActualWallNormal, shift);	// save normal
			doublevec_mult(shift, myWallOffset);		// shift distance
			
			doublevec_set(myActualWallPos, p.pos);	// new position of the wall
			doublevec_add(myActualWallPos, shift);	// shift along normal
		}
	}

	

	void vdWShape::addSphere(vdWPos vp, const position &moleculeCenter){
		doublevec_set(vp.dr, vp.pos);
		doublevec_sub(vp.dr, moleculeCenter);
		spheres.push_back(vp);
		actual_spheres.push_back(vp);
	}
	
	
	void vdWShape::addBlock(vdWBlock vb) {
		blocks.push_back(vb);
		actual_blocks.push_back(vb);
	}
	
	
	void vdWShape::calculateImpact(fullForce &myForce, fullForce &otherForce, vdWShape *vs, const bool isPeriodic[3], const doublevec &cellSize)
	{
		position dr, ri, dr_temp;
		double prefactor, B12, R0;
		double r12, dist, rr1, rr3, rr4, rr6; //, rr7;
		force df;
		torque dt;
		long int clr1, clr2;
		
		// first interact all spheres with each other -- nothing happens, if there are none
		// Therefore, begin with the other's spheres (could be a wall..) -- we have at least one vdw sphere, or we would not be here.
		for(vector<vdWPos>::iterator mlc2 = vs->actual_spheres.begin(); mlc2 != vs->actual_spheres.end(); ++mlc2)
		{
			clr2 = (*mlc2).color;
			
			for(vector<vdWPos>::iterator mlc1 = actual_spheres.begin(); mlc1 != actual_spheres.end(); ++mlc1)
			{
				clr1 = (*mlc1).color;
				
				// calculate relative distances etc
				doublevec_set(dr, (*mlc2).pos);
				doublevec_sub(dr, (*mlc1).pos);
				B12 = (*mlc1).r + (*mlc2).r;
				
				double dx = 0.0;
				double dy = 0.0;
				double dz = 0.0;
				int ix = 2, iy = 2, iz = 2;
				
				if(isPeriodic[0]) {
					dx = -cellSize[0];
					ix = 0;
				}
				if(isPeriodic[1]) {
					dy = -cellSize[1];
					iy = 0;
				}
				if(isPeriodic[2]) {
					dz = -cellSize[2];
					iz = 0;
				}
				
				while(true)
				{
					dr_temp[0] = dr[0] + dx;
					dr_temp[1] = dr[1] + dy;
					dr_temp[2] = dr[2] + dz;
					
					r12 = doublevec_length(dr_temp);
					
					dist = r12 - B12;		// empty space between the spheres <- determines the interaction
					
					// check whether dist > cut-off  =>  skip further caluclation
					if((dist <= paramsArray[clr1][clr2].Rmax))
					{
						
						// check for minimal distance for numerical stability
						if(dist < paramsArray[clr1][clr2].Rmin) {
							dist = paramsArray[clr1][clr2].Rmin;
						}
						
						// --- here comes the interaction itself
						dist += paramsArray[clr1][clr2].DR;
						R0 = paramsArray[clr1][clr2].R_0;
						
						rr1 = R0 / dist;		// powers of the inverse distance
						rr3 = rr1 * rr1 * rr1;
						rr4 = rr3 * rr1;
						rr6 = rr3 * rr3;
						// rr7 = rr6 * rr1;
						
						// prefactor = -(12.0 * paramsArray[clr1][clr2].C_12 * rr6 + 6.0 * paramsArray[clr1][clr2].C_6) * rr7 / R0;
						prefactor = -((12.0*paramsArray[clr1][clr2].C_12 * rr6 + 6.0*paramsArray[clr1][clr2].C_6) * rr3 + 3.0*paramsArray[clr1][clr2].C_3) * rr4 / R0;
						
						doublevec_set(df, dr_temp);					// for the direction
						doublevec_mult(df, (prefactor/r12));	// normalize dr by its length
						
						// add force to previous results
						doublevec_add(myForce.f, df);
						doublevec_sub(otherForce.f, df);
						
						// calculate and add torques
						doublevec_kreuzprod(dt, (*mlc1).dr, df);
						doublevec_add(myForce.t, dt);
						doublevec_kreuzprod(dt, (*mlc2).dr, df);
						doublevec_sub(otherForce.t, dt);
					}
					
					// for the periodic boundary conditions
					if(ix < 2) {
						ix++;
						dx += cellSize[0];
					}
					else
					{
						if(isPeriodic[0]) {
							ix = 0;
							dx = -cellSize[0];
						}
						
						if(iy < 2) {
							iy++;
							dy += cellSize[1];
						}
						else
						{
							if(isPeriodic[1]) {
								iy = 0;
								dy = -cellSize[1];
							}
							
							if(iz < 2) {
								iz++;
								dz += cellSize[2];
							}
							else
							{
								break;
							}
						}
					}
				}
			}
		}
		
		
		// calculate interaction with vdw blocks
		// As blocks do not have an interaction off the side walls, there is no reason to apply periodic conditions
		// => blocks are supposed to span the complete simulation volume from wall to wall
		for(vector<vdWBlock>::iterator mlc2 = vs->actual_blocks.begin(); mlc2 != vs->actual_blocks.end(); ++mlc2)
		{
			clr2 = (*mlc2).color;
			doublevec distV;
			double abstand;
			
			for(vector<vdWPos>::iterator mlc1 = actual_spheres.begin(); mlc1 != actual_spheres.end(); ++mlc1)
			{
				clr1 = (*mlc1).color;
				
				// calculate distance and direction, do nothing if outside of bounds
				if((*mlc2).getDistance(distV, abstand, (*mlc1).pos, (*mlc1).r))
				{
					// check whether dist > cut-off  =>  skip further caluclation
					if((abstand <= paramsArray[clr1][clr2].Rmax))
					{
						// check for minimal distance
						if(abstand < paramsArray[clr1][clr2].Rmin) {
							abstand = paramsArray[clr1][clr2].Rmin;
						}
						
						// --- here comes the interaction itself
						abstand += paramsArray[clr1][clr2].DR;
						R0 = paramsArray[clr1][clr2].R_0;
						
						rr1 = R0 / abstand;		// powers of the inverse distance
						rr3 = rr1 * rr1 * rr1;
						rr4 = rr3 * rr1;
						rr6 = rr3 * rr3;
						// rr7 = rr6 * rr1;
						
						// prefactor = (12.0 * paramsArray[clr1][clr2].C_12 * rr6 + 6.0 * paramsArray[clr1][clr2].C_6) * rr7 / R0;
						prefactor = -((12.0*paramsArray[clr1][clr2].C_12 * rr6 + 6.0*paramsArray[clr1][clr2].C_6) * rr3 + 3.0*paramsArray[clr1][clr2].C_3) * rr4 / R0;
						
						// calculate the force normal to the wall (normal is normalized -- should be :-))
						doublevec_set(df, distV);
						doublevec_mult(df, prefactor);
						
						doublevec_add(myForce.f, df);
						doublevec_kreuzprod(dt, (*mlc1).dr, df);
						doublevec_add(myForce.t, dt);
					}
				}
			}
		}
			
			
		// now -- for walls -- check, if the global vdWColor is >= 0 (should remain unset (-1) for molecules built from vdw spheres)
		// Again, periodicity does not make much sense parallel to an infinitely extending wall with an interaction directly off it
		if((clr2 = vs->getWallColor()) >= 0)
		{
			position *wallPos;
			doublevec *wallNormal;
			position abstand;
			
			wallPos = vs->getWallPosition();
			wallNormal = vs->getWallNormal();
			
			int i = 1;
			for(vector<vdWPos>::iterator mlc1 = actual_spheres.begin(); mlc1 != actual_spheres.end(); ++mlc1, ++i)
			{
				clr1 = (*mlc1).color;
				
				doublevec_set(abstand, (*mlc1).pos);
				doublevec_sub(abstand, (*wallPos));
				
				// empty space between the mol and the wall, measured normal to the wall (= shortest distance)
				dist = doublevec_scalprod(abstand, (*wallNormal)) - (*mlc1).r;
				
				// check whether dist > cut-off  =>  skip further caluclation
				if((dist <= paramsArray[clr1][clr2].Rmax))
				{
					// check for minimal distance
					if(dist < paramsArray[clr1][clr2].Rmin) {
						dist = paramsArray[clr1][clr2].Rmin;
					}
					
					// --- here comes the interaction itself
					dist += paramsArray[clr1][clr2].DR;
					R0 = paramsArray[clr1][clr2].R_0;
					
					rr1 = R0 / dist;		// powers of the inverse distance
					rr3 = rr1 * rr1 * rr1;
					rr4 = rr3 * rr1;
					rr6 = rr3 * rr3;
					// rr7 = rr6 * rr1;
					
					// prefactor = (12.0 * paramsArray[clr1][clr2].C_12 * rr6 + 6.0 * paramsArray[clr1][clr2].C_6) * rr7 / R0;
					prefactor = ((12.0*paramsArray[clr1][clr2].C_12 * rr6 + 6.0*paramsArray[clr1][clr2].C_6) * rr3 + 3.0*paramsArray[clr1][clr2].C_3) * rr4 / R0;

					// calculate the force normal to the wall (normal is normalized -- should be :-))
					doublevec_set(df, (*wallNormal));
					doublevec_mult(df, prefactor);
					
					doublevec_add(myForce.f, df);
					doublevec_kreuzprod(dt, (*mlc1).dr, df);
					doublevec_add(myForce.t, dt);
				}
			}
		}
	}
	
	
	int vdWShape::size() {
		return(spheres.size());
	}
	
	
	void vdWShape::addMyCM(position &theCM, double &theMass)
	{
		doublevec tmppos;
		double m;
		
		for(vector<vdWPos>::iterator it = actual_spheres.begin(); it != actual_spheres.end(); it++)
		{
			m = (*it).r;
			m *= (m*m);		// mass propto r^3

			doublevec_set(tmppos, (*it).pos);
			// std::cerr << "vdWShape::addMyCM: pos = " << tmppos[0]<<", "<<tmppos[1]<<", "<<tmppos[2]<< std::endl;
			
			doublevec_mult(tmppos, m);
			doublevec_add(theCM, tmppos);
			theMass += m;
			
			// std::cerr << "vdWShape::addMyCM: r = " << (*it).r << "  m = " << m << "   tmppos = " << tmppos[0]<<", "<<tmppos[1]<<", "<<tmppos[2]<< std::endl;
		}
	}
	
	
	void vdWShape::dumpStructure() {
		std::cout << "\n\t\tVdwShape mySpheres" << std::endl;
		
		for(std::vector<vdWPos>::iterator it = spheres.begin(); it != spheres.end(); it++) {
			std::cout << "\t\t\t"<< (*it).pos[0]<<" "<<(*it).pos[1]<<" "<<(*it).pos[2]<<"   "<< (*it).r << "   " << (*it).color << std::endl;
		}
		
		std::cout << "\t\tEnd" << std::endl;
	}

	
	void vdWShape::dumpSpheres(char *prefix) 
	{
		int i = 1;
		for(std::vector<vdWPos>::iterator it = actual_spheres.begin(); it != actual_spheres.end(); it++, i++) {
			std::cout << prefix <<" " << i << " "<< (*it).pos[0]<<" "<<(*it).pos[1]<<" "<<(*it).pos[2]<<" "<< (*it).r << " " << (*it).color << std::endl;
		}

	}

	
	void vdWShape::dumpVMD(char *prefix) 
	{
		int i = 1;
		for(std::vector<vdWPos>::iterator it = actual_spheres.begin(); it != actual_spheres.end(); it++, i++) {
			// std::cout << prefix <<" " << i << " "<< (*it).pos[0]<<" "<<(*it).pos[1]<<" "<<(*it).pos[2]<<" "<< (*it).r << " " << (*it).color << std::endl;
			std::cout << "graphics 1 sphere {" << 10.0*(*it).pos[0] <<" "<<10.0*(*it).pos[1]<<" "<<10.0*(*it).pos[2]<<"} radius " << (*it).r*10.0 << " resolution 28" << std::endl;
		}

	}

}