TG Atmos simulation

import turtle
from math import floor
from time import perf_counter

G = dict(P=0.0,O2=0.0,CO2=0.0)
C = dict(P=200.0,O2=20.0,CO2=30.0)
T = 0.0
last_results = []

MHC = 0.0003
PUT = 1370+273.15
PMBT = 100+273.15
POFB = 10.0
FPER = 3000000.0

def heat_capacity():
    return G["P"]*C["P"]+G["O2"]*C["O2"]+G["CO2"]*C["CO2"]

def energy():
    return heat_capacity()*T

tick = 0
plotting = turtle.Turtle(visible=False)
turtle.tracer(None, None)

def zburn():
    global G,T,tick

    tick = 0
    the_process = []

    while (True):
        fuel_burnt = 0
        energy_released = 0
        old_heat_capacity = heat_capacity()
        if (G["P"] > MHC):
            PBR = 0
            OBR = 0
            TS = 0
            if (T > PUT):
                TS = 1
            else:
                TS = (T-PMBT)/(PUT-PMBT)
            if (TS >= 0.0):
                OBR = 1.4 - TS
                if (G["O2"] > G["P"]*POFB):
                    PBR = 0.25 * G["P"] * TS
                else:
                    PBR = 0.25 * TS * (G["O2"] / POFB)
                if (PBR > MHC):
                    G["P"] -= PBR
                    G["O2"] -= PBR*OBR
                    G["CO2"] += PBR
                    energy_released += FPER*PBR
                    fuel_burnt += PBR*(1+OBR)
            if (energy_released > 0):
                NHC = heat_capacity()
                if (NHC > MHC):
                    T = (T*old_heat_capacity + energy_released) / NHC
                    the_process.append([G["P"], G["O2"], G["CO2"], T, energy_released])
        if (fuel_burnt > 0):
            tick += 1
            continue
        return the_process


PFLAG = 1<<0
OFLAG = 1<<1
CFLAG = 1<<2
TFLAG = 1<<3
E0FLAG = 1<<4
E1FLAG = 1<<5
DFLAG = 1<<6
PLOTFLAG = 1<<7

def make_sim(p0,o0,t,modes, give_process=0,c2=0):
    global T,G,tick,last_results
    G = dict(P=p0,O2=o0,CO2=c2)
    T = float(t)
    energy0 = energy()

    if (modes & PFLAG):
        print(round(G["P"], 3), 'mol', end='')
    if (modes & OFLAG):
        print(' :', end=' ')
        print(round(G["O2"], 3), 'mol', end='')
    if (modes & E0FLAG):
        print(' :', end=' ')
        print(floor(energy0), 'J', end='')
    last_results += G["P"], G["O2"], energy0
    the_process = zburn()
    if (modes and modes != PLOTFLAG): print(' --> ', end='')
    if (modes & CFLAG):
        print(round(G["CO2"], 3), 'mol', end='')
    if (modes & TFLAG):
        print(' :', end=' ')
        print(round(T, 3), 'K', end='')
    if (modes & E1FLAG):
        print(' :', end=' ')
        print(floor(energy()-energy0), 'J', end='')
    if (modes & DFLAG):
        print(' :', end=' ')
        print(tick/2, 'seconds (normal conditions)', end='')
    if (modes & PLOTFLAG):
        maxTemp = 0
        for tick in range(len(the_process)):
            if the_process[tick][3] > maxTemp:
                maxTemp = the_process[tick][3]
        turtle.setworldcoordinates(0, 0, len(the_process), maxTemp)
        turtle.penup()
        turtle.goto(0, t)
        turtle.pendown()
        for tick in range(len(the_process)):
            turtle.goto(tick, the_process[tick][3])
    last_results += G["CO2"], T, energy()-energy0, tick
    if (give_process):
        return the_process

TMODE = 1<<0
EMODE = 1<<1
SPMODE = 1<<2
SNMODE = 1<<3

def find_ideal_ratio_for (mode, t=393.15, pconst=10, oupconst=1000):
    global last_results
    des_val = 0
    ideal_ratio = 0
    if (mode in (TMODE, EMODE)):
        for x in range(oupconst):
            curO2 = pconst*x/oupconst
            make_sim(pconst, curO2, t, 0)
            cur_val = 0
            if (mode == TMODE): cur_val = last_results[4]
            if (mode == EMODE): cur_val = last_results[5]
            if (cur_val > des_val):
                ideal_ratio = curO2 / (curO2+pconst)
                des_val = cur_val
            last_results = []
    if (mode in (SPMODE, SNMODE)):
        if (mode == SNMODE): des_val = 9e40
        for x in range(15*pconst): #Getting higher than 1:15 is silly
            curO2 = pconst*x/15
            make_sim(pconst, curO2, t, 0)
            cur_val = last_results[6]
            if (mode == SNMODE):
                if (cur_val < des_val and cur_val != 0):
                    ideal_ratio = curO2 / (curO2+pconst)
                    des_val = cur_val
            if (mode == SPMODE):
                if (cur_val > des_val and cur_val != 0):
                    ideal_ratio = curO2 / (curO2+pconst)
                    des_val = cur_val
            last_results = []
    print ('Ideal ratio is ', round(100*(1-ideal_ratio), 3), '%:', round(100*ideal_ratio, 3),'%', sep='')
    print ('Desired value achieved:', floor(des_val))

def find_var(var):
    if (var.upper() == "N"):
        # 8.31 = pV / nT
        # 1/8.31 = nT / pV
        # pV / 8.31T
        # p = 8.3nT / V
        print ("p V T")
    if (var.upper() == "P"):
        print ("n T V")

    args = input()
    args = args.split(",")

    for x in range(len(args)):
        args[x] = float(args[x])

    if (var.upper() == "N"):
        return (args[0] * args[1]) / (8.31 * args[2])
    if (var.upper() == "P"):
        return (8.31 * args[0] * args[1]) / args[2]

