%%file roots.py

def quad_roots(a=1.0, b=2.0, c=0.0): 
    """Returns the roots of a quadratic equation: ax^2 + bx + c = 0.
    
    INPUTS
    =======
    a: float, optional, default value is 1
        Coefficient of quadratic term
    b: float, optional, default value is 2
        Coefficient of linear term
    c: float, optional, default value is 0
        Constant term
    
    RETURNS
    ========
    roots: 2-tuple of complex floats
        Has the form (root1, root2) unless a = 0 
        in which case a ValueError exception is raised
    
    EXAMPLES
    =========
    >>> quad_roots(1.0, 1.0, -12.0)
    ((3+0j), (-4+0j))
    """
    import cmath # Can return complex numbers from square roots
    if a == 0:
        raise ValueError("The quadratic coefficient is zero.  This is not a quadratic equation.")
    else:
        sqrtdisc = cmath.sqrt(b * b - 4.0 * a * c)
        r1 = -b + sqrtdisc
        r2 = -b - sqrtdisc
        return (r1 / 2.0 / a, r2 / 2.0 / a)


%%file test_roots.py
import roots

def test_quadroots_result():
    assert roots.quad_roots(1.0, 1.0, -12.0) == ((3+0j), (-4+0j))

def test_quadroots_types():
    try:
        roots.quad_roots("", "green", "hi")
    except TypeError as err:
        assert(type(err) == TypeError)

def test_quadroots_zerocoeff():
    try:
        roots.quad_roots(a=0.0)
    except ValueError as err:
        assert(type(err) == ValueError)

!pytest --cov

"""
Explanation: Homework 5
Due Date:  Tuesday, October 3rd at 11:59 PM
Problem 1
We discussed documentation and testing in lecture and also briefly touched on code coverage.  You must write tests for your code for your final project (and in life).  There is a nice way to automate the testing process called continuous integration (CI).
This problem will walk you through the basics of CI and show you how to get up and running with some CI software.
Continuous Integration
The idea behind continuous integration is to automate away the testing of your code.
We will be using it for our projects.
The basic workflow goes something like this:

You work on your part of the code in your own branch or fork
On every commit you make and push to GitHub, your code is automatically tested on a fresh machine on Travis CI. This ensures that there are no specific dependencies on the structure of your machine that your code needs to run and also ensures that your changes are sane
Now you submit a pull request to master in the main repo (the one you're hoping to contribute to). The repo manager creates a branch off master. 
This branch is also set to run tests on Travis. If all tests pass, then the pull request is accepted and your code becomes part of master.

We use GitHub to integrate our roots library with Travis CI and Coveralls. Note that this is not the only workflow people use. Google git..github..workflow and feel free to choose another one for your group.
Part 1:  Create a repo
Create a public GitHub repo called cs207test and clone it to your local machine.
Note: No need to do this in Jupyter.
Part 2:  Create a roots library
Use the example from lecture 7 to create a file called roots.py, which contains the quad_roots and linear_roots functions (along with their documentation).
Also create a file called test_roots.py, which contains the tests from lecture.
All of these files should be in your newly created cs207test repo.  Don't push yet!!!
End of explanation
"""


%%file .travis.yml
language: python
python:
    - "3.5"
before_install:
    - pip install pytest pytest-cov
script:
    - pytest

%%file setup.cfg
[tool:pytest]
addopts = --doctest-modules --cov-report term-missing --cov roots

"""
Explanation: Part 3:  Create an account on Travis CI and Start Building
Part A:
Create an account on Travis CI and set your cs207test repo up for continuous integration once this repo can be seen on Travis.
Part B:
Create an instruction to Travis to make sure that

python is installed
its python 3.5
pytest is installed

The file should be called .travis.yml and should have the contents:
yml
language: python
python:
    - "3.5"
before_install:
    - pip install pytest pytest-cov
script:
    - pytest
You should also create a configuration file called setup.cfg:
cfg
[tool:pytest]
addopts = --doctest-modules --cov-report term-missing --cov roots
Part C:
Push the new changes to your cs207test repo.
At this point you should be able to see your build on Travis and if and how your tests pass.
End of explanation
"""


