
Thanks for the report.
using SharpDX;
namespace SharpNVE
{
public sealed class Gear5<T> : IPredictorCorrector<T> where T : ISimulationParticle
{
private Vector3[] ThirdDerivatives { get; set; }
private Vector3[] FourthDerivatives { get; set; }
public Gear5(int particleCount)
{
ThirdDerivatives = new Vector3[particleCount];
FourthDerivatives = new Vector3[particleCount];
}
public override void Predict(float timeStep, SimulationBox<T> simulationBox)
{
var c1 = timeStep;
var c2 = c1 * timeStep / 2.0f;
var c3 = c2 * timeStep / 3.0f;
var c4 = c3 * timeStep / 4.0f;
for (var i = 0; i < simulationBox.Particles.Length; ++i)
{
var thirdDerivative = ThirdDerivatives[i];
var fourthDerivative = FourthDerivatives[i];
simulationBox.Particles[i].Position += c1 * simulationBox.Particles[i].Velocity +
c2 * simulationBox.Particles[i].Acceleration +
c3 * thirdDerivative + c4 * fourthDerivative;
simulationBox.Particles[i].Velocity += c1 * simulationBox.Particles[i].Acceleration +
c2 * thirdDerivative + c3 * fourthDerivative;
simulationBox.Particles[i].Acceleration += c1 * thirdDerivative + c2 * fourthDerivative;
ThirdDerivatives[i] += c1 * fourthDerivative;
}
}
public override void Correct(float timeStep, SimulationBox<T> simulationBox,
float atomicMass, out float kineticEnergyPerAtom)
{
var c1 = timeStep;
var c2 = c1 * timeStep / 2.0f;
var c3 = c2 * timeStep / 3.0f;
var c4 = c3 * timeStep / 4.0f;
const float gear0 = 19.0f / 120.0f;
const float gear1 = 3.0f / 4.0f;
const float gear3 = 1.0f / 2.0f;
const float gear4 = 1.0f / 12.0f;
var cr = gear0 * c2;
var cv = gear1 * c2 / c1;
var cb = gear3 * c2 / c3;
var cc = gear4 * c2 / c4;
kineticEnergyPerAtom = 0.0f;
for (var i = 0; i < simulationBox.Particles.Length; ++i)
{
var ax = simulationBox.Particles[i].Forces / atomicMass;
var corr = ax - simulationBox.Particles[i].Acceleration;
simulationBox.Particles[i].Position += cr * corr;
simulationBox.Particles[i].Velocity += cv * corr;
simulationBox.Particles[i].Acceleration = ax;
ThirdDerivatives[i] += cb * corr;
FourthDerivatives[i] += cc * corr;
kineticEnergyPerAtom += simulationBox.Particles[i].Velocity.LengthSquared();
}
kineticEnergyPerAtom *= 0.5f * atomicMass;
}
}
}
using SharpDX;
namespace SharpNVE
{
public sealed class LennardJones<T> : IInteractionCalculator<T> where T : ISimulationParticle
{
public override void ComputeForces(SimulationBox<T> simulationBox)
{
Potential = 0.0f;
Virial = 0.0f;
foreach (var simulationParticle in simulationBox.Particles)
simulationParticle.Forces = Vector3.Zero;
for (var i = 0; i < simulationBox.Particles.Length - 1; ++i)
{
var iPosition = simulationBox.Particles[i].Position;
for (var j = 1; j < simulationBox.Particles.Length; ++j)
{
var jPosition = simulationBox.Particles[j].Position;
var displacementVector = jPosition - iPosition;
var distance = Vector3.DistanceSquared(iPosition, jPosition);
if (distance >= SquaredCutoffRadius)
continue;
var invDistSquare = 1.0f / distance;
var invDistSixtic = invDistSquare * invDistSquare * invDistSquare;
var invDistTwelve = invDistSixtic * invDistSixtic;
var vDiff = invDistTwelve - invDistSixtic;
Potential += vDiff;
var wDiff = vDiff + invDistTwelve;
Virial += wDiff;
var forcesDiff = wDiff * invDistSquare;
displacementVector *= forcesDiff;
simulationBox.Particles[i].Forces += displacementVector;
simulationBox.Particles[j].Forces -= displacementVector;
}
}
foreach (var simulationParticle in simulationBox.Particles)
simulationParticle.Forces *= 24.0f;
Potential *= 4.0f;
Virial *= 24.0f / 3.0f;
}
}
}
F5BQP is trying to use copper-loaded filament to #3Dprint a 24 GHz horn antenna. Measurements expected. pic.twitter.com/aO3kPQ0QXf
— S☰bastien F4GRX (@f4grx) March 21, 2017