<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Hydrodynamics Systems with a General EOS
This notebook explores the eigensystem of the Euler equations augmented with an additional thermodynamic variable to describe a general equation of state in the reconstruction of interface states
Step1:  The two routines below simplify the analysis of the eigensystem by simultaneously finding the set of orthonormal left and right eigenvectors of a matrix from the primitive form of the hydro equations
Step2:  Euler Equations with $(\rho e)$
The Euler equations in primitive variable form, $q = (\rho, u, p, (\rho e))^\intercal$ appear as
Step3:  The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
Step4:  We see that there are 2 eigenvalues $u$ -- the addition of $(\rho e)$ to the system adds this degeneracy.
Eigenvectors
The right eigenvectors are defined for a given eigenvalue, $\lambda$, as
Step5:  $\beta$'s and final update
The final interface state is writen by projecting the jump in primitive variables, $\Delta q$, into characteristic variables (as $l \cdot \Delta q$), and then adding up all the jumps that reach the interface.
The convention is to write $\beta^\nu = l^\nu \cdot \Delta q$, where the superscript identifies which eigenvalue (and corresponding eigenvectors) we are considering.  Note, that often a reference state is used, and the jump, $\Delta q$, will be the difference with respect to this reference state.  For PPM, the $\Delta q$ will take the form of the integral under the parabola over the range that each wave can reach.
The final interface state is then
Step6:  Now compute the $\beta$s
Step7:  and now the final interface states
Step8:  Euler Equations with $(\gamma_e)$
We can define $\gamma_e = p/(\rho e) + 1$ and differentiate it to get
Step9:  The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
Step10:  We see that there are 2 eigenvalues $u$ -- the addition of $\gamma_e$ to the system adds this degeneracy.
Eigenvectors
Step11:  $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$$
\Delta q = \left ( \begin{array}{c} \Delta \tau \ \Delta u \ \Delta p \ \Delta \gamma_e \end{array} \right )
$$
and then
$\beta^\nu = l^\nu \cdot \Delta q$
Step12:  and now the final interface states
Step13:  Gray FLD Radiation Euler Equations with $(\rho e)$
Following Zhang et al. (2011), the equations of gray radiation hydrodynamics with primitive
variables $q = (\rho, u, p, (\rho e)_g, E_r)^\intercal$ are
Step14:  The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
Step15:  We see that there are 3 eigenvalues $u$.  We identify the total sound speed (radiation + gas) as
Step16:  We see that these are the same form of the eigenvalues for the pure hydrodynamics system
Eigenvectors
Step17:  $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$\beta^\nu = l^\nu \cdot \Delta q$
Step18:  and now the final interface states
Step19:  Gray FLD Radiation Euler Equations with $(\gamma_e)$
We now look at the same system with a different auxillary thermodynamic variable (as we did with pure hydro), using $q = (\tau, u, p, {\gamma_e}_g, E_r)^\intercal$
Step20:  The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
Step21:  We see that there are 3 eigenvalues $u$.  We identify the total sound speed (radiation + gas) as
Step22:  We see that these are the same form of the eigenvalues for the pure hydrodynamics system
Eigenvectors
Step23:  $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$\beta^\nu = l^\nu \cdot \Delta q$
Step24:  and now the final interface states
<ASSISTANT_TASK:>
Python Code:
from sympy import init_session
init_session()
from sympy.abc import rho, tau, alpha
rho, tau, c, h, p = symbols("rho tau c h p", real=True, positive=True)
re = symbols(r"(\rho{}e)", real=True, positive=True)
ge = symbols(r"\gamma_e", real=True, positive=True)
alpha, u = symbols("alpha u", real=True)
Explanation: Hydrodynamics Systems with a General EOS
This notebook explores the eigensystem of the Euler equations augmented with an additional thermodynamic variable to describe a general equation of state in the reconstruction of interface states
End of explanation
plus = u + c
zero = u
minus = u - c
class Eigenvector(object):
    def __init__(self, name, ev, r, l=None):
        self.name = name
        if name == "minus":
            self.d = 0
        elif name == "zero":
            self.d = 1
        elif name == "plus":
            self.d = 2
        else:
            self.d = None
        self.ev = ev
        self.l = l
        self.r = r
    def __lt__(self, other):
        return self.d < other.d
    def __str__(self):
        return "{} wave, r = {}, l = {}".format(self.eigenvalue, self.r, self.l)
