from IPython.display import Image
Image(filename='SHMCircularMotion.png',width=300) #Size is in pixels

phi = 30.0*pi/180. # Phase in degrees
rad = 1.0
x = (0,rad*cos(phi))
y = (0,rad*sin(phi))
axis('scaled')
axis([-1.1,1.1,-1.1,1.1])
plot(x,y)
# Add arc to indicate phi
theta = linspace(0,phi,1000)
plot(0.5*rad*cos(theta),0.5*rad*sin(theta),'b-')
# Label
text(0.5, 0.18, r"$\phi$", fontsize=20, color="blue")
# Now draw full circle
theta = linspace(0,2.0*pi,1000)
plot(rad*cos(theta),rad*sin(theta),'g-')
# Plot the phasor
omega = 2.0
t = 2.5
psix = (0,rad*cos(omega*t+phi))
psiy = (0,rad*sin(omega*t+phi))
plot(psix,psiy,'r-')
# Projection on the x-axis
plot(psix,(0.,0.),'r--')
# Label
text(0.4,-0.2, r"$\psi$", fontsize=20, color="red")
# Arc showing omega t
theta = linspace(phi,omega*t+phi,1000)
plot(0.6*rad*cos(theta),0.6*rad*sin(theta),'r:')
# Label
text(-0.5, 0.6, r"$\omega t$", fontsize=20, color="red")

from IPython.display import Image
Image(filename='SHMDispVelAcc.png',width=300) #Size is in pixels

t = linspace(0,10.0,1000)
m = 1.0
k = 2.0
b = 5.0
omega0 = sqrt(k/m)
gamma = b/(2.0*m)
print "omega0 and gamma are ",omega0,gamma
A = 1.0
if gamma<omega0 :
    print "Light damping"
    omega = sqrt(omega0*omega0 - gamma*gamma)
    plot(t,A*exp(-gamma*t)*cos(omega*t))
elif gamma>omega0 : 
    print "Heavy damping"
    mum = gamma + sqrt(gamma*gamma - omega0*omega0)
    mup = gamma - sqrt(gamma*gamma - omega0*omega0)
    B = 1.0
    plot(t,A*exp(-mum*t) + B*exp(-mup*t),label='sum')
    plot(t,A*exp(-mum*t),'r--',label=r"$\mu_{-}$")
    plot(t,B*exp(-mup*t),'g--',label=r"$\mu_{+}$")
    legend()

# Set up parameters of oscillator
m = 1.0
k = 2.0
b = 0.2
gamma = b/(2.0*m)
# Find resonant frequency
omega0 = sqrt(k/m)
print "omega0 and gamma are ",omega0,gamma
F0 = 1.0
# Create array of values of omega from omega0-1 to omega0+1
omega = linspace(omega0-1.0,omega0+1.0,200)
# Calculate phi
phi = arctan(-2.0*gamma*omega/(omega0*omega0 - omega*omega))
# arctan in NumPy is defined to lie between -pi/2 and pi/2 so we correct
# by -pi to give physical reasonableness
phi[100:200] -=pi
# Calculate amplitude
fac = (omega0*omega0 - omega*omega)**2 + 4*gamma*gamma*omega*omega
A = F0/(sqrt(fac)*m)
plot(omega,phi,label='phi')
plot(omega,A,label='A')
legend()

# We use the values of b, m, k, omega and A defined above
power = 0.5*b*omega*omega*A*A
impi = m*omega - k/omega
Z = sqrt(b*b+impi*impi)
plot(omega,power,label='Power')
plot(omega,impi,label='Imaginary part of impedance')
plot(omega,Z,label=r'$\vert Z\vert$')
legend(loc=4)

t = linspace(0,50.0,1000)
# Driving frequency
omegaD = 7.0
# Amplitude of driving (F0) and response (F1)
F0 = 10.0
F1 = 0.10
fac = (omega0*omega0 - omegaD*omegaD)**2 + 4*gamma*gamma*omegaD*omegaD
A = F0/(sqrt(fac)*m)
# Define transient
transient = 2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2)
steady = A*cos(omegaD*t+pi/2)
plot(t,transient+steady)
plot(t,transient)
#plot(t,A*cos(omegaD*t+pi/2)+2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2))
#plot(t,2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2))

from IPython.display import Image
Image(filename='PhasorSameFreqAmp.png',width=200)

from IPython.display import Image
Image(filename='Phasors.png',width=200)

fig = figure(figsize=[12,3])
# Add subplots: number in y, x, index number
ax = fig.add_subplot(121,autoscale_on=False,xlim=(0,65),ylim=(-2,2))
ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(0,65),ylim=(-2,2))
t = linspace(0.0,65.0,1000)
o1 = 1.3
o2 = 1.0
omega = 0.5*(o1+o2)
domega = 0.5*(o1-o2)
# Plot the resulting pattern
ax.plot(t,2.0*cos(omega*t)*cos(domega*t))
# Show the envelope in a dashed red line - we need two lines
ax.plot(t,2.0*cos(domega*t),'r--')
ax.plot(t,2.0*cos(domega*t+pi),'r--')
# On a separate plot, show the two individual oscillations
ax2.plot(t,cos(o1*t))
ax2.plot(t,cos(o2*t))

fig = figure(figsize=[12,3])
# Add subplots: number in y, x, index number
xmax = 15.0
ax = fig.add_subplot(121,autoscale_on=False,xlim=(0,xmax),ylim=(-2,2))
ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(0,xmax),ylim=(-2,2))
t = linspace(0.0,xmax,1000)
o1 = 25./3
o2 = 7.1
omega = 0.5*(o1+o2)
domega = 0.5*(o1-o2)
# Plot the resulting pattern
ax.plot(t,2.0*cos(omega*t)*cos(domega*t))
# Show the envelope in a dashed red line - we need two lines
ax.plot(t,2.0*cos(domega*t),'r--')
ax.plot(t,2.0*cos(domega*t+pi),'r--')
# On a separate plot, show the two individual oscillations
ax2.plot(t,cos(o1*t))
ax2.plot(t,cos(o2*t))

from IPython.display import Image
Image(filename='PhasorsDiffFreq2.png',width=200) #Size is in pixels

# Define figure size to get square
fig = figure(figsize=[5,5])
# The range of t may need to change with angular frequencies
t = linspace(0,10*pi,200)
# Original frequencies: 1.2 and 0.6
o1 = 1.2
o2 = 0.6
# Original amplitudes: 1.5 and 1.3
A1 = 1.5
A2 = 1.3
# Calculate x and y projections of the two oscillations when added
x = A1*cos(o1*t) + A2*cos(o2*t)
y = A1*sin(o1*t) + A2*sin(o2*t)
# Plot path of the tip of the vector
plot(x,y)