#include "AdunKernel/AdunBerendsenThermostat.h"

@implementation AdBerendsenThermostat

- (id) init
{
	return [self initWithTargetTemperature: 300.0
		couplingFactor: 100.0];
}

- (id) initWithTargetTemperature: (double) doubleOne
	couplingFactor: (double) doubleTwo
{
	if((self = [super init]))
	{
		[self setTargetTemperature: doubleOne];
		couplingFactor = doubleTwo;
		timeStep = 1;
		timePerCouplingFactor = (double)timeStep/couplingFactor;
	}

	return self;
}

- (NSString*) description
{
	NSMutableString* description = [NSMutableString string];
	
	[description appendFormat: @"%@. Target temperature %5.2lf. Coupling factor %5.2lf", 
		NSStringFromClass([self class]), targetTemperature, couplingFactor];
	
	return description;	
}

- (void) simulator: (AdSimulator*) aSimulator 
		willBeginProductionWithSystems: (AdSystemCollection*) aSystemCollection 
		forceFields: (AdForceFieldCollection*) aForceFieldCollection
{
	timeStep = [aSimulator timeStep];
	timePerCouplingFactor = (double)timeStep/couplingFactor;
}

- (void) simulatorDidPerformSecondVelocityUpdateForSystem: (AdSystem*) system
{
	int j,k;
	double temperature, factor;
	AdMatrix* velocities;

	temperature = [system temperature];

	if(temperature == 0)
	{
		factor = 0;
	}	
	else
	{
		factor = 1 - timePerCouplingFactor*(1 - (targetTemperature/temperature));
		factor = sqrt(factor);
	}	
	
	velocities = [system velocities];

	[system object: self willBeginWritingToMatrix: velocities]; 

	for(j=0; j < velocities->no_rows; j++)
		for(k=0; k<3; k++)
			velocities->matrix[j][k] *= factor;
	
	[system object: self didFinishWritingToMatrix: velocities]; 
}

- (void) setTargetTemperature: (double) aDouble
{
	if(aDouble < 0)
		[NSException raise: NSInvalidArgumentException
			format: @"Temperature cannot be less than 0"];

	targetTemperature = aDouble;
}

- (double) targetTemperature
{
	return targetTemperature;
}

- (void) setCouplingFactor: (double) aDouble
{
	couplingFactor = aDouble;
	timePerCouplingFactor = (double)timeStep/couplingFactor;
}

- (double) couplingFactor
{
	return couplingFactor;
}
	
- (void) simulatorWillPerformFirstVelocityUpdateForSystem: (AdSystem*) aSystem
{
	//Does nothing here
}

- (void) simulatorWillPerformPositionUpdateForSystem: (AdSystem*) aSystem
{
	//Does nothing here
}

- (void) simulatorDidPerformPositionUpdateForSystem: (AdSystem*) aSystem
{
	//Does nothing here
}

- (void) simulatorWillPerformSecondVelocityUpdateForSystem: (AdSystem*) aSystem
{
	//Does nothing here
}

- (void) simulatorDidFinishProduction: (AdSimulator*) aSimulator
{
	//Does nothing here
}

@end

