local min = math.min
local max = math.max
local abs = math.abs
local physics = {}
physics.__index = physics
local function sign(x)
if x > 0 then return 1 end
if x < 0 then return -1 end
return 0
end
function physics.new(trainConfig, trackSampler)
return setmetatable({
config = trainConfig,
sampler = trackSampler,
}, physics)
end
function physics:getEquivalentMass()
local c = self.config
return (c.serviceMass or c.mass or 112000.0) * (c.rotatingMassFactor or 1.0)
end
function physics:computeRunningResistance(speed, state)
local c = self.config
local v = abs(speed)
local mass = c.serviceMass or c.mass or 112000.0
local rollingCoefficient = c.rollingResistanceCoefficient or 0.0018
local airDensity = c.airDensity or 1.225
local dragCoefficient = c.dragCoefficient or 0.20
local frontalArea = c.frontalArea or 10.9
local rolling = rollingCoefficient * mass * 9.81
local aerodynamic = 0.5 * airDensity * dragCoefficient * frontalArea * v * v
local davis = rolling + aerodynamic
local startingResistance = 0.0
local startSpeed = c.startingResistanceFadeSpeed or 1.4
if v < startSpeed then
local normalizedSpeed = v / max(startSpeed, 0.001)
local fade = 1.0 - normalizedSpeed
fade = max(0.0, min(1.0, fade))
local requestedForce = abs(state.requestedTractiveForce or 0.0)
startingResistance = min(requestedForce, (c.startingResistanceForce or 12000.0) * fade)
end
local drivetrain = c.drivetrainDragForce or 1200.0
local rollingAndMechanical = rolling
return davis + startingResistance + drivetrain, rollingAndMechanical, aerodynamic, startingResistance, drivetrain
end
function physics:computeCurveResistance(position)
local c = self.config
local radius = self.sampler:getMeanRadiusOverTrain(position)
if abs(radius) < 0.001 then
return 0.0, radius
end
local k = c.curveResistanceK or 650.0
local resistancePerMille = k / abs(radius)
local force = (c.serviceMass or c.mass or 112000.0) * 9.81 * (resistancePerMille / 1000.0)
return force, radius
end
function physics:computeGradeForce(position)
local c = self.config
local grade = self.sampler:getMeanGradeOverTrain(position)
local force = (c.serviceMass or c.mass or 112000.0) * 9.81 * grade
return force, grade
end
function physics:computeTractiveForce(state, dt)
local c = self.config
local v = abs(state.speed or 0.0)
local dir = state.driveDirection or 0
if dir == 0 then
return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
end
local tractionFraction = state.tractionFraction or 0.0
local steps = c.tractionForceSteps or c.tractionSteps or {0.00, 0.20, 0.38, 0.58, 0.78, 1.00}
local forceFraction = steps[(state.tractionActualStep or 0) + 1] or 0.0
local powerLimitedForce = ((c.power or 2160000.0) * tractionFraction) / max(v, c.minPowerSpeed or 2.5)
local constantForceLimit = (c.startingTractiveEffort or 135000.0) * forceFraction
local continuousForceLimit = (c.continuousTractiveEffort or 90000.0) + ((c.power or 2160000.0) / max(v, c.minPowerSpeed or 2.5) - (c.continuousPower or 1500000.0) / max(v, c.minPowerSpeed or 2.5))
local adhesionLimitedForce = (c.adhesionCoefficient or 0.185) * (c.serviceMass or c.mass or 112000.0) * 9.81 * (c.creepPeakFactor or 1.02)
local rawRequestedForce = min(constantForceLimit, powerLimitedForce, max(0.0, continuousForceLimit))
local previousRequestedForce = state.requestedTractiveForce or 0.0
local riseTime = c.tractiveForceRiseTime or 1.35
local releaseTime = c.tractiveForceReleaseTime or 0.70
local tau = rawRequestedForce > previousRequestedForce and max(riseTime, 0.001) or max(releaseTime, 0.001)
local blend = 1.0 - math.exp(-dt / tau)
local requestedForce = previousRequestedForce + (rawRequestedForce - previousRequestedForce) * blend
local tractiveMagnitude = min(requestedForce, adhesionLimitedForce)
local slipRatio = 0.0
if adhesionLimitedForce > 0.001 then
slipRatio = max(0.0, requestedForce - adhesionLimitedForce) / adhesionLimitedForce
end
return tractiveMagnitude * dir, requestedForce, rawRequestedForce, powerLimitedForce, adhesionLimitedForce, constantForceLimit, slipRatio
end
function physics:computeBrakeForce(state)
local c = self.config
local v = abs(state.speed or 0.0)
local refSign = sign(state.speed or 0.0)
local brakeMax = state.emergencyBrake and (c.emergencyBrakeForceMax or 165000.0) or (c.serviceBrakeForceMax or 125000.0)
local baseMagnitude = brakeMax * (state.brakeFraction or 0.0)
local lowSpeedBoost = 1.0
local threshold = 20.0 / 3.6
local standstillBoost = state.emergencyBrake and 1.18 or 1.50
if v < threshold then
local t = v / max(threshold, 0.001)
lowSpeedBoost = standstillBoost + (1.0 - standstillBoost) * t
end
local magnitude = baseMagnitude * lowSpeedBoost
if refSign == 0 then
local gradeForce = state.gradeForce
if gradeForce == nil and self.sampler and state.position ~= nil then
local grade = self.sampler:getMeanGradeOverTrain(state.position or 0.0)
gradeForce = (self.config.serviceMass or self.config.mass or 112000.0) * 9.81 * grade
end
local tractiveForce = state.tractiveForce or 0.0
local tendency = sign(tractiveForce + (gradeForce or 0.0))
if tendency == 0 then
return 0.0, magnitude
end
return -tendency * magnitude, magnitude
end
return -refSign * magnitude, magnitude
end
function physics:getAccelerationFor(state, dt)
local eqMass = self:getEquivalentMass()
local tractiveForce, requestedForce, rawRequestedForce, powerLimitedForce, adhesionLimitedForce, constantForceLimit, slipRatio = self:computeTractiveForce(state, dt)
local runningResistance, rollingResistance, aerodynamicDrag, startingResistance, drivetrainDrag = self:computeRunningResistance(state.speed or 0.0, state)
local curveResistance, radius = self:computeCurveResistance(state.position or 0.0)
local gradeForce, grade = self:computeGradeForce(state.position or 0.0)
local brakeForce, brakeForceMagnitude = self:computeBrakeForce(state)
local speedSign = sign(state.speed or 0.0)
local resistanceMagnitude = runningResistance + curveResistance
local resistanceForce = 0.0
if speedSign ~= 0 then
resistanceForce = -speedSign * resistanceMagnitude
else
local tendency = sign(tractiveForce + gradeForce)
resistanceForce = -tendency * resistanceMagnitude
end
local netForce = tractiveForce + resistanceForce + gradeForce + brakeForce
local acceleration = netForce / eqMass
return {
acceleration = acceleration,
tractiveForce = tractiveForce,
requestedTractiveForce = requestedForce,
rawRequestedTractiveForce = rawRequestedForce,
powerLimitedForce = powerLimitedForce,
adhesionLimitedForce = adhesionLimitedForce,
constantForceLimit = constantForceLimit,
slipRatio = slipRatio,
runningResistance = runningResistance,
rollingResistance = rollingResistance,
aerodynamicDrag = aerodynamicDrag,
startingResistance = startingResistance,
drivetrainDrag = drivetrainDrag,
curveResistance = curveResistance,
gradeForce = gradeForce,
brakeForce = brakeForce,
brakeForceMagnitude = brakeForceMagnitude,
resistanceForce = resistanceForce,
resistanceMagnitude = resistanceMagnitude,
netForce = netForce,
meanGrade = grade,
meanRadius = radius,
}
end
function physics:evaluateDerivative(sampleState, dt)
local forces = self:getAccelerationFor(sampleState, dt)
return {
dx = sampleState.speed or 0.0,
dv = forces.acceleration,
forces = forces,
}
end
function physics:step(state, dt)
local function makeSample(base, dx, dv, scale)
return {
position = (base.position or 0.0) + dx * scale,
speed = (base.speed or 0.0) + dv * scale,
driveDirection = base.driveDirection,
tractionActualStep = base.tractionActualStep,
tractionFraction = base.tractionFraction,
brakeFraction = base.brakeFraction,
emergencyBrake = base.emergencyBrake,
requestedTractiveForce = base.requestedTractiveForce,
}
end
local s0 = {
position = state.position or 0.0,
speed = state.speed or 0.0,
driveDirection = state.driveDirection,
tractionActualStep = state.tractionActualStep,
tractionFraction = state.tractionFraction,
brakeFraction = state.brakeFraction,
emergencyBrake = state.emergencyBrake,
requestedTractiveForce = state.requestedTractiveForce or 0.0,
}
local k1 = self:evaluateDerivative(s0, dt)
local k2 = self:evaluateDerivative(makeSample(s0, k1.dx, k1.dv, dt * 0.5), dt)
local k3 = self:evaluateDerivative(makeSample(s0, k2.dx, k2.dv, dt * 0.5), dt)
local k4 = self:evaluateDerivative(makeSample(s0, k3.dx, k3.dv, dt), dt)
state.position = s0.position + (dt / 6.0) * (k1.dx + 2.0 * k2.dx + 2.0 * k3.dx + k4.dx)
state.speed = s0.speed + (dt / 6.0) * (k1.dv + 2.0 * k2.dv + 2.0 * k3.dv + k4.dv)
local finalForces = self:getAccelerationFor(state, dt)
local snapSpeed = self.config.neutralSnapSpeed or 0.04
local snapAccel = self.config.neutralSnapAcceleration or 0.05
local shouldZeroKinematics = false
local holdingBrakeActive = (state.brakeFraction or 0.0) > 0.01 or state.emergencyBrake
local nearStandstill = abs(state.speed) < snapSpeed
local brakeCanHoldTrain = abs((finalForces.tractiveForce or 0.0) + (finalForces.gradeForce or 0.0)) <= (finalForces.brakeForceMagnitude or 0.0)
if nearStandstill and holdingBrakeActive and brakeCanHoldTrain then
state.speed = 0.0
shouldZeroKinematics = true
elseif abs(state.speed) < snapSpeed and abs(finalForces.acceleration) < snapAccel then
state.speed = 0.0
shouldZeroKinematics = true
end
if shouldZeroKinematics then
finalForces.acceleration = 0.0
finalForces.netForce = 0.0
finalForces.brakeForce = 0.0
finalForces.resistanceForce = 0.0
end
state.acceleration = finalForces.acceleration
state.tractiveForce = finalForces.tractiveForce
state.rawRequestedTractiveForce = finalForces.rawRequestedTractiveForce
state.requestedTractiveForce = finalForces.requestedTractiveForce
state.powerLimitedForce = finalForces.powerLimitedForce
state.adhesionLimitedForce = finalForces.adhesionLimitedForce
state.constantForceLimit = finalForces.constantForceLimit
state.slipRatio = finalForces.slipRatio or 0.0
state.slipState = min(1.0, state.slipRatio / (self.config.fullSlipRatio or 0.04))
state.resistanceForce = finalForces.resistanceForce
state.resistanceForceMagnitude = finalForces.resistanceMagnitude
state.rollingResistanceForce = finalForces.rollingResistance
state.aerodynamicDragForce = finalForces.aerodynamicDrag
state.startingResistanceForce = finalForces.startingResistance
state.drivetrainDragForce = finalForces.drivetrainDrag
state.curveResistanceForce = finalForces.curveResistance
state.gradeForce = finalForces.gradeForce
state.brakeForce = finalForces.brakeForce
state.brakeForceMagnitude = finalForces.brakeForceMagnitude
state.netForce = finalForces.netForce
state.meanGrade = finalForces.meanGrade
state.currentGrade = finalForces.meanGrade
state.curveRadius = finalForces.meanRadius
state.visualPitch = -finalForces.meanGrade
state.tractiveForceMagnitude = abs(finalForces.tractiveForce)
end
return physics