from numpy import *
import matplotlib.pyplot as plt


om0 = 1.
max_freq = 2

om = linspace(0,max_freq*om0,500)

def imampl(om,gam):
    return 2.*gam*om/((om0**2 - om**2)**2 + (2*gam*om)**2)

def reampl(om,gam):
    return (om0**2 - om**2)/((om0**2 - om**2)**2 + (2*gam*om)**2)

#imamplitude = imampl(om,0.1)
#reamplitude = reampl(om,0.1)

#maxampl = amax(amplitude)

fig, (ax1,ax2) = plt.subplots(2,sharex=True)
for gam, c in zip([0.05,0.1, 0.25, 0.5],['k','b','r','m']):
    amplitude = imampl(om,gam)
    ax1.plot(om,amplitude,lw=2,color=c,label=r'$\gamma = {}$'.format(gam))
    ind = argmax(amplitude)
    ax1.plot([om[ind],om[ind]],[0,amax(amplitude)],'--',color=c,lw=2)
    ax2.plot(om,reampl(om,gam),lw=2,color=c)
ax2.set_xlabel(r'normierte Frequenz $\omega/\omega_0$',fontsize=16)
ax2.plot([0,max_freq*om0],[0,0],'--',lw=2,c='k')
ax1.set_ylim([0.05,12])
ax1.set_ylabel('Im(A)',fontsize=16)
ax2.set_ylabel('Re(A)',fontsize=16)
ax1.tick_params(labelsize=15)
ax2.tick_params(labelsize=15)
ax1.legend(loc='upper right')
for ax in (ax1,ax2):
    ax.label_outer()
plt.subplots_adjust(hspace=0)
plt.savefig('Im_Amplitude.eps',dpi=100)
plt.show()