ONEATMOS = 101.325

def simulate_explosion(p0, t0, cp0, p1, t1, give_results=0, verbose=1, debug=0):
    results = []
    RAmt = 70*p0 / (8.31*t0)
    PAmt = RAmt*cp0
    CAmt = RAmt-PAmt
    OAmt = 70*p1 / (8.31*t1)
    # Merging this shit together
    C0 = 200*PAmt + 30*CAmt
    C1 = 20*OAmt
    C2 = C0+C1
    E0 = C0*t0
    E1 = C1*t1
    E2 = E0+E1
    simul = list()
    # E = CnT
    T = E2 / C2
    amt1 = RAmt+OAmt
    p1 = 8.31*amt1*T / 140
    if verbose:
        print("Merging temp:     ", T)
        print("Merging pressure: ", p1)
    if (T > PMBT):
        simul = make_sim(PAmt, OAmt, T, 0, 1, CAmt)
    ft = 0
    td = 0
    p2 = p1
    integrity = 3
    for data in simul:
        amt2 = data[0]+data[1]+data[2]
        t2   = data[3]
        p2   = 8.31*amt2*t2 / 140
        if (p2 > 50*ONEATMOS):
            ft += 3
            break
        elif (p2 > 40*ONEATMOS):
            if (integrity <= 0):
                if verbose:
                    print('Tank will rupture without explosion after {0} seconds'.format(td))
                return False
            else:
                integrity -= 1
        elif (p2 > 30*ONEATMOS):
            if (integrity <= 0):
                if verbose:
                    print('Tank will leak without explosion after {0} seconds'.format(td))
                return False
            else:
                integrity -= 1
        elif (integrity < 3):
            integrity += 1
        ft  += 1
        td  += 0.5
    ft = max(min(len(simul)-1, ft), 0)
    if (ft):
        EMD = simul[ft] # Explosion Moment Data
        if verbose:
            print("EMD:")
            print("Plasma:           ", EMD[0])
            print("Oxygen:           ", EMD[1])
            print("CO2:              ", EMD[2])
            print("T:                ", EMD[3])
        amt2 = EMD[0]+EMD[1]+EMD[2]
        t2  = EMD[3]
        p2 = 8.31*amt2*t2 / 140
        radius = (p2-50*ONEATMOS)/(10*ONEATMOS)
        results.append(radius)
        results.append(floor(td))
        if verbose:
            print(floor(0.25*radius), floor(0.5*radius), floor(radius), sep=', ')
            print("Will explode after {0} seconds".format(td))
        if (give_results):
            return results
        return True
    if verbose:
        print("Won\'t explode")
    return False

def perform_cross_simulation(pw, pt, o2t, res=10):
    st = perf_counter()

    cold_gas_maxamt = 70*10*ONEATMOS / (8.31*293.15)
    hot_plsm_maxamt  = 70*10*ONEATMOS / (8.31*pt)
    cold_gas_mult   = cold_gas_maxamt / res
    hot_plsm_mult    = hot_plsm_maxamt / res

    co2mult = 30*(1-pw)
    plsmmult = 200*pw
    mixmult = co2mult+plsmmult

    cold_energy_mult = 293.15*200
    hot_energy_mult = pt*mixmult

    bHPAmt = 0
    bCPAmt = 0
    bO2Amt = 0
    bTemp  = 0
    bRad   = 0
    bTime  = 0

    for cold_plsm_iter in range(1, res+1):
        cold_plsm_amt = cold_gas_mult*cold_plsm_iter
        actual_cold_energy = cold_plsm_amt * cold_energy_mult
        for hot_plsm_iter in range(1, res+1):
            hot_plsm_amt = hot_plsm_mult*hot_plsm_iter
            actual_hot_energy = hot_plsm_amt * hot_energy_mult
            merging_temp = (actual_hot_energy + actual_cold_energy) / (hot_plsm_amt * mixmult + 200 * cold_plsm_amt)
            merging_pres = 8.31 * (hot_plsm_amt+cold_plsm_amt) * merging_temp / 70
            if (merging_pres > 10*ONEATMOS):
                break
            for cold_oxygen_iter in range(1, res+1):
                cold_oxygen_amt = cold_gas_mult*cold_oxygen_iter
                oxygen_pres = 8.31*cold_oxygen_amt*o2t / 70
                actual_pw = (cold_plsm_amt+(hot_plsm_amt*pw)) / (cold_plsm_amt+hot_plsm_amt)
                result = simulate_explosion(merging_pres, merging_temp, actual_pw, oxygen_pres, o2t, 1, 0)
                if result:
                    if result[0] > bRad:
                        bHPAmt = hot_plsm_amt
                        bTemp  = merging_temp
                        bCPAmt = cold_plsm_amt
                        bO2Amt = cold_oxygen_amt
                        bRad = result[0]
                        bTime = result[1]
    bHPAmt = 8.31 * (bCPAmt+bHPAmt)*bTemp / 70
    bCPAmt = 8.31 * bCPAmt * 293.15 / 70
    bO2Amt = 8.31 * bO2Amt * o2t / 70
    print("Statistics:")
    print("Cold Plasma Pressure: {0}\nHot Plasma Pressure: {1}\nO2 Pressure: {2}\nTemp.: {3}".format(bCPAmt, bHPAmt, bO2Amt, bTemp))
    print("Expected radius:")
    print("{0} {1} {2}".format(floor(0.25*bRad), floor(0.5*bRad), floor(bRad)))
    print("True radius:")
    print(bRad)
    print("Will take {0} seconds before explosion".format(bTime))
    print(perf_counter() - st)