%%file .travis.yml
language: python
python:
    - "3.5"
before_install:
    - pip install pytest pytest-cov
    - pip install coveralls
script:
    - py.test
after_success:
    - coveralls


"""
Explanation: Part 4:  Coveralls Integration
In class, we also discussed code coverage.  Just like Travis CI runs tests automatically for you, Coveralls automatically checks your code coverage.  One minor drawback of Coveralls is that it can only work with public GitHub accounts.  However, this isn't too big of a problem since your projects will be public.
Part A:
Create an account on Coveralls, connect your GitHub, and turn Coveralls integration on.
Part B:
Update your the .travis.yml file as follows:
yml
language: python
python:
    - "3.5"
before_install:
    - pip install pytest pytest-cov
    - pip install coveralls
script:
    - py.test
after_success:
    - coveralls
Be sure to push the latest changes to your new repo.
End of explanation
"""


%%file README.md
[![Build Status](https://travis-ci.org/crystalzhaizhai/cs207testing.svg?branch=master)](https://travis-ci.org/crystalzhaizhai/cs207testing.svg?branch=master)

[![Coverage Status](https://coveralls.io/repos/github/crystalzhaizhai/cs207testing/badge.svg?branch=master)](https://coveralls.io/github/crystalzhaizhai/cs207testing?branch=master)

"""
Explanation: Part 5:  Update README.md in repo
You can have your GitHub repo reflect the build status on Travis CI and the code coverage status from Coveralls.  To do this, you should modify the README.md file in your repo to include some badges.  Put the following at the top of your README.md file:
```


```
Of course, you need to make sure that the links are to your repo and not mine.  You can find embed code on the Coveralls and Travis CI sites.
End of explanation
"""


%%file reaction_coeffs.py
from math import exp
def k_constant(k):
    return k
def k_arr(a,e,t):
    r=8.314
    if a<= 0:
        raise ValueError("a<0!")
    if t<= 0:
        raise ValueError("t<0!")

    return (a*exp(-e/(r*t)))
def k_mod_arr(a,e,t,b):
    r=8.314
    if a<= 0:
        raise ValueError("a<0!")
    if t<= 0:
        raise ValueError("t<0!")
    return (a*(t**b)*exp(-e/(r*t)))


import reaction_coeffs as rc
    


a=10**7
b=0.5
e=10**3
t=10**2
print(rc.k_arr(a,e,t))
print(rc.k_mod_arr(a,e,t,b))
print(rc.k_constant(rc.k_arr(a,e,t)))


"""
Explanation: Problem 2
Write a Python module for reaction rate coefficients.  Your module should include functions for constant reaction rate coefficients, Arrhenius reaction rate coefficients, and modified Arrhenius reaction rate coefficients.  Here are their mathematical forms:
\begin{align}
  &k_{\textrm{const}}   = k \tag{constant} \
  &k_{\textrm{arr}}     = A \exp\left(-\frac{E}{RT}\right) \tag{Arrhenius} \
  &k_{\textrm{mod arr}} = A T^{b} \exp\left(-\frac{E}{RT}\right) \tag{Modified Arrhenius}
\end{align}
Test your functions with the following paramters:  $A = 10^7$, $b=0.5$, $E=10^3$.  Use $T=10^2$.
A few additional comments / suggestions:
* The Arrhenius prefactor $A$ is strictly positive
* The modified Arrhenius parameter $b$ must be real 
* $R = 8.314$ is the ideal gas constant.  It should never be changed (except to convert units)
* The temperature $T$ must be positive (assuming a Kelvin scale)
* You may assume that units are consistent
* Document each function!
* You might want to check for overflows and underflows
Recall: A Python module is a .py file which is not part of the main execution script.  The module contains several functions which may be related to each other (like in this problem).  Your module will be importable via the execution script.  For example, suppose you have called your module reaction_coeffs.py and your execution script kinetics.py.  Inside of kinetics.py you will write something like:
```python
import reaction_coeffs
Some code to do some things
:
:
:
Time to use a reaction rate coefficient:
reaction_coeffs.const() # Need appropriate arguments, etc
Continue on...
:
:
:
```
Be sure to include your module in the same directory as your execution script.
End of explanation
"""


%%file progress_rate.py

import numpy as np
import copy
def progress_rate(x,v,vv,kk):
    """Returns progression rate of chemical reactoins
    
    INPUTS
    =======
    x: array of float
        concentration of molecule species
    v: array of integers
        Coefficient of molecule species on the left side of equation
    vv: array of integers
        Coefficient of melecule species on the right side of equation
    kk: array of float
        reaction rate coefficient
    
    RETURNS
    ========
    roots: array of floats
    
    EXAMPLES
    =========
    >>> progress_rate([[1.0],[2.0],[3.0]],[[2.0],[1.0],[0.0]],[[0],[0],[1]],10)
    20.0
    """
    if (kk<0):
        raise ValueError("k<0")
    if all([i[0]<=0 for i in v]):
        raise ValueError("no reactants")
    if all([i[0]<=0 for i in vv]):
        raise ValueError("no products")
    tmp=[x[i][0]**v[i][0] for i in range(len(x))]
    w=copy.copy(kk)
    for i in range(len(x)):
        w*=tmp[i]
#     return (np.transpose([np.array([w*(vv[i][0]-v[i][0]) for i in range(len(x))])]))
    return(w)

%%file progress_test.py
import progress_rate as pr
def test_normal():
    x=[[1.0],[2.0],[3.0]]
    v_i_prime=[[2.0],[1.0],[0.0]]
    v_i_prime_prime=[[0],[0],[1]]
    k=10
    assert pr.progress_rate(x,v_i_prime,v_i_prime_prime,k)==20.0
def test_k_positive():
    try:
        x=[[1.0],[2.0],[3.0]]
        v_i_prime=[[2.0],[1.0],[0.0]]
        v_i_prime_prime=[[0],[0],[1]]
        k=-10
        pr.progress_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)
def test_reactants():
    try:
        x=[[1.0],[2.0],[3.0]]
        v_i_prime=[[0],[0],[0]]
        v_i_prime_prime=[[0],[0],[1]]
        k=10
        pr.progress_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(type(err)==ValueError)

def test_products():
    try:
        x=[[1.0],[2.0],[3.0]]
        v_i_prime=[[2],[1],[0.0]]
        v_i_prime_prime=[[0],[0],[0]]
        k=10
        pr.progress_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(type(err)==ValueError)

import doctest
doctest.testmod(verbose=True)

!pytest

!pytest --cov

"""
Explanation: Problem 3
Write a function that returns the progress rate for a reaction of the following form:
\begin{align}
  \nu_{A} A + \nu_{B} B \longrightarrow \nu_{C} C.
\end{align}
Order your concentration vector so that 
\begin{align}
  \mathbf{x} = 
  \begin{bmatrix}
    \left[A\right] \
    \left[B\right] \
    \left[C\right]
  \end{bmatrix}
\end{align}
Test your function with
\begin{align}
  \nu_{i}^{\prime} = 
  \begin{bmatrix}
    2.0 \
    1.0 \
    0.0
  \end{bmatrix}
  \qquad 
  \mathbf{x} = 
  \begin{bmatrix}
    1.0 \ 
    2.0 \ 
    3.0
  \end{bmatrix}
  \qquad 
  k = 10.
\end{align}
You must document your function and write some tests in addition to the one suggested.  You choose the additional tests, but you must have at least one doctest in addition to a suite of unit tests.
End of explanation
"""


%%file progress1_rate.py
import numpy as np
import copy
def progress1_rate(x,v,vv,k):
    """Returns progress rate of chemical reactoins
    
    INPUTS
    =======
    x: array of float
        concentration of molecule species
    v: array of integers
        Coefficient of molecule species on the left side of equation
    vv: array of integers
        Coefficient of melecule species on the right side of equation
    kk: array of float
        reaction rate coefficient
    
    RETURNS
    ========
    roots: array of floats
    
    EXAMPLES
    =========
    >>> progress1_rate(np.array([[1.0],[2.0],[1]]),np.array([[1,2],[2,0],[0,2]]),np.array([[0,0],[0,1],[2,1]]),[10,10])
    [40.0, 10.0]
    """
    if (any([kk<0 for kk in k])):
        raise ValueError("k<0")
    flag=True
    for j in range(v.shape[1]):
        if all([v[i][j]<=0 for i in range(v.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no reactants")
    flag=True
    for j in range(vv.shape[1]):
        if all([vv[i][j]<=0 for i in range(vv.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no products")
    if x.shape[0]!=v.shape[0] or v.shape!=vv.shape or v.shape[1]!=len(k):
        raise ValueError("dimensions not match")
    tmp=np.array([x[i][0]**v[i][j] for i in range(v.shape[0]) for j in range(v.shape[1])]).reshape(v.shape)
    w=copy.copy(k)
    for i in range(len(x)):
        for j in range(v.shape[1]):
            a=tmp[i][j]
            w[j]*=a
            
    return (w)



import progress1_rate as pr1
import numpy as np
x=np.array([[1.0],[2.0],[1]])
v_i_prime=np.array([[1,2],[2,0],[0,2]])
v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
k=[10,10]
print(pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k))

%%file progress1_test.py
import progress1_rate as pr1
import numpy as np
def test_normal1():
    x=np.array([[1.0],[2.0],[1]])
    v_i_prime=np.array([[1,2],[2,0],[0,2]])
    v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
    k=[10,10]
    assert pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k)==[40.0, 10.0]

def test_nagative_k1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
        k=[-10,10]
        pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(type(err)==ValueError)

def test_reactants1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[0,2],[0,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
        k=[10,10]
        pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)

def test_products1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,0],[2,0]])
        k=[10,10]
        pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)
def test_dimensions1():
    try:
        x=np.array([[1.0],[2.0]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[1,0],[0,1],[2,0]])
        k=[10,10]
        pr1.progress1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)

"""
Explanation: Problem 4
Write a function that returns the progress rate for a system of reactions of the following form:
\begin{align}
  \nu_{11}^{\prime} A + \nu_{21}^{\prime} B \longrightarrow \nu_{31}^{\prime\prime} C \
  \nu_{12}^{\prime} A + \nu_{32}^{\prime} C \longrightarrow \nu_{22}^{\prime\prime} B + \nu_{32}^{\prime\prime} C
\end{align}
Note that $\nu_{ij}^{\prime}$ represents the stoichiometric coefficient of reactant $i$ in reaction $j$ and $\nu_{ij}^{\prime\prime}$ represents the stoichiometric coefficient of product $i$ in reaction $j$.  Therefore, in this convention, I have ordered my vector of concentrations as 
\begin{align}
  \mathbf{x} = 
  \begin{bmatrix}
    \left[A\right] \
    \left[B\right] \
    \left[C\right]
  \end{bmatrix}.
\end{align}
Test your function with 
\begin{align}
  \nu_{ij}^{\prime} = 
  \begin{bmatrix}
    1.0 & 2.0 \
    2.0 & 0.0 \
    0.0 & 2.0
  \end{bmatrix}
  \qquad
  \nu_{ij}^{\prime\prime} = 
  \begin{bmatrix}
    0.0 & 0.0 \
    0.0 & 1.0 \
    2.0 & 1.0
  \end{bmatrix}
  \qquad
  \mathbf{x} = 
  \begin{bmatrix}
    1.0 \
    2.0 \
    1.0
  \end{bmatrix}
  \qquad
  k_{j} = 10, \quad j=1,2.
\end{align}
You must document your function and write some tests in addition to the one suggested.  You choose the additional tests, but you must have at least one doctest in addition to a suite of unit tests.
End of explanation
"""


%%file reaction1_rate.py
import numpy as np
import copy
def reaction1_rate(x,v,vv,k):
    """Returns reaction rate of chemical reactoins
    
    INPUTS
    =======
    x: array of float
        concentration of molecule species
    v: array of integers
        Coefficient of molecule species on the left side of equation
    vv: array of integers
        Coefficient of melecule species on the right side of equation
    kk: array of float
        reaction rate coefficient
    
    RETURNS
    ========
    roots: array of floats
    
    EXAMPLES
    =========
    >>> reaction1_rate(np.array([[1.0],[2.0],[1]]),np.array([[1,0],[2,0],[0,2]]),np.array([[0,1],[0,2],[1,0]]),[10,10])
    array([[-40.],
           [ 10.],
           [-80.],
           [ 20.],
           [ 40.],
           [-20.]])
    """
    if (any([kk<0 for kk in k])):
        raise ValueError("k<0")
    flag=True
    for j in range(v.shape[1]):
        if all([v[i][j]<=0 for i in range(v.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no reactants")
    flag=True
    for j in range(vv.shape[1]):
        if all([vv[i][j]<=0 for i in range(vv.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no products")
    if x.shape[0]!=v.shape[0] or v.shape!=vv.shape or v.shape[1]!=len(k):
        raise ValueError("dimensions not match")
    tmp=np.array([x[i][0]**v[i][j] for i in range(v.shape[0]) for j in range(v.shape[1])]).reshape(v.shape)
    w=copy.copy(k)
    for i in range(len(x)):
        for j in range(v.shape[1]):
            a=tmp[i][j]
            w[j]*=a
            
    return (np.transpose([np.array([w[j]*(vv[i][j]-v[i][j]) for i in range(len(x)) for j in range(v.shape[1])])]))



%%file reaction1_test.py
import reaction1_rate as rr1
import numpy as np
def test_normal1():
    x=np.array([[1.0],[2.0],[1]])
    v_i_prime=np.array([[1,0],[2,0],[0,2]])
    v_i_prime_prime=np.array([[0,1],[0,2],[1,0]])
    k=[10,10]
    assert all(rr1.reaction1_rate(x,v_i_prime,v_i_prime_prime,k)==[[-40],[10],[-80],[20],[40],[-20]])

def test_nagative_k1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
        k=[-10,10]
        rr1.reaction1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(type(err)==ValueError)

def test_reactants1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[0,2],[0,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,1],[2,1]])
        k=[10,10]
        rr1.reaction1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)

def test_products1():
    try:
        x=np.array([[1.0],[2.0],[1]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[0,0],[0,0],[2,0]])
        k=[10,10]
        rr1.reaction1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)
def test_dimensions1():
    try:
        x=np.array([[1.0],[2.0]])
        v_i_prime=np.array([[1,2],[2,0],[0,2]])
        v_i_prime_prime=np.array([[1,0],[0,1],[2,0]])
        k=[10,10]
        rr1.reaction1_rate(x,v_i_prime,v_i_prime_prime,k)
    except ValueError as err:
        assert(True)

import doctest
doctest.testmod(verbose=True)

!pytest --cov

"""
Explanation: Problem 5
Write a function that returns the reaction rate of a system of irreversible reactions of the form:
\begin{align}
  \nu_{11}^{\prime} A + \nu_{21}^{\prime} B &\longrightarrow \nu_{31}^{\prime\prime} C \
  \nu_{32}^{\prime} C &\longrightarrow \nu_{12}^{\prime\prime} A + \nu_{22}^{\prime\prime} B
\end{align}
Once again $\nu_{ij}^{\prime}$ represents the stoichiometric coefficient of reactant $i$ in reaction $j$ and $\nu_{ij}^{\prime\prime}$ represents the stoichiometric coefficient of product $i$ in reaction $j$.  In this convention, I have ordered my vector of concentrations as
\begin{align}
  \mathbf{x} = 
  \begin{bmatrix}
    \left[A\right] \
    \left[B\right] \
    \left[C\right]
  \end{bmatrix}
\end{align}
Test your function with 
\begin{align}
  \nu_{ij}^{\prime} = 
  \begin{bmatrix}
    1.0 & 0.0 \
    2.0 & 0.0 \
    0.0 & 2.0
  \end{bmatrix}
  \qquad
  \nu_{ij}^{\prime\prime} = 
  \begin{bmatrix}
    0.0 & 1.0 \
    0.0 & 2.0 \
    1.0 & 0.0
  \end{bmatrix}
  \qquad
  \mathbf{x} = 
  \begin{bmatrix}
    1.0 \
    2.0 \
    1.0
  \end{bmatrix}
  \qquad
  k_{j} = 10, \quad j = 1,2.
\end{align}
You must document your function and write some tests in addition to the one suggested.  You choose the additional tests, but you must have at least one doctest in addition to a suite of unit tests.
End of explanation
"""


%%file chemkin.py
import numpy as np
import copy
def reaction1_rate(x,v,vv,k):
    """Returns reaction rate of chemical reactoins
    
    INPUTS
    =======
    x: array of float
        concentration of molecule species
    v: array of integers
        Coefficient of molecule species on the left side of equation
    vv: array of integers
        Coefficient of melecule species on the right side of equation
    kk: array of float
        reaction rate coefficient
    
    RETURNS
    ========
    roots: array of floats
    
    EXAMPLES
    =========
    >>> reaction1_rate(np.array([[1.0],[2.0],[1]]),np.array([[1,0],[2,0],[0,2]]),np.array([[0,1],[0,2],[1,0]]),[10,10])
    array([[-40.],
           [ 10.],
           [-80.],
           [ 20.],
           [ 40.],
           [-20.]])
    """
    if (any([kk<0 for kk in k])):
        raise ValueError("k<0")
    flag=True
    for j in range(v.shape[1]):
        if all([v[i][j]<=0 for i in range(v.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no reactants")
    flag=True
    for j in range(vv.shape[1]):
        if all([vv[i][j]<=0 for i in range(vv.shape[0])]):
            flag=False
            break 
    if (flag==False):
        raise ValueError("no products")
    if x.shape[0]!=v.shape[0] or v.shape!=vv.shape or v.shape[1]!=len(k):
        raise ValueError("dimensions not match")
    tmp=np.array([x[i][0]**v[i][j] for i in range(v.shape[0]) for j in range(v.shape[1])]).reshape(v.shape)
    w=copy.copy(k)
    for i in range(len(x)):
        for j in range(v.shape[1]):
            a=tmp[i][j]
            w[j]*=a
            
    return (np.transpose([np.array([w[j]*(vv[i][j]-v[i][j]) for i in range(len(x)) for j in range(v.shape[1])])]))




import numpy as np
import chemkin as ck
x=np.array([[2],[1],[0.5],[1],[1]])
v=np.array([[2,0,0],[1,0,1],[0,1,0],[0,1,0],[0,0,1]])
vv=np.array([[1,0,0],[0,1,0],[2,0,1],[0,0,1],[0,1,0]])
import chemkin as ck
import reaction_coeffs as rc
t=[750,1500,2500]

a1=10**8
b1=0.5
e1=5*10**4

k2=10**4

a3=10**7
e3=10**4


for i in range(len(t)):
    k=[rc.k_mod_arr(a1,e1,t[i],b1),rc.k_constant(k2),rc.k_arr(a3,e3,t[i])]
    print(ck.reaction1_rate(x,v,vv,k))

"""
Explanation: Problem 6
Put parts 3, 4, and 5 in a module called chemkin.
Next, pretend you're a client who needs to compute the reaction rates at three different temperatures ($T = \left{750, 1500, 2500\right}$) of the following system of irreversible reactions:
\begin{align}
  2H_{2} + O_{2} \longrightarrow 2OH + H_{2} \
  OH + HO_{2} \longrightarrow H_{2}O + O_{2} \
  H_{2}O + O_{2} \longrightarrow HO_{2} + OH
\end{align}
The client also happens to know that reaction 1 is a modified Arrhenius reaction with $A_{1} = 10^{8}$, $b_{1} = 0.5$, $E_{1} = 5\times 10^{4}$, reaction 2 has a constant reaction rate parameter $k = 10^{4}$, and reaction 3 is an Arrhenius reaction with $A_{3} = 10^{7}$ and $E_{3} = 10^{4}$.
You should write a script that imports your chemkin module and returns the reaction rates of the species at each temperature of interest given the following species concentrations:
\begin{align}
  \mathbf{x} = 
  \begin{bmatrix}
    H_{2}  \
    O_{2}  \
    OH     \
    HO_{2} \
    H_{2}O
  \end{bmatrix} = 
  \begin{bmatrix}
    2.0 \
    1.0 \
    0.5 \
    1.0 \
    1.0
  \end{bmatrix}
\end{align}
You may assume that these are elementary reactions.
End of explanation
"""


