#ifndef _INCLUDED_L3_SPAR_H
#define _INCLUDED_L3_SPAR_H

// *Smoothed particles* Copyright (C) Krzysztof Bosak, 2000-03-27...2000-06-26
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include "l3_array.h"
#include "l3_vec.h"
#include "l3_queue.h"

class Force
{
public:
	virtual void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& r,
		const basearray<Vector3D>& v,
		const basearray<double>& m_recip,
		double dt
		) =0;
	virtual void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const =0;
	virtual ~Force()
	{
	}
};

class ImposedGravity: public Force
{
private:
	double _g;
	Vector3D _gvec;

public:
	inline ImposedGravity()
		: _g(9.81),
		_gvec(0, -9.81, 0)
	{
	}
	inline explicit ImposedGravity(double g)
		: _g(g),
		_gvec(0, -g, 0)
	{
		assert(g>0);
	}
	virtual inline void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& ,
		const basearray<Vector3D>& ,
		const basearray<double>& ,
		double
		)
	{
		const int size=a.size();
		for(int i=0; i<size; i++)
		{
			a[i].y-=_g;
		}
	}
	virtual inline void ShowForces(growing_arrayqueue<Line3D>& , const basearray<Vector3D>& ) const
	//inline void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const
	{
		/*
		for(int i=0; i<r.size(); i++)
		{
			forces.push(Line3D(r[i], r[i]+_gvec));
		}*/
	}
};

class ElasticEarth: public Force
{
private:
	double _k;// elasticity of the ground - Captain Hook's constant
	double _ground_level;

public:
	inline ElasticEarth(double k, double zlevel)
		: _k(k),
		_ground_level(zlevel)
	{
		assert(k>0);
	}
	virtual inline void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& r,
		const basearray<Vector3D>& ,
		const basearray<double>& m_recip,
		double
		)
	{
		const int size=a.size();
		for(int i=0; i<size; i++)
		{
			const double deep=_ground_level-r[i].y;
			if(deep>0)
			{
				a[i].y+=_k*deep*m_recip[i];
			}
		}
	}
	virtual inline void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const
	{
		for(int i=0; i<r.size(); i++)
		{
			const double deep=_ground_level-r[i].y;
			if(deep>0)
			{
				forces.push(Line3D(r[i], r[i]+Vector3D(0, deep, 0)));
			}
		}
	}
	inline double GroundLevel() const
	{
		return _ground_level;
	}
	inline double GroundHookConstant() const
	{
		return _ground_level;
	}
};

class ViscuousElasticEarth: public ElasticEarth
{
private:
	double _viscosity;// viscosity while touching the ground

public:
	inline ViscuousElasticEarth(double k, double zlevel, double viscosity)
		: ElasticEarth(k, zlevel),
		_viscosity(viscosity)
	{
		assert(viscosity>0);
	}
	virtual inline void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& r,
		const basearray<Vector3D>& v,
		const basearray<double>& m_recip,
		double dt
		)
	{
		for(int i=0; i<a.size(); i++)
		{
			const double deep=ElasticEarth::GroundLevel()-r[i].y;
			if(deep>0)
			{
				const Vector3D da_viscuous=_viscosity*v[i]*m_recip[i];
				const double da_deep=ElasticEarth::GroundHookConstant()*deep*m_recip[i];
				if(da_viscuous.Length()*dt >= (v[i]+Vector3D(0, da_deep*dt, 0)).Length())
				{
					a[i]-=v[i]/dt;// vicosity is strong, velocity contribution is ommited
				}
				else
				{
					a[i].y+=da_deep;
					a[i]-=da_viscuous;
				}
			}
		}
	}
	virtual inline void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const
	{
		for(int i=0; i<r.size(); i++)
		{
			const double deep=ElasticEarth::GroundLevel()-r[i].y;
			if(deep>0)
			{
				forces.push(Line3D(r[i], r[i]+Vector3D(0, deep, 0)));
			}
		}
	}
};

