#!/usr/bin/env python
# coding: utf-8

# ## Exercise for the course [Python for MATLAB users](http://sese.nu/python-for-matlab-users-ht15/)

# Original exercise by Claus Führer, modified by Olivier Verdier

# In[ ]:


get_ipython().run_line_magic('pylab', '')
get_ipython().run_line_magic('matplotlib', 'inline')


# -----

# ## System matrix

# Consider the matrix
# \\[
#  A=\begin{bmatrix}
#   0 & I \\
# K & D
#  \end{bmatrix}
# \\]
# where $0$ and $I$ are the $2 \times 2$ zero and identity matrices and $K$ and $D$ are $2 \times 2$
# matrices of the following form:
# \\[
#  K=\begin{bmatrix}
#     -k & 0.5 \\ 0.5 & -k
#    \end{bmatrix}
# \qquad
# D=\begin{bmatrix}
#     -d & 1.0 \\ 1.0 & -d
#    \end{bmatrix}
# \\]
# with $k$ and $d$ being real parameters.
# 
# 

# ### Write a function `stiffness` which constructs the matrix $K$ above.

# In[ ]:


def stiffness(k):
    return zeros([2,2]) # implement this!


# In[ ]:


assert(allclose(stiffness(1.), array([[-1.,.5],[.5,-1.]])))


# ### Write a function `damping` which constructs the matrix $D$ above.

# In[ ]:


def damping(d):
    return zeros([2,2]) # implement this!


# In[ ]:


assert(allclose(damping(1.), array([[-1.,1.],[1.,-1.]])))


# ### Write a function `system_matrix` which takes $k$ and $d$ as input and which generates the matrix $A$.

# 
# Hint: use the function `concatenate`. Check its documentation by running:

# In[ ]:


get_ipython().run_line_magic('pinfo', 'concatenate')


# Use also `identity` (or `eye`), `zeros` (or `zeros_like`).

# In[ ]:


def system_matrix(d, k):
    return zeros([4,4]) # implement this!


# In[ ]:


A = system_matrix(10.,20.)
assert(allclose(A[:2,:2], zeros([2,2])))
assert(allclose(A[:2,2:4], identity(2)))
assert(allclose(A[2:4,:2], stiffness(20.)))
assert(allclose(A[2:4,2:4], damping(10.)))


# ### For samples of the values $d \in [0,100]$ and the fixed value $k=1000$, plot the four eigenvalues on the complex plane.

# In[ ]:


for d in linspace(0,100,200):
    pass


# ### Bonus question: there is a bifurcation in the diagram above. Can you find the bifurcation point programmatically?

# In[ ]:





# ## Frequency Response Plot

# In technical applications there occurs often linear systems of the form
# \\[
# \dot x(t) = A x(t) + B u(t)
# \\]
# where $u$ is an given input signal. $x$ is called the state. From the state some quantities $y(t)$ can be
# measured, this is decribed by the equation
# \\[
#  y(t)=C x(t).
# \\]
# We assume here that the input signal is an harmonic oscillation $u(t)=\sin(\omega t)$ with a given frequency $\omega$ and amplitude one. Then, $y(t)$ is again a harmonic oscillation with the same frequency, but another amplitude. The amplitude depends on the frequency. 
# 
# The relationship between the in- and out-amplitude is given by the formula
# \\[
#  \mathrm{amplitude}(\omega)=\\|(G(i\omega))\\|\quad\text{where}
# \quad G(i\omega)=C \cdot (i\omega I -A)^{-1} \cdot B
# \\]
# and $i$ is the imaginary unit.
# 
# 

# In[ ]:


get_ipython().run_line_magic('pinfo', 'inv')


# In[ ]:


get_ipython().run_line_magic('pinfo', 'norm')


# In[ ]:


def amplitude(A, B, C, omega):
    pass


# Plot the amplitude versus omega, for $\omega \in [0, 160]$, with $A$ being the system matrix above with $d=20$ and $k=500$, and
# \\[
#  C=\begin{bmatrix} 1 & 0 & 0 & 0  \end{bmatrix} \qquad B=\begin{bmatrix}0 \\ 0 \\ 0\\ 1 \end{bmatrix} .
# \\]
# 

# In[ ]:





# Find out the relationship between $A$'s eigenvalues and the peak(s) in the figure.
# 

# In[ ]: