from scipy.optimize import root
import numpy as np

def model(x_list):
    d_ck, d_ch = x_list

    c0_kl = 0.1
    c0_hl = 0.1
    c0_cl = 0.2

    c0_kr = 0.1
    c0_hr = 0
    c0_cr = 0.1

    mu_dh = 0.727

    c_kl = c0_kl - d_ck
    c_hl = c0_hl - d_ch

    c_kr = c0_kr + d_ck
    c_hr = c0_hr + d_ch

    #mu_l = 0.5 * (c_kl + c_hl + c0_cl)
    #mu_r = 0.5 * (c_kr + c_hr + c0_cr)
    #eq1 = (np.sqrt(mu_r) - np.sqrt(mu_l))/np.sqrt(mu_dh) - np.log((c0_kr + d_ck)/(c0_kl - d_ck))
    #eq1 = (np.sqrt(mu_r) - np.sqrt(mu_l))/np.sqrt(mu_dh) - np.log((c0_hr + d_ch)/(c0_hl - d_ch))
    #eq1 = np.log((c0_hr + d_ch)/(c0_hl - d_ch)) -  np.log((c0_kr + d_ck)/(c0_kl - d_ck))
    
    eq1 = ((c0_hr + d_ch) * (c0_kl - d_ck)) / ((c0_hl - d_ch) * (c0_kr + d_ck)) - 1
    eq2 = c_hl + c_kl - 0.2

    return [eq1, eq2]

guess = [-0.05, 0.05]
ans = root(model, guess)

print(ans)

ck_l = 0.1 - ans.x[0]
ck_r = 0.1 + ans.x[0]
ch_l = 0.1 - ans.x[1]
ch_r = ans.x[1]

mu_l = 0.5 * (ck_l + ch_l + 0.2)
mu_r = 0.5 * (ck_r + ch_r + 0.1)
mu_dh = 0.727

print('d_mu_k =', np.exp(-np.sqrt(mu_l/mu_dh))*ck_l - np.exp(-np.sqrt(mu_r/mu_dh))*ck_r)
print('d_mu_h =', np.exp(-np.sqrt(mu_l/mu_dh))*ch_l - np.exp(-np.sqrt(mu_r/mu_dh))*ch_r)
print('d_mu_cl =', np.exp(-np.sqrt(mu_l/mu_dh))*0.2 - np.exp(-np.sqrt(mu_r/mu_dh))*0.1)

print('C_l =', ck_l + ch_l)
print('C_r =', ck_r + ch_r)

print('End concentrations: ')
print('|       L       |       R       |')
print('| C_K+ =', round(0.1 - ans.x[0], 4), '| C_K+ =', round(0.1 + ans.x[0], 4), '|')
print('| C_H+ =', round(0.1 - ans.x[1], 4), '| C_H+ =', round(ans.x[1], 4), '|')