def eigensystem(A, suba=None, subb=None):
    # get the left and right eigenvectors that diagonalize the system.
    # it is best to use sympy diagonalize() for this purpose than getting
    # the left and right eigenvectors independently.
    e = []
    R, D = A.diagonalize()
    # the columns of R are the right eigenvectors and the diagonal
    # element of D is the corresponding eigenvalues
    for n in range(A.shape[0]):
        r = R.col(n)
        ev = D[n,n]
                    
        #print("here", r, ev)
        if suba is not None and subb is not None:
            ev = ev.subs(suba, subb)
            
        # which eigenvalue are we?
        if simplify(ev - minus) == 0:
            name = "minus"
        elif simplify(ev - plus) == 0:
            name = "plus"
        elif simplify(ev - zero) == 0:
            name = "zero"
        else:
            return None
                                                                                                            
        # normalize the right eigenvector
        v = r[0]
        if v != 0:
            r = r/v
        if suba is not None and subb is not None:
            r = simplify(r.subs(suba, subb))
            
        e.append(Eigenvector(name=name, ev=ev, r=r))
    # now sort the system from smallest (u-c) to largest (u+c)
    e.sort()
    # now let's construct the R with this sorting
    for n in range(A.shape[0]):
        R[:,n] = e[n].r
    # the left eigenvector matrix, L, is just the inverse
    L = R**-1
    for n in range(A.shape[0]):
        e[n].l = L.row(n)
    return e
Explanation: The two routines below simplify the analysis of the eigensystem by simultaneously finding the set of orthonormal left and right eigenvectors of a matrix from the primitive form of the hydro equations
End of explanation
q = Matrix([rho, u, p, re]).transpose()
A = Matrix([[u, rho, 0, 0], [0, u, rho**-1, 0], [0, c**2 * rho, u, 0], [0, rho*h, 0, u]])
A
Explanation: Euler Equations with $(\rho e)$
The Euler equations in primitive variable form, $q = (\rho, u, p, (\rho e))^\intercal$ appear as:
\begin{align}
\frac{\partial \rho}{\partial t} &= -u \frac{\partial \rho}{\partial x} - \rho \frac{\partial u}{\partial x} \
\frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \
\frac{\partial (\rho e)}{\partial t} &= -u \frac{\partial (\rho e)}{\partial x} - \rho h \frac{\partial u}{\partial x} 
\end{align}
In vector form, we have:
$$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{cccc} u  & \rho     & 0    & 0\                        
                                  0  &  u       & 1/\rho  & 0\                   
                                  0  & \rho c^2 & u       & 0\
                                  0  & \rho h & 0       & u\end{array} \right )
$$
The sound speed is related to the adiabatic index, $\Gamma_1$, as $c^2 = \Gamma_1 p /\rho$.
We can represent this matrix symbolically in SymPy and explore its eigensystem.
End of explanation
A.eigenvals()
Explanation: The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
End of explanation
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
    print(e.name)
    display(e.r, e.l)
Explanation: We see that there are 2 eigenvalues $u$ -- the addition of $(\rho e)$ to the system adds this degeneracy.
Eigenvectors
The right eigenvectors are defined for a given eigenvalue, $\lambda$, as:
$$A r = \lambda r$$
and the left eigenvectors satisfy:
$$l A = \lambda l$$
Note that the left and right eigenvectors are orthogonal to those corresponding to a different eigenvalue, and usually normalized so:
$$l^i \cdot r^j = \delta_{ij}$$
End of explanation
from sympy.abc import delta
dr = symbols(r"\Delta\rho")
du = symbols(r"\Delta{}u")
dp = symbols(r"\Delta{}p")
dre = symbols(r"\Delta(\rho{}e)")
rhoi = symbols(r"\rho_\mathrm{int}")
ui = symbols(r"u_\mathrm{int}")
pri = symbols(r"p_\mathrm{int}")
rei = symbols(r"(\rho{}e)_\mathrm{int}")
# this is the jump
dq = Matrix([[dr, du, dp, dre]]).transpose()
# this is the interface state vector
qint = Matrix([[rhoi, ui, pri, rei]]).transpose()
Explanation: $\beta$'s and final update
The final interface state is writen by projecting the jump in primitive variables, $\Delta q$, into characteristic variables (as $l \cdot \Delta q$), and then adding up all the jumps that reach the interface.
The convention is to write $\beta^\nu = l^\nu \cdot \Delta q$, where the superscript identifies which eigenvalue (and corresponding eigenvectors) we are considering.  Note, that often a reference state is used, and the jump, $\Delta q$, will be the difference with respect to this reference state.  For PPM, the $\Delta q$ will take the form of the integral under the parabola over the range that each wave can reach.
The final interface state is then:
$$
q_\mathrm{int} = q - \sum_\nu \beta^\nu r^\nu
$$
The tracing projects the primitive variables into characteristic variables by defining
$$
\Delta q = \left ( \begin{array}{c} \Delta \rho \ \Delta u \ \Delta p \ \Delta (\rho e) \end{array} \right )
$$
and then
$\beta^\nu = l^\nu \cdot \Delta q$
End of explanation
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))
Explanation: Now compute the $\beta$s
End of explanation
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))
Explanation: and now the final interface states
End of explanation
q = Matrix([[tau, u, p, ge]]).transpose()
A = Matrix([[u, -tau, 0, 0], [0, u, tau, 0], [0, c**2/tau, u, 0], [0, -alpha, 0, u]])
A
Explanation: Euler Equations with $(\gamma_e)$
We can define $\gamma_e = p/(\rho e) + 1$ and differentiate it to get:
$$
\frac{\partial \gamma_e}{\partial t} = -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x}
$$
The original CG paper used $\tau = 1/\rho$ in place of density.  With this, the continuity equation becomes:
$$
\frac{\partial \tau}{\partial t} = -u\frac{\partial \tau}{\partial x}  + \tau \frac{\partial u}{\partial x}
$$
The Euler equations with this set of primitive variables, $q = (\tau, u, p, \gamma_e)^\intercal$ appear as:
\begin{align}
\frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x}  + \tau \frac{\partial u}{\partial x} \
\frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \
\frac{\partial \gamma_e}{\partial t} &= -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x}
\end{align}
For convenience, we define
$$
\alpha = (\gamma_e - 1)(\gamma_e - \Gamma_1)
$$
and then in vector form, we have:
$$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{cccc} u  & -\tau     & 0    & 0\                        
                                  0  &  u       & \tau  & 0\                   
                                  0  & c^2/\tau & u       & 0\
                                  0  & -\alpha & 0       & u\end{array} \right )