class ElasticRepulsion: public Force
{
private:
	double _hook_repulsive;
	double _stable_radius;
	double _viscosity;
	int _particle1_index;
	int _particle2_index;

	double _temp_stable_radius_pow2;

public:
	inline ElasticRepulsion(
		double hook_repulsive,
		double stable_radius,
		double viscosity,
		int particle1_index,
		int particle2_index
		)
		:_hook_repulsive(hook_repulsive),
		_stable_radius(stable_radius),
		_viscosity(viscosity),
		_particle1_index(particle1_index),
		_particle2_index(particle2_index),
		_temp_stable_radius_pow2(stable_radius*stable_radius)
	{
		assert(hook_repulsive>0);
		assert(stable_radius>0);
		assert(viscosity>=0);
		assert(particle1_index!=particle2_index);
		assert(particle1_index>=0);
		assert(particle2_index>=0);
	}
	virtual void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& r,
		const basearray<Vector3D>& v,
		const basearray<double>& m_recip,
		double
		)
	{
		const Vector3D distance_vector=r[_particle1_index]-r[_particle2_index];
		const double distance_vector_length=distance_vector.Length();
		if(distance_vector_length<_stable_radius)
		{
			const double hook_len=(_stable_radius-distance_vector_length)*_hook_repulsive;
			const double visc_len1=_viscosity*v[_particle1_index].Length();
			if(hook_len>visc_len1)
			{
				a[_particle1_index]+=(distance_vector*(_hook_repulsive-visc_len1))*m_recip[_particle1_index];
			}
			const double visc_len2=_viscosity*v[_particle2_index].Length();
			if(hook_len>visc_len2)
			{
				a[_particle2_index]-=(distance_vector*(_hook_repulsive-visc_len2))*m_recip[_particle2_index];
			}
		}
	}
	virtual inline void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const
	{
		Vector3D fdist=r[_particle1_index]-r[_particle2_index];
		if(fdist.LengthPow2()<=_temp_stable_radius_pow2*1.1)
		{
			forces.push(Line3D(r[_particle1_index], r[_particle2_index]));
		}
	}
};

class Elasticity: public Force
{
private:
	double _hook_attractive;
	double _hook_repulsive;

	double _stable_radius;
	double _min_elongation;
	double _max_elongation;
	double _min_radius;
	double _max_radius;
	double _temp_hasr;

	double _viscosity;
	int _particle1_index;
	int _particle2_index;
	bool _broken;

public:
	inline Elasticity(
		double hook_attractive,
		double hook_repulsive,
		double stable_radius,
		double min_elongation,
		double max_elongation,
		double viscosity,
		int particle1_index,
		int particle2_index
		)
		:_hook_attractive(hook_attractive),
		_hook_repulsive(hook_repulsive),
		_stable_radius(stable_radius),
		_min_elongation(min_elongation),
		_max_elongation(max_elongation),
		_min_radius(stable_radius*min_elongation),
		_max_radius(stable_radius*max_elongation),
		_temp_hasr(hook_attractive*stable_radius),
		_viscosity(viscosity),
		_particle1_index(particle1_index),
		_particle2_index(particle2_index),
		_broken(false)
	{
		assert(hook_attractive>=0);
		assert(hook_repulsive>=0);
		assert(stable_radius>0);
		assert(min_elongation>=0);
		assert(min_elongation<1);
		assert(max_elongation>1);
		assert(min_elongation<max_elongation);
		assert(viscosity>=0);
		assert(particle1_index!=particle2_index);
		assert(particle1_index>=0);
		assert(particle2_index>=0);
	}
	virtual void Accelerate(
		basearray<Vector3D>& a,
		const basearray<Vector3D>& r,
		const basearray<Vector3D>& v,
		const basearray<double>& m_recip,
		double
		)
	{
		Vector3D fdist=(r[_particle1_index]-r[_particle2_index]);
		const double fl=fdist.Length();
		if(fl>_stable_radius)
		{
			if(fl>_max_radius)
			{
				_broken=true;
			}
			else if(_broken==false)
			{
				fdist*=_hook_attractive-_temp_hasr/fl;
				a[_particle1_index]+=(-fdist-_viscosity*v[_particle1_index])*m_recip[_particle1_index];
				a[_particle2_index]+=(fdist-_viscosity*v[_particle2_index])*m_recip[_particle2_index];
			}
		}
		else
		{
			fdist*=_hook_repulsive;
			a[_particle1_index]+=(fdist-_viscosity*v[_particle1_index])*m_recip[_particle1_index];
			a[_particle2_index]+=(-fdist-_viscosity*v[_particle2_index])*m_recip[_particle2_index];
		}
	}
	virtual inline void ShowForces(growing_arrayqueue<Line3D>& forces, const basearray<Vector3D>& r) const
	{
		if(_broken==false)
		{
			forces.push(Line3D(r[_particle1_index], r[_particle2_index]));
		}
		else
		{
			const Vector3D distance_vector=(r[_particle1_index]-r[_particle2_index]);
			if(distance_vector.Length()<=_stable_radius*1.1)
			{
				forces.push(Line3D(r[_particle1_index], r[_particle2_index]));
			}
		}
	}
};

