#Vegard Jervell, 26.01-17
#Dette scriptet løser oppgave1 i øving 1 numerisk ved å bruke scipy.optimize.root() metoden

import scipy.optimize as sciopt

def modellPropen(x_list):

    #definerer kjennte variabler

    n1_c2h8 = 100       #fødestrøm

    x3_c2h8 = 0.01      #molbrøker ut av separator (toppstrøm)
    x3_c2h6 = 0.9
    x3_h2 = 0.99

    x4_c2h8 = 0.99      #molbrøker ut av separator (bunnstrøm)
    x4_c2h6 = 0.1
    x4_h2 = 0.01

    X = 0.8     #reaksjonsomfang

    #definerer ukjente variabler

    n2_c2h8 = x_list[0]     #ut av reaktor
    n2_c2h6 = x_list[1]
    n2_h2 = x_list[2]

    n3_c2h8 = x_list[3]     #toppstrøm
    n3_c2h6 = x_list[4]
    n3_h2 = x_list[5]

    n4_c2h8 = x_list[6]     #bunnstrøm
    n4_c2h6 = x_list[7]
    n4_h2 = x_list[8]

    #setter opp likninger over reaktor
    eq1 = (1-X)*n1_c2h8 - n2_c2h8
    eq2 = X*n1_c2h8 - n2_c2h6
    eq3 = X*n1_c2h8 - n2_h2

    #setter opp likninger over separator (toppstrøm)
    eq4 = n3_c2h8 - x3_c2h8*n2_c2h8
    eq5 = n3_c2h6 - x3_c2h6*n2_c2h6
    eq6 = n3_h2 - x3_h2*n2_h2

    #bunnstrøm
    eq7 = n4_c2h8 - x4_c2h8*n2_c2h8
    eq8 = n4_c2h6 - x4_c2h6*n2_c2h6
    eq9 = n4_h2 - x4_h2*n2_h2

    balance = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9]

    return balance

#gjettninger på forskjellige variable
Gn2_c2h8 = 20  # ut av reaktor
Gn2_c2h6 = 80
Gn2_h2 = 80

Gn3_c2h8 = 0.2  # toppstrøm
Gn3_c2h6 = 70
Gn3_h2 = 79

Gn4_c2h8 = 19  # bunnstrøm
Gn4_c2h6 = 10
Gn4_h2 = 1

Glist = [Gn2_c2h8,Gn2_c2h6,Gn2_h2,Gn3_c2h8,Gn3_c2h6,Gn3_h2,Gn4_c2h8,Gn4_c2h6,Gn4_h2]

#Løser likningssett
ans = sciopt.root(modellPropen,Glist)

print(ans)