hh52.py

from math import *

import numpy as np

def initializeConstants(states, variables):
    states[0] = 0.05
    states[1] = 0.6
    states[2] = 0.325
    states[3] = 0.0
    variables[0] = 0.3
    variables[1] = 1.0
    variables[2] = 0.0
    variables[3] = 36.0
    variables[4] = 120.0

def computeComputedConstants(variables):
    variables[6] = variables[2]-10.613
    variables[8] = variables[2]-115.0
    variables[14] = variables[2]+12.0

def computeRates(voi, states, rates, variables):
    variables[10] = 0.1*(states[3]+25.0)/(exp((states[3]+25.0)/10.0)-1.0)
    variables[11] = 4.0*exp(states[3]/18.0)
    rates[0] = variables[10]*(1.0-states[0])-variables[11]*states[0]
    variables[12] = 0.07*exp(states[3]/20.0)
    variables[13] = 1.0/(exp((states[3]+30.0)/10.0)+1.0)
    rates[1] = variables[12]*(1.0-states[1])-variables[13]*states[1]
    variables[16] = 0.01*(states[3]+10.0)/(exp((states[3]+10.0)/10.0)-1.0)
    variables[17] = 0.125*exp(states[3]/80.0)
    rates[2] = variables[16]*(1.0-states[2])-variables[17]*states[2]
    variables[5] = -20.0 if ((voi >= 10.0) & (voi <= 10.5)) else 0.0
    variables[7] = variables[0]*(states[3]-variables[6])
    variables[15] = variables[3]*pow(states[2], 4.0)*(states[3]-variables[14])
    variables[9] = variables[4]*pow(states[0], 3.0)*states[1]*(states[3]-variables[8])
    rates[3] = -(-variables[5]+variables[9]+variables[15]+variables[7])/variables[1]

def computeVariables(voi, states, rates, variables):
    variables[7] = variables[0]*(states[3]-variables[6])
    variables[9] = variables[4]*pow(states[0], 3.0)*states[1]*(states[3]-variables[8])
    variables[10] = 0.1*(states[3]+25.0)/(exp((states[3]+25.0)/10.0)-1.0)
    variables[11] = 4.0*exp(states[3]/18.0)
    variables[12] = 0.07*exp(states[3]/20.0)
    variables[13] = 1.0/(exp((states[3]+30.0)/10.0)+1.0)
    variables[15] = variables[3]*pow(states[2], 4.0)*(states[3]-variables[14])
    variables[16] = 0.01*(states[3]+10.0)/(exp((states[3]+10.0)/10.0)-1.0)
    variables[17] = 0.125*exp(states[3]/80.0)

NB_OF_STATES = 4
NB_OF_VARIABLES = 18
DELTA_T = 0.01

def outputData(voi, states, variables):
    print(f"{voi}", end="")

    for state in states:
        print(f",{state}", end="")

    for variable in variables:
        print(f",{variable}", end="")

    print("")

states = np.zeros(NB_OF_STATES)
rates = np.zeros(NB_OF_STATES)
variables = np.zeros(NB_OF_VARIABLES)

initializeConstants(states, variables)
computeComputedConstants(variables)
computeRates(0.0, states, rates, variables)
computeVariables(0.0, states, rates, variables)

outputData(0.0, states, variables)

for i in range(5000+1):
    voi = i*DELTA_T

    computeRates(voi, states, rates, variables)

    for j in range(NB_OF_STATES):
        states[j] += DELTA_T*rates[j]

    if fmod(i+1, 10) == 0:
        computeVariables(voi, states, rates, variables)

        outputData((i+1)*DELTA_T, states, variables)