class SmoothedParticles
{
private:
	int _total;
	double _total_mass;
	double _dt, _dt2, _1_dt, _1_2dt;
	int _numforces;
	bool _to_be_reheated;
	int _max_untouchables;
	basearray<int> _untouchable_particle;
	basearray<double> _m;
	basearray<double> _m_recip;
	basearray<Vector3D> _r;
	basearray<Vector3D> _r1;
	basearray<Vector3D> _rnew;
	basearray<Vector3D> _v;
	basearray<Vector3D> _vnew;
	basearray<Vector3D> _a;
	basearray<Vector3D> _anew;
	basearray<Force*> _force;// xxx it is up to used to allocate and deallocate forces outside of container.
	inline void UpdateForces(basearray<Vector3D>& a, basearray<Vector3D>& r, basearray<Vector3D>& v)
	{
		a.erase();
		int i;
		for(i=0; i<_numforces; i++)
		{
			_force[i]->Accelerate(a, r, v, _m_recip, _dt);
		}

		for(i=0; i<_max_untouchables; i++)
		{
			a[_untouchable_particle[i]].Clear();
		}
	}

public:
	inline SmoothedParticles()
		:_total(0),
		_total_mass(0),
		_numforces(0),
		_to_be_reheated(true),
		_max_untouchables(0)
	{
		SetStepSize(1);
	}
	inline void AddParticle(
		double mass,
		const Vector3D& position,
		const Vector3D& velocity,
		bool untouchable)
	{
		assert(mass>0);
		const int tp1=_total+1;
		if(_m.size()<tp1)
		{
			const int setsize=tp1<<1;
			_untouchable_particle.resize(setsize);
			_m.resize(setsize);
			_m_recip.resize(setsize);
			_r.resize(setsize);
			_r1.resize(setsize);
			_rnew.resize(setsize);
			_v.resize(setsize);
			_vnew.resize(setsize);
			_a.resize(setsize);
			_anew.resize(setsize);
		}
		if(untouchable==true)
		{
			_untouchable_particle[_max_untouchables]=_total;
			_max_untouchables++;
		}
		_m[_total]=mass;
		_m_recip[_total]=1./mass;
		_total_mass+=mass;
		_r[_total]=position;
		_v[_total]=velocity;
		_total++;
	}
	inline void Reheat()
	{
		assert(_to_be_reheated==true);
		_to_be_reheated=false;
		_untouchable_particle.setsize(_total);
		_m.setsize(_total);
		_m_recip.setsize(_total);
		_r.setsize(_total);
		_r1.setsize(_total);
		_rnew.setsize(_total);
		_v.setsize(_total);
		_vnew.setsize(_total);
		_a.setsize(_total);
		_anew.setsize(_total);
		_r1=_r;
		for(int j=0; j<_total; j++)
		{
			for(int i=j+1; i<_total; i++)
			{// Any system with particles placed too close would be extremely unstable.
				assert( Distance(_r[i], _r[j]) > 1e-14 );
			}
		}
		_dt/=10;
		for(int s=0; s<10; s++)
			DoEulerTriplePCStep();
		_dt*=10;
	}
	inline void SetStepSize(double stepsize)
	{
		assert(stepsize>0);
		_dt=stepsize;
		_dt2=stepsize*stepsize;
		_1_2dt=0.5/stepsize;
		_1_dt=1/stepsize;
	}
	inline void AddForce(Force * const forceobject)
	{
		if(_force.size()<_numforces+1)
		{
			_force.resize((_numforces+1)<<1);
		}
		_force[_numforces]=forceobject;
		_numforces++;
	}
	inline void DoEulerStep()
	{
		UpdateForces(_a, _r, _v);
		int i;
		for(i=0; i<_total; i++)
		{
			_r[i]+=_v[i]*_dt;
			_v[i]+=_a[i]*_dt;
		}
	}
	inline void DoEulerPCStep()
	{
		UpdateForces(_a, _r, _v);
		int i;
		for(i=0; i<_total; i++)
		{
			_r[i]+=_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _r, _vnew);
		for(i=0; i<_total; i++)
		{
			_v[i]+=_a[i]*_dt;
		}
	}
	inline void DoEulerDoublePCStep()
	{
		int i;
		UpdateForces(_a, _r, _v);
		for(i=0; i<_total; i++)
		{
			_rnew[i]=_r[i]+_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _rnew, _vnew);
		for(i=0; i<_total; i++)
		{
			_r[i]+=_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _r, _vnew);
		for(i=0; i<_total; i++)
		{
			_v[i]+=_a[i]*_dt;
		}
	}
	inline void DoEulerTriplePCStep()
	{
		int i;
		UpdateForces(_a, _r, _v);
		for(i=0; i<_total; i++)
		{
			_rnew[i]=_r[i]+_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _rnew, _vnew);
		for(i=0; i<_total; i++)
		{
			_rnew[i]=_r[i]+_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _rnew, _vnew);
		for(i=0; i<_total; i++)
		{
			_r[i]+=_v[i]*_dt;
			_vnew[i]=_v[i]+_a[i]*_dt;
		}
		UpdateForces(_a, _r, _vnew);
		for(i=0; i<_total; i++)
		{
			_v[i]+=_a[i]*_dt;
		}
	}
	inline void DoVerletStep()
	{
		assert(_to_be_reheated==false);
		UpdateForces(_a, _r, _v);
		for(int i=0; i<_total; i++)
		{
			_rnew[i]=_dt2*_a[i]+_r[i]+_r[i]-_r1[i];
			_v[i]=(_rnew[i]-_r1[i])*_1_2dt;
		}
		_r.swap(_r1);
		_rnew.swap(_r);
	}
	inline void DoMidpointStep()
	{	// also known as MidEuler
		int i;
		UpdateForces(_a, _r, _v);
		for(i=0; i<_total; i++)
		{
			_rnew[i]=_r[i]+_v[i]*_dt*.5;
			_vnew[i]=_v[i]+_a[i]*_dt*.5;
		}
		UpdateForces(_a, _rnew, _vnew);
		for(i=0; i<_total; i++)
		{
			_r[i]+=_v[i]*_dt;
			_v[i]+=_a[i]*_dt;
		}
	}
	/*
	inline void DoHeunStep()
	{
		static basearray<Vector3D> k1v(_total);
		static basearray<Vector3D> k1a(_total);
		int i;
		UpdateForces(_a, _r, _v);
		for(i=0; i<_total; i++)
		{
			k1v[i]=_v[i]*_dt;
			k1a[i]=_a[i]*_dt;
			_rnew[i]=_r[i]+k1v[i]*2/3;
			_vnew[i]=_v[i]+k1a[i]*2/3;
		}
		UpdateForces(_a, _rnew, _vnew);
		for(i=0; i<_total; i++)
		{
			_r[i]+=_vnew[i]*_dt*.25+k1v[i]*.75;
			_v[i]+=_anew[i]*_dt*.25+k1a[i]*.75;
		}
	}
	*/
	inline Vector3D CenterOfMass_Position() const
	{// Position of center of mass
		Vector3D center(0, 0, 0);
		for(int i=0; i<_total; i++)
		{
			center+=_r[i]*_m[i];
		}
		return center/_total_mass;
	}
	inline Vector3D CenterOfMass_Velocity() const
	{// Velocity of center of mass
		Vector3D center(0, 0, 0);
		for(int i=0; i<_total; i++)
		{
			center+=_v[i]*_m[i];
		}
		return center/_total_mass;
	}
	inline Vector3D Position(int i) const
	{// Single real object position, passed by value: position array may grow and deplace during system configuration phase
		assert(i>=0);
		assert(i<_total);
		return _r[i];
	}
	inline Vector3D Velocity(int i) const
	{// Single real object velocity, passed by value: velocity array may grow and deplace during system configuration phase
		assert(i>=0);
		assert(i<_total);
		return _v[i];
	}
	inline double Mass(int i) const
	{// Single real object mass, passed by value: mass array may grow and deplace during system configuration phase
		assert(i>=0);
		assert(i<_total);
		return _m[i];
	}
	inline const basearray<Vector3D> & Positions() const
	{// Real objects positions, can no be used suring assemblying (due to reference and growing array of particles)
		assert(_to_be_reheated==false);
		assert(_r.size()==_total);
		return _r;
	}
	inline void ShowForces(growing_arrayqueue<Line3D>& forces) const
	{// ShowForces imposed on real objects
		forces.clear();
		for(int i=0; i<_numforces; i++)
		{
			_force[i]->ShowForces(forces, _r);
		}
	}
	inline int Total() const
	{
		return _total;
	}
	inline double KineticEnergy() const
	{
		double kenergy=0;
		for(int i=0; i<_total; i++)
		{
			kenergy+=_m[i]*_v[i].LengthPow2();
		}
		return kenergy*.5;
	}
	inline double KineticEnergy_PerParticle() const
	{
		return KineticEnergy()/_total;
	}
};

