A simulator for the Mekanism Fusion reactor. Goes in SageMath and makes pretty graphs of temperature and fuel level comparing two different sets of initial conditions.

Reactor simulator

import numpy

class Reactor:
    def __init__(self, alpha, beta, gamma, kappa, Tb, steps):
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.kappa = kappa
        self.Tb = Tb
        self.steps = steps

    def iterate_T(self, T, F):
        if T >= self.Tb:
            dTdt = (min(F + self.alpha, self.beta * (T - self.Tb)))* self.kappa - self.gamma * T
        else:
            dTdt = -self.gamma * T
        return T + dTdt / self.steps

    def iterate_F(self, T, F):
        if T >= self.Tb:
            dFdt = self.alpha - min(F + self.alpha, self.beta * (T - self.Tb))
        else:
            dFdt = 0
        return F + dFdt / self.steps

    def simulate(self, T0, F0, t):
        i = 0
        T = [T0]
        F = [F0]
        F_full = [F0 + self.alpha]
        while i <= t*self.steps:
            T_here = self.iterate_T(T[i], F[i])
            F_here = self.iterate_F(T[i], F[i])

            T = T + [T_here]
            F = F + [F_here]
            F_full = F_full + [F_here + self.alpha]
            i = i + 1
        return T, F_full

    def simulate_enhanced(self, T0, F0, t):
        i = 0
        T = [T0]
        F = [F0]
        F_full = [F0 + self.alpha]
        while i <= t*self.steps:
            T_here = self.iterate_T(T[i], F[i])
            F_here = self.iterate_F(T[i], F[i])

            T_next = self.iterate_T(T_here, F_here)
            F_next = self.iterate_T(T_here, F_here)

            T = T + [1/2 * (T_here + T_next)]
            F = F + [1/2 * (T_here + T_next)]
            F_full = F_full + [1/2 * (F_here + F_next) + self.alpha]
            i = i + 1
        return T, F_full


#Parameters of reactor
#Fuel-injection rate, how many units of fuel are provided per iteration
a=20
#Fuel burn rate, how many units of fuel are burned per iteration per degree above Tb
b=1
#Cooling coefficient. How much of the reactor's temperature is lost at the end of each iteration
g=0.2
#Fuel temperature rate, how many units of temperature are gained per unit fuel burned
k=5000000
#Fuel burn temperature, the temperature at which the fuel starts to burn
Tb=100000000


#Ignition temperature, above which reaction becomes self-sustaining.
Ti = Tb * b*k/(b*k-g)
#Minimum fuel-injection rate for reaction to become stable.
amin = Ti * g/k
#Stable temperature at given injection rate, given a>=amax
Tmax = a*k/g

print Ti, Tmax, amin


steps = 1
reactor = Reactor(a, b, g, k, Tb, steps)


tMax = 100
t = numpy.arange(0, 30, 1/steps)

#This one should *nearly* ignite
T_0 = Ti - 1
F_0 = 0

T1, F1 = reactor.simulate(T_0, F_0, tMax)

#We want this one to ignite
T_0 = Ti + 1
F_0 = 0


T2, F2 = reactor.simulate(T_0, F_0, tMax)


show(list_plot(zip(t, T1), plotjoined=True, rgbcolor=hue(0))+list_plot(zip(t, T2), plotjoined=True))


show(list_plot(zip(t, F1), plotjoined=True, rgbcolor=hue(0))+list_plot(zip(t, F2), plotjoined=True))