$$
We can represent this matrix symbolically in SymPy and explore its eigensystem.
End of explanation
A.eigenvals()
Explanation: The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
End of explanation
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
    print(e.name)
    display(e.r, e.l)
Explanation: We see that there are 2 eigenvalues $u$ -- the addition of $\gamma_e$ to the system adds this degeneracy.
Eigenvectors
End of explanation
dtau = symbols(r"\Delta\tau")
dge = symbols(r"\Delta\gamma_e")
dq = Matrix([[dtau, du, dp, dge]]).transpose()
taui = symbols(r"\tau_\mathrm{int}")
gei = symbols(r"(\gamma_e)_\mathrm{int}")
qint = Matrix([[taui, ui, pri, gei]]).transpose()
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))
Explanation: $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$$
\Delta q = \left ( \begin{array}{c} \Delta \tau \ \Delta u \ \Delta p \ \Delta \gamma_e \end{array} \right )
$$
and then
$\beta^\nu = l^\nu \cdot \Delta q$
End of explanation
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))
Explanation: and now the final interface states
End of explanation
cg = symbols(r"c_g", real=True, positive=True)
hg = symbols(r"h_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True)
f = symbols("f", real=True, positive=True)
q = Matrix([[rho, u, p, re, Er]]).transpose()
A = Matrix([[u, rho, 0, 0, 0], 
            [0, u, rho**-1, 0, lf/rho], 
            [0, rho*cg**2, u, 0, 0], 
            [0, rho*hg, 0, u, 0],
            [0, (lf+1)*Er, 0, 0, u]])
A
Explanation: Gray FLD Radiation Euler Equations with $(\rho e)$
Following Zhang et al. (2011), the equations of gray radiation hydrodynamics with primitive
variables $q = (\rho, u, p, (\rho e)_g, E_r)^\intercal$ are:
\begin{align}
\frac{\partial \rho}{\partial t} &= -u\frac{\partial \rho}{\partial x}  - \rho \frac{\partial u}{\partial x} \
\frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x}
    - \frac{\lambda_f}{\rho} \frac{\partial E_r}{\partial x}\
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c_g^2 \frac{\partial u}{\partial x} \
\frac{\partial (\rho e)_g}{\partial t} &= -u \frac{\partial (\rho e)_g}{\partial x} - \rho h \frac{\partial u}{\partial x} \
\frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x}
\end{align}
where $(\rho e)_g$ is the gas internal energy density, $h_g$ is the gas specific enthalpy, $c_g$ is the gas sound speed (obeying $c_g^2 = \Gamma_1 p /\rho$), $E_r$ is the radiation energy density, $f$ is the Eddington factor, and $\lambda_f$ is the flux limiter.
Following Zhang et al., we make the approximation that 
$$
\frac{3-f}{2} = \lambda_f + 1
$$
and then in vector form, we have:
$$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{ccccc} u  & \rho     & 0    & 0    & 0 \                        
                                      0  &  u       & 1/\rho & 0   & \lambda_f/\rho\                   
                                      0  & \rho c_g^2 & u       & 0 & 0 \
                                      0  & \rho h_g   & 0       & u & 0 \
                                      0  & (3-f)E_r/2   & 0       & 0 & u\end{array} \right )
$$
We can represent this matrix symbolically in SymPy and explore its eigensystem.
End of explanation
A.eigenvals()
Explanation: The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
End of explanation
cc = c**2 - (lf +1)*lf*Er/rho
evs = A.eigenvals()
for e in evs.keys():
    display(powsimp(simplify(e.subs(cg**2, cc))))