class SmoothedLine
{
private:
	basearray<int> _nodes;

public:
	inline SmoothedLine(SmoothedParticles& sp, int p1, int p2, int sections,
		double k1, double k2, double elongation, double viscosity)
		:_nodes(sections+1)
	{
		assert(sections>=2);
		assert(p1!=p2);
		const double m_start=sp.Mass(p1);
		const double m_end=sp.Mass(p2);
		const double m_grad=m_end-m_start;
		const Vector3D r_start=sp.Position(p1);
		const Vector3D r_end=sp.Position(p2);
		const Vector3D r_vector=r_end-r_start;
		const Vector3D v_start=sp.Velocity(p1);
		const Vector3D v_end=sp.Velocity(p2);
		const Vector3D v_vector=v_end-v_start;

		int i;

		_nodes[0]=p1;
		int curnod=sp.Total();
		for(i=1; i<sections; i++)
		{
			_nodes[i]=curnod;
			curnod++;
		}
		_nodes[i]=p2;

		for(i=1; i<sections; i++)
		{
			const double factor=static_cast<double>(i)/sections;
			const double mass=m_start+m_grad*factor;
			const Vector3D position=r_start+r_vector*factor;
			const Vector3D velocity=v_start+v_vector*factor;
			sp.AddParticle(mass, position, velocity, false);
		}
		const double r_stable=r_vector.Length()/sections;
		for(i=1; i<=sections; i++)
		{
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, _nodes[i-1], _nodes[i]));
		}
	}
	inline const basearray<int> & Nodes() const
	{
		return _nodes;
	}
};

