import math


def get_r(x, t, xso2):
    R = 1.987
    Keff = 0

    if 693.15 <= t < 748.15:
        Keff = 7.6915 * (10 ** 18) * math.exp(-76062 / (R * t))
    if 748.15 <= t <= 873.15:
        Keff = 1.5128 * (10 ** 7) * math.exp(-35992 / (R * t))

    K = 2.3 * (10 ** (-8)) * math.exp(27200 / (R * t))
    Kp = 2.26203 * (10 ** (-5)) * math.exp(11295.3 / t)

    Pso2 = (xso2 - xso2 * x) / (1 - xso2 * x / 2)
    Pso3 = (xso2 * x) / (1 - xso2 * x / 2)
    Po2 = (0.17 - xso2 - xso2 * x / 2) / (1 - xso2 * x / 2)

    r1 = Po2 * Pso2 / Pso3
    r2 = Pso3 / (Pso2 * math.sqrt(Po2) * Kp)
    B = 48148 * math.exp(-7355.5 / t)

    r3 = math.sqrt(B + (B - 1) * (1 - x) / x) + math.sqrt(K * (1 - x) / x)

    r = Keff * K * r1 * (1 - r2 * r2) / (r3 * r3)

    return r


def get_y(x, t, xso2):
    h = 0.0001
    y = (get_r(x, t + h, xso2) - get_r(x, t - h, xso2)) / (2 * h)
    return y


def get_t(t, x0, x):
    H = -23135
    Cp = 254.9
    rou = 0.500
    c = 1.282

    l = -H * c / (rou * Cp)
    y = t + l * (x - x0)
    return y


def fun1(x, t, xso2):

    y = - get_y(x, t, xso2) / (get_r(x, t, xso2) * get_r(x, t, xso2))
    return y


def jifen(x, t0, xso2):
    sum = 0.0
    x1 = x
    t2 = 693.15
    h = 0.0001
    xout = 0

    t1 = get_t(t0, x, x1)
    x2 = x1 + h / 10
    t2 = get_t(t0, x, x2)

    sum += h * (fun1(x1, t1, xso2) + fun1(x2, t2, xso2)) / 20
    x1 = x2
    while sum < 0:
        t1 = get_t(t0, x, x1)
        x2 = x1 + h / 10
        t2 = get_t(t0, x, x2)

        if t2 > 873.15:
            xout = x1
            break
        sum += h * (fun1(x1, t1, xso2) + fun1(x2, t2, xso2)) / 20
        x1 = x2
        xout = x1 - h / 10
    return xout


def weijifen(xin, xou, tin, xso2):
    x1 = xin
    sum = 0.0
    h = 0.0001

    t1 = get_t(tin, xin, x1)
    x2 = x1 + h
    t2 = get_t(tin, xin, x2)

    sum += (1 / get_r(x1, t1, xso2) + 1 / get_r(x2, t2, xso2)) * h / 1000
    x1 = x2

    while x2 <= xou:
        t1 = get_t(tin, xin, x1)
        x2 = x1 + h
        t2 = get_t(tin, xin, x2)
        if t2 >= 873.15:
            break
        else:
            sum += (1 / get_r(x1, t1, xso2) + 1 / get_r(x2, t2, xso2)) * h / 1000
            x1 = x2
    wcat = sum * 131 * 1000 / 3600
    return wcat


def write_2_file():
    xso2 = 0.07
    x0 = 0.0001
    t00 = 717
    wsum = 0.0

    xou = [0] * 5
    tou = [0] * 5
    xin = [0] * 5
    tin = [0] * 5

    with open('data-3.txt', 'w') as f:

        while True:
            h = 0.0001
            x0 = 0.0001
            t0 = t00
            i = 0

            xin[i] = x0
            tin[i] = t00

            while True:
                xout = jifen(x0, t0, xso2)
                tout = get_t(t0, x0, xout)
                xou[i] = xout
                tou[i] = tout
                t1 = 693.15
                t1 += 0.01

                while math.fabs((10 ** 5) * get_r(xout, t1, xso2) - (10 ** 5) * get_r(xout, tout, xso2)) > h:
                    t1 += 0.01

                x0 = xout
                t0 = t1
                i += 1

                xin[i] = xout
                tin[i] = t0

                if i > 3:
                    break
            t00 -= 0.1
            if x0 > 0.98:
                break
        f.write('i    tin        xin        tout        xout         wcat\n')

        for j in range(4):
            wcat = weijifen(xin[j], xou[j], tin[j], xso2)
            f.write('%d   %5.5f   %5.5f   %5.5f   %5.5f    %5.5f\n' % (j + 1, tin[j], xin[j], tou[j], xou[j], wcat))
            wsum += wcat

        f.write('wsum = %f (KG) %f (Ton)\n' % (wsum, wsum / 1000))


def main():
    write_2_file()


if __name__ == '__main__':
    main()

+++++++++++++++++++++
运行结果：
+++++++++++++++++++++

i    tin        xin        tout        xout         wcat
1   717.00000   0.00010   873.14936   0.67110    4648.31851
2   738.42000   0.67110   790.77774   0.89609    4167.40302
3   732.77000   0.89609   747.91253   0.96116    5352.48512
4   708.46000   0.96116   712.88384   0.98017    12312.46358
wsum = 26480.670240 (KG) 26.480670 (Ton)