Explanation: We see that there are 3 eigenvalues $u$.  We identify the total sound speed (radiation + gas) as:
$$c^2 = c_g^2 + (\lambda_f + 1)\frac{\lambda_f E_r}{\rho}$$
We can simplify these by substituting in that relationship
End of explanation
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))
# printing them out for inspection
for e in eigen:
    print(e.name)
    display(e.r)
    display(simplify(e.l))
Explanation: We see that these are the same form of the eigenvalues for the pure hydrodynamics system
Eigenvectors
End of explanation
dEr = symbols(r"\Delta{}E_r")
Eri = symbols(r"{E_r}_\mathrm{int}")
dq = Matrix([[dr, du, dp, dre, dEr]]).transpose()
qint = Matrix([[rhoi, ui, pri, rei, Eri]])
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))
Explanation: $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$\beta^\nu = l^\nu \cdot \Delta q$
End of explanation
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))
Explanation: and now the final interface states
End of explanation
geg = symbols(r"{\gamma_e}_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True, positive=True)
cg = symbols(r"c_g", real=True, positive=True)
q = Matrix([[tau, u, p, geg, Er]]).transpose()
A = Matrix([[u, -tau, 0, 0, 0], 
            [0, u, tau, 0, tau*lf], 
            [0, cg**2/tau, u, 0, 0], 
            [0, -alpha, 0, u, 0],
            [0, (lf+1)*Er, 0, 0, u]])
A
Explanation: Gray FLD Radiation Euler Equations with $(\gamma_e)$
We now look at the same system with a different auxillary thermodynamic variable (as we did with pure hydro), using $q = (\tau, u, p, {\gamma_e}_g, E_r)^\intercal$:
\begin{align}
\frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x}  + \tau \frac{\partial u}{\partial x} \
\frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \tau \frac{\partial p}{\partial x}
    - \tau \lambda_f \frac{\partial E_r}{\partial x}\
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \frac{c_g^2}{\tau} \frac{\partial u}{\partial x} \
\frac{\partial {\gamma_e}_g}{\partial t} &= -u \frac{\partial {\gamma_e}_g}{\partial x} + \alpha \frac{\partial u}{\partial x} \
\frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x}
\end{align}
here, ${\gamma_e}_g$ is defined solely in terms of the gas pressure and energy, the remaining variables have the same meaning as above.  We again make the approximation that 
$$
\frac{3-f}{2} = \lambda_f + 1
$$
and then in vector form, we have:
$$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{ccccc} u  & -\tau    & 0    & 0    & 0 \                        
                                      0  &  u       & \tau & 0   & \tau \lambda_f\                   
                                      0  & c_g^2/\tau & u       & 0 & 0 \
                                      0  & -\alpha   & 0       & u & 0 \
                                      0  & (3-f)E_r/2   & 0       & 0 & u\end{array} \right )
$$
We can represent this matrix symbolically in SymPy and explore its eigensystem.
End of explanation
A.eigenvals()
Explanation: The eigenvalues are the speeds at which information propagates with.  SymPy returns them as a
dictionary, giving the multiplicity for each eigenvalue.
End of explanation
cc = c**2 - (lf +1)*lf*Er*tau
evs = A.eigenvals()
for e in evs.keys():
    display(powsimp(simplify(e.subs(cg**2, cc))))
Explanation: We see that there are 3 eigenvalues $u$.  We identify the total sound speed (radiation + gas) as:
$$c^2 = c_g^2 + (\lambda_f + 1)\lambda_f E_r\tau$$
We can simplify these by substituting in that relationship
End of explanation
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))
# printing them out for inspection
for e in eigen:
    print(e.name)
    display(e.r)
    display(simplify(e.l))
Explanation: We see that these are the same form of the eigenvalues for the pure hydrodynamics system
Eigenvectors
End of explanation
dgeg = symbols(r"\Delta{\gamma_e}_g")
dEr = symbols(r"\Delta{}E_r")
dq = Matrix([[dtau, du, dp, dgeg, dEr]]).transpose()
Eri = symbols(r"{E_r}_\mathrm{int}")
qint = Matrix([[taui, ui, pri, gei, Eri]])
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
    betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
    print(e.name)
    beta = e.l.dot(dq)
    display(Eq(betas[n], simplify(beta)))
Explanation: $\beta$'s and final update
The tracing projects the primitive variables into characteristic variables by defining
$\beta^\nu = l^\nu \cdot \Delta q$
End of explanation
for n in range(len(eigen)):
    rhs = q[n]
    for m in range(len(eigen)):
        rhs -= betas[m]*eigen[m].r[n]
    display(Eq(qint[n],rhs))
Explanation: and now the final interface states
End of explanation
<END_TASK>