class SmoothedHatch
{
private:
	const SmoothedLine& _smoothed_line_1;
	const SmoothedLine& _smoothed_line_2;

public:
	SmoothedHatch(SmoothedParticles& sp, const SmoothedLine& smoothed_line_1, const SmoothedLine& smoothed_line_2,
		double k1, double k2, double elongation, double viscosity)
		:_smoothed_line_1(smoothed_line_1),
		_smoothed_line_2(smoothed_line_2)
	{
		assert(smoothed_line_1.Nodes().size()==smoothed_line_2.Nodes().size());
		assert(&smoothed_line_1!=&smoothed_line_2);
		const basearray<int> p1=smoothed_line_1.Nodes();
		const basearray<int> p2=smoothed_line_2.Nodes();
		for(int i=0; i<p1.size(); i++)
		{
			const double r_stable=( sp.Position(p1[i])-sp.Position(p2[i]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i], p2[i]));
		}
	}

private:
	SmoothedHatch& operator=(const SmoothedHatch&);// Forbidden.
	explicit SmoothedHatch(const SmoothedHatch&);// Forbidden.
};

class SmoothedXHatch
{
private:
	const SmoothedLine& _smoothed_line_1;
	const SmoothedLine& _smoothed_line_2;

public:
	SmoothedXHatch(SmoothedParticles& sp, const SmoothedLine& smoothed_line_1, const SmoothedLine& smoothed_line_2,
		double k1, double k2, double elongation, double viscosity)
		:_smoothed_line_1(smoothed_line_1),
		_smoothed_line_2(smoothed_line_2)
	{
		assert(smoothed_line_1.Nodes().size()==smoothed_line_2.Nodes().size());
		assert(&smoothed_line_1!=&smoothed_line_2);
		const basearray<int> p1=smoothed_line_1.Nodes();
		const basearray<int> p2=smoothed_line_2.Nodes();
		for(int i=1; i<p1.size(); i++)
		{
			double r_stable=( sp.Position(p1[i-1])-sp.Position(p2[i]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i-1], p2[i]));
			r_stable=( sp.Position(p1[i])-sp.Position(p2[i-1]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i], p2[i-1]));
		}
	}

