import math
import matplotlib.pyplot as plt
import numpy as np

def mu1(C1):
    return ((5.5*C1*C1) - (14.7*C1) + 14.6)

def mu2(C2):
    return ((0.55 - (0.3*C2))*C2 + 1.0)

def k1(s):
    return (pow(((s - 0.26) / 0.715),3))

def k2(s):
    return (pow(((0.8 - s)/0.81),3))

def F(s, C1, C2):
    mu = (mu2(C2)/mu1(C1))
    return (k2(s)/(k2(s) + (mu*k1(s))))

K = 2.0
h = 0.01
l = 0.001
n1 = 200
n2 = 1000
C1 = [[0 for j in range(n2)] for i in range(n1)]
C2 = [[0 for j in range(n2)] for i in range(n1)]
s = [[0 for j in range(n2)] for i in range(n1)]

for j in range(n2):
    s[0][j] = 0.0
    C1[0][j] = 1.0
    C2[0][j] = C1[0][j] / K

for i in range(n1):
    s[i][0] = 0.2
    C1[i][0] = 0.0
    C2[i][0] = 0.0

for j in range(0, n2-1):
    for i in range(1, n1):
        C2[i][j] = C1[i][j] / K
        a = ((K - ((K - 1) * F(s[i][j], C1[i][j], C2[i][j])))/(1.0 + (s[i][j]*(K - 1))))
        C1[i][j+1] =  C1[i][j] - ((l / h) * a * (C1[i][j] - C1[i-1][j]))
        s[i][j+1] = s[i][j] + ((l / h)*(F(s[i][j], C1[i][j], C2[i][j]) - F(s[i-1][j], C1[i-1][j], C2[i-1][j])))

for i in range(0, n1):
    print (1-s[i][100])