for i = 1,N_SCsim do
		-- Set up iteration wheel-rotation points:
		local PartialWheelRot = prevWheelRot + (i/N_SCsim)*(wheelRot - prevWheelRot)

		if PartialWheelRot > 2*math.pi then
			PartialWheelRot = PartialWheelRot - 2*math.pi
		elseif PartialWheelRot < 0 then
			PartialWheelRot = PartialWheelRot + 2*math.pi
		end

		-- Subdivide wheelRotspeed
		prevPartialWheelRotspeed = PartialWheelRotspeed
		PartialWheelRotspeed = prevWheelRotspeed + (i/N_SCsim)*(wheelRotspeed - prevWheelRotspeed)

		-- Cranks
		Cranks[1] = PartialWheelRot + 1.0*math.pi
		Cranks[2] = PartialWheelRot + 0.0*math.pi
		Cranks[3] = PartialWheelRot + 1.5*math.pi
		Cranks[4] = PartialWheelRot + 0.5*math.pi

		for j = 1,4 do
			-- Reverse indices if reverser < 0:
			if revdir < 0 then
				Cranks[j] = Cranks[j] + math.pi
			end

			-- Check out-of-bounds:
			if Cranks[j] > 2*math.pi then
				Cranks[j] = Cranks[j] - 2*math.pi
			end
		end

		-- Obtain Piston Postitions
		PistPos[1] = WheelRot2PistPos(PartialWheelRot)						-- Left, Front.
		PistPos[2] = WheelRot2PistPos(PartialWheelRot + math.pi)			-- Left, Rear.
		PistPos[3] = WheelRot2PistPos(PartialWheelRot + 0.5*math.pi)		-- Right, Front.
		PistPos[4] = WheelRot2PistPos(PartialWheelRot + 1.5*math.pi)		-- Right, Rear.

		-- Piston Speeds:
		-- This checks out according to wikipedia.
		PistSpeed[1] = PartialWheelRotspeed * WheelRot2PistSpeed(PartialWheelRot)
		PistSpeed[2] = PartialWheelRotspeed * WheelRot2PistSpeed(PartialWheelRot + math.pi)
		PistSpeed[3] = PartialWheelRotspeed * WheelRot2PistSpeed(PartialWheelRot + 0.5*math.pi)
		PistSpeed[4] = PartialWheelRotspeed * WheelRot2PistSpeed(PartialWheelRot + 1.5*math.pi)

		-- Piston Accelerations:
		PistAccel[1] = PartialWheelRotspeed * WheelRot2PistAccel(PartialWheelRot)
		PistAccel[2] = PartialWheelRotspeed * WheelRot2PistAccel(PartialWheelRot + math.pi)
		PistAccel[3] = PartialWheelRotspeed * WheelRot2PistAccel(PartialWheelRot + 0.5*math.pi)
		PistAccel[4] = PartialWheelRotspeed * WheelRot2PistAccel(PartialWheelRot + 1.5*math.pi)




		-- Check the valves
		Inlet = InValves(Cranks, PistPos, Valves, cutoff)		-- Returns a vector with 1 for open, 0 for closed.
		Exhaust = ExhValves(Cranks, PistPos, Valves)

		-- Outflow term due to cylinders filling up, using PV = const.
		V_tot = V_SC + A_cyl * (Inlet[1]*PistPos[1] + Inlet[2]*PistPos[2] + Inlet[3]*PistPos[3] + Inlet[4]*PistPos[4]) + sum(Inlet) * V_lost
		dV_tot = A_cyl * (Inlet[1]*PistSpeed[1] + Inlet[2]*PistSpeed[2] + Inlet[3]*PistSpeed[3] + Inlet[4]*PistSpeed[4]) * (dTime/N_SCsim)
		dP_SC_out 	= (P_SC * V_tot)/(V_tot + dV_tot)

		-- Inflow over regulator
		-- This is just a randomly tuned thing.
		dP_SC_in	= 100 * (((P_Boiler + 1) - (P_REG))/(P_Boiler + 1))^0.3 * virtualRegulator^1.5 * (dTime/N_SCsim)

		-- Euler-forward integration.
		P_SC = dP_SC_out + math.max(0, dP_SC_in)

		if P_SC > P_Boiler + 0.9 then
			P_SC = P_Boiler + 0.8
		end



		-- Back Pressure
		-- Some constants:
		A_EXH = 0.01828		-- m^2 Exhaust passage area near valves
		A_BP = 0.0109		-- m^2 Blastpipe area.

		D_EXH =	0.075		-- m Exhaust passage hydraulic radius.
		D_BP = 0.118		-- m Blastpipe Radius.

		L_EXH = 5			-- m Exhaust passage equivalent length.
		L_BP = 0.96			-- m Blaspipe equivalent length.

		-- Exhaust Steam Pressure:
		P_EXH = (P_SC * A_cyl * (cutoff*0.66)) * (A_cyl * 0.6)^-1

		-- Speed of flow through exhaust passage:
		v_EXH =((2 * cutoff* 0.66 * A_cyl)/(6.75 / bound(1, math.abs(speed), 100))) * A_EXH^-1		-- for one passage!

		-- Speed of flow through blastpipe:
		v_BP = ((4 * cutoff * A_cyl)/(6.75 / bound(1 , math.abs(speed), 100))) * A_BP^-1

		-- Reynolds Numbers:
		Re_EXH	= math.max(0.1, (v_EXH * SteamTableDensity(P_EXH) * D_EXH)/0.000015)
		Re_BP	= math.max(0.1, (v_BP * SteamTableDensity(P_EXH) * D_BP)/0.000015)

		-- Friction Coefficients from Moody:
		if Re_EXH > 3000 then
			f_EXH = 0.0055 * (1 + (800 + 1e6/Re_EXH)^(1/3))
		else
			f_EXH = 64/Re_EXH
		end

		if Re_BP > 3000 then
			f_BP = 0.0055 * (1 + (250 + 1e6/Re_BP)^(1/3))
		else
			f_BP = 64/Re_BP
		end

		-- Resulting pressure drop = back pressure:
		-- I'm a huge cheat - it should be 1e-5 but that didn't quite work.
		P_BackPress = 1 + 1e-4 * (L_BP*f_BP*0.5*SteamTableDensity(P_EXH)*(v_BP^2/D_BP)  +  2*L_EXH*f_EXH*0.5*SteamTableDensity(P_EXH)*(v_EXH^2/D_EXH))



		-- Cylinder pressures
		-- Iterate over each cylinder:
		for j = 1,4 do
			if Cranks[j] < math.pi then
				-- Filling up
				if PistPos[j] < cutoff*0.66 then
					-- Intake
					PistPress[j] = P_SC
					PV[j] = PistPress[j] * (V_lost + A_cyl * (cutoff*0.66))
				elseif PistPos[j] > cutoff*0.66 and PistPos[j] < Valves[1]*0.66 then
					-- Adiabatic Expansion
					PistPress[j] = PV[j] * (V_lost + A_cyl * PistPos[j])^-1
					PistReleasePress[j] = PistPress[j]
				else
					-- Exhaust to back pressure.
					PistPress[j] = PistReleasePress[j] + (PistPos[j] - Valves[1]*0.66) * (P_BackPress - PistReleasePress[j])/(0.66 - 0.66*Valves[1])
				end
			else
				-- Exhausting
				if PistPos[j] > Valves[2]*0.66 then
					-- Exhaust under back pressure
					PistPress[j] = P_BackPress
					PV[j] = PistPress[j] * (V_lost + A_cyl * (0.66 - Valves[2]*0.66))
				elseif PistPos[j] < 0.66 - Valves[2]*0.66 and PistPos[j] > 0.66 - Valves[3]*0.66 then
					-- Adiabatic  Compression
					PistPress[j] = PV[j] * (V_lost + A_cyl * PistPos[j])^-1
					PistAdmitPress[j] = PistPress[j]
				elseif PistPos[j] < 0.66 - Valves[3]*0.66 then
					-- Admission
					PistPress[j] = math.min(P_SC, PistAdmitPress[j] + (PistPos[j] - (0.66 - 0.66*Valves[3])) * (PistAdmitPress[j] - P_SC)/(0.66 - 0.66*Valves[3]))
				end
			end
		end



		-- Figure out piston forces and axle torque
		-- Some Constants:
		L_CR = 1.9		-- [m] Length of Connecting Rod
		R_Crank = 0.33	-- [m] Crank
		M_Pist = 250	-- [kg] Weight of piston/crosshead assembly

		-- Forces [kN], positive towards rear.
		PistForce[1] = 1e2 * A_cyl * (PistPress[1] - PistPress[2]) - M_Pist * PistAccel[2] * 1e-3
		PistForce[2] = 1e2 * A_cyl * (PistPress[3] - PistPress[4]) - M_Pist * PistAccel[4] * 1e-3

		-- Torques [kNm]:
		PistTorque[1] = (PistForce[1] / math.cos(math.asin((R_Crank * math.cos(PartialWheelRot + math.pi))/L_CR)))
							* (R_Crank * math.sin((PartialWheelRot + math.pi + math.asin((R_Crank * math.cos(PartialWheelRot + math.pi))/L_CR))))
		PistTorque[2] = (PistForce[2] / math.cos(math.asin((R_Crank * math.cos(PartialWheelRot + 1.5*math.pi))/L_CR)))
							* (R_Crank * math.sin((PartialWheelRot + 1.5*math.pi + math.asin((R_Crank * math.cos(PartialWheelRot + 1.5*math.pi))/L_CR))))

		TrekKracht = (PistTorque[1] + PistTorque[2])/(0.5*2.15)



		-- Average the parameters for stack emitters over the loop
		local iStackEmitVelocity	= 0.5 + 3.00 * ((sum(PistReleasePress) - 4) / 15) * (((PistReleasePress[1]*Exhaust[1] + PistReleasePress[2]*Exhaust[2] + PistReleasePress[3]*Exhaust[3] + PistReleasePress[4]*Exhaust[4])) / 20)^0.7
		local iStackEmitAlpha		= 3 * ((PistPress[1]*Exhaust[1] + PistPress[2]*Exhaust[2] + PistPress[3]*Exhaust[3] + PistPress[4]*Exhaust[4]) / 20)^1.4

		if i == 1 then
			StackEmitVelocity = 0
			StackEmitAlpha = 0
		end

		StackEmitVelocity 	= StackEmitVelocity + (1/N_SCsim)*iStackEmitVelocity
		StackEmitAlpha 		= StackEmitAlpha + (1/N_SCsim)*iStackEmitAlpha



		-- Slip simulation
		-- Weight on Drivers:
		local W_Drivers = 280.566			-- [kN]
		local I_Drivers = 5138.34			-- [kg/m^2]

		-- Factor of adhesion, according to Kofmann:
		local mu_start = 0.24
		local mu = mu_start * (1 - 0.021*speed + 0.0003*speed^2);

		-- Adhesion disappears when slipping:
		if slip then
			mu = (0.90/(50 * (wheelspeed - speed) + 1) + 0.10) * mu 	-- Adhesion drops further the more we slip.
		end

		-- Detect slip:
		-- If either TE > adhesion, or already in a slip:
		if TrekKracht >= mu*W_Drivers or (slip and wheelspeed > 1.05*speed) then
			Call("*:ActivateNode", "LocLantaarnL", 0)
			slip = true
		else
			Call("*:ActivateNode", "LocLantaarnL", 1)
			slip = false
		end

		-- If slipping do a couple of things:
		if slip then
			FricTorque = -(mu*W_Drivers)*(2.15/2)

			dWheelRotSpeed = 1e3*(sum(PistTorque) + FricTorque - (0.7*wheelspeed)^1.0) / (2 * I_Drivers)
			wheelRotspeed = wheelRotspeed + dWheelRotSpeed*(dTime/N_SCsim)
		else
			wheelRotspeed = speed / (2.15/2)
		end

		if Debug then
			CurrentTime = CurrentTime + dTime/N_SCsim
			file:write(CurrentTime .. "," .. PartialWheelRot .. "," .. TrekKracht .. "," .. PistPress[1] .. "," .. PistForce[1] .. "," .. P_SC .. "\n")
		end
	end