private:
	SmoothedXHatch& operator=(const SmoothedXHatch&);// Forbidden.
	explicit SmoothedXHatch(const SmoothedXHatch&);// Forbidden.
};

class SmoothedXBoxHatch
{
private:
	const SmoothedLine& _smoothed_line_1;
	const SmoothedLine& _smoothed_line_2;

public:
	SmoothedXBoxHatch(SmoothedParticles& sp, const SmoothedLine& smoothed_line_1, const SmoothedLine& smoothed_line_2,
		double k1, double k2, double elongation, double viscosity)
		:_smoothed_line_1(smoothed_line_1),
		_smoothed_line_2(smoothed_line_2)
	{
		assert(smoothed_line_1.Nodes().size()==smoothed_line_2.Nodes().size());
		assert(&smoothed_line_1!=&smoothed_line_2);
		const basearray<int> p1=smoothed_line_1.Nodes();
		const basearray<int> p2=smoothed_line_2.Nodes();
		double r_stable=( sp.Position(p1[0])-sp.Position(p2[0]) ).Length();
		sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[0], p2[0]));
		for(int i=1; i<p1.size(); i++)
		{
			r_stable=( sp.Position(p1[i-1])-sp.Position(p2[i]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i-1], p2[i]));
			r_stable=( sp.Position(p1[i])-sp.Position(p2[i-1]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i], p2[i-1]));
			r_stable=( sp.Position(p1[i])-sp.Position(p2[i]) ).Length();
			sp.AddForce(new Elasticity(k1, k2, r_stable, 0, elongation, viscosity, p1[i], p2[i]));
		}
	}

private:
	SmoothedXBoxHatch& operator=(const SmoothedXBoxHatch&);// Forbidden.
	explicit SmoothedXBoxHatch(const SmoothedXBoxHatch&);// Forbidden.
};

#endif //_INCLUDED_L3_SPAR_H