#!/usr/bin/env python
from __future__ import division
from numpy import *
import matplotlib.pyplot as plt

m = 1
k = 1
gamma = .1
omega = linspace(-4,4,500)

A = 1/(k-m*(omega**2) + (gamma*omega)*1j)

plt.figure()
plt.plot(omega, angle(A), label='phase')
plt.xlabel('$\omega$ (rad/s)')
plt.ylabel('phase of response (rad)')
plt.grid()
plt.title("Phase of response vs. driving frequency of $m\ddot{x} + \gamma \dot{x} + kx = e^{i\omega t}$")

plt.figure()
plt.plot(omega, absolute(A), label='magnitude')
plt.xlabel('$\omega$ (rad/s)')
plt.ylabel('amplitude of response')
plt.grid()
plt.title("Amplitude of response vs. driving frequency of $m\ddot{x} + \gamma \dot{x} + kx = e^{i\omega t}$")

plt.show()