Source

NBody / NBody.DomainModel / Integrators / HermiteIntegrator.cs

Full commit
Ade Miller b47cfa5 










Ade Miller f2153ee 
Ade Miller b47cfa5 
Ade Miller f2153ee 


Ade Miller b47cfa5 


Ade Miller f2153ee 






Ade Miller b47cfa5 
Ade Miller f2153ee 

Ade Miller b47cfa5 
Ade Miller f2153ee 


Ade Miller b47cfa5 
Ade Miller f2153ee 





Ade Miller b47cfa5 
Ade Miller f2153ee 

Ade Miller b47cfa5 
Ade Miller f2153ee 


Ade Miller b47cfa5 
Ade Miller f2153ee 
Ade Miller b47cfa5 
Ade Miller f2153ee 







Ade Miller b47cfa5 
Ade Miller f2153ee 







Ade Miller b47cfa5 
Ade Miller f2153ee 


Ade Miller b47cfa5 
Ade Miller f2153ee 
Ade Miller b47cfa5 
Ade Miller f2153ee 

Ade Miller b47cfa5 
Ade Miller f2153ee 





Ade Miller b47cfa5 
// Copyright (c) Ade Miller.  All Rights Reserved.
// This code released under the terms of the
// Microsoft Public License (MS-PL, http://opensource.org/licenses/ms-pl.html.)

#if HERMITE_SUPPORTED
using System;
using System.ComponentModel;
#endif

namespace NBody.DomainModel.Integrators
{
	// Example: http://www.artcompsci.org/vol_1/v1_web/node46.html
#if HERMITE_SUPPORTED
	[Description("Hermite (C#)")]
	public class HermiteIntegrator : IIntegrate
	{
		[ReadOnly(true)]
		public double SofteningLength { get; set; }

		public void Initialize(Body[] bodies)
		{
			for (int i = 0; i < bodies.Length; i++)
			{
				bodies[i].Acceleration = new Vector();
				bodies[i].Jerk = new Vector();
			}

			CalculateAccellerationAndJerk(bodies);
		}

		public void Integrate(double dT, Body[] bodies)
		{
			double dT2 = dT * dT;

			for (int i = 0; i < bodies.Length; i++)
			{
				bodies[i].OldPosition = bodies[i].Position;
				bodies[i].OldVelocity = bodies[i].Velocity;
				bodies[i].OldAcceleration = bodies[i].Acceleration;
				bodies[i].OldJerk = bodies[i].Jerk;

				bodies[i].Position += bodies[i].Velocity * dT + bodies[i].Acceleration * dT2 / 2 + bodies[i].Jerk * dT2 * dT / 6;
				bodies[i].Velocity += bodies[i].Acceleration * dT + bodies[i].Jerk * dT2 / 2;

				bodies[i].Acceleration = new Vector();
				bodies[i].Jerk = new Vector();
			}

			CalculateAccellerationAndJerk(bodies);

			for (int i = 0; i < bodies.Length; i++)
			{
				bodies[i].Velocity = bodies[i].OldVelocity + (bodies[i].OldAcceleration + bodies[i].Acceleration) * dT / 2   // Eq. 6.2
							  + (bodies[i].OldJerk - bodies[i].Jerk) * dT2 / 12;
				bodies[i].Position = bodies[i].OldPosition + (bodies[i].OldVelocity + bodies[i].Velocity) * dT / 2             // Eq. 6.1
							  + (bodies[i].OldAcceleration - bodies[i].Acceleration) * dT2 / 12;
			}
		}

		private void CalculateAccellerationAndJerk(Body[] bodies)
		{
			for (int i = 0; i < bodies.Length; i++)
			{
				for (int j = (i + 1); j < bodies.Length; j++)
				{
					Vector r = bodies[j].Position - bodies[i].Position;
					Vector v = bodies[j].Velocity - bodies[i].Velocity;

					double r2 = r.MagnitudeSquared;
					double r3 = Math.Sqrt(r2) * r2;
					double rv = Vector.ScalarProduct(r, v) / r2;

					Vector f = r * (bodies[i].Mass * bodies[j].Mass / r3);        // Don't use Body.Force because that would mean recalculation of r3 later.

					bodies[i].Acceleration += f / bodies[i].Mass;
					bodies[j].Acceleration -= f / bodies[j].Mass;

					bodies[i].Jerk += bodies[j].Mass * (v - (3 * rv * r)) / r3;   // Eq. 6.4
					bodies[j].Jerk -= bodies[i].Mass * (v - (3 * rv * r)) / r3;
				}
			}
		}
	}
#endif
}