fig = figure(figsize=[12,3])
# Add subplots: number in y, x, index number
ax = fig.add_subplot(121,autoscale_on=False,xlim=(-10,10),ylim=(-2,2))
ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(-10,10),ylim=(-2,2))
# Set up limits on the graph: we will plot from -10 to 10, with the boundary at 0
xmin = -10
xmax = 10
# Create x arrays for x<0 and x>0
xL = linspace(xmin,0.0,1000)
xR = linspace(0.0,xmax,1000)
# System parameters: tension, masses per unit length, frequency, amplitude of incoming wave
T = 1.0
mu1 = 1.0
mu2 = 2.0
omega = 1.0
A = 1.0
# The time at which we visualise can have an important effect, and produce surprising results
time = 0.0
# Derive speeds, wavenumbers, impedances
c1 = sqrt(T/mu1)
c2 = sqrt(T/mu2)
k1 = omega/c1
k2 = omega/c2
Z1 = sqrt(T*mu1)
Z2 = sqrt(T*mu2)
# Ra and Ta are reflection and transmission coefficients (to avoid confusion with tension T above)
Ra = (Z1-Z2)/(Z1+Z2)
Ta = 1+Ra
print "Z1 and Z2 are: ",Z1,Z2
print "R and T are: ",Ra,Ta
# Create psi_i, psi_r and psi_t (include time - note that psi_r travels in negative x direction)
left_i = A * cos(k1*xL-omega*time)
left_r = A * Ra*cos(k1*xL+omega*time)
# left is the total wave for x<0
left = left_i + left_r
# right is the wave for x>0
right = A*Ta*cos(k2*xR-omega*time)
# Plot total waves in left plot
ax.plot(xL,left,label=r"$\psi_i+ \psi_r$")
ax.plot(xR,right,label=r"$\psi_t$")
ax.legend()
# Plot individual waves in right plot
ax2.plot(xL,left_i,label=r"$\psi_i$")
ax2.plot(xR,right,label=r"$\psi_t$")
ax2.plot(xL,left_r,label=r"$\psi_r$")
ax2.legend()

fig = figure(figsize=[12,7])
# Set length scale
L = 5.0
# We will put the subplots into an array, ax
ax = []
# Define x coordinates and omega
x = linspace(0,L,1000)
omega = 1.0
# Loop over first five modes, n
for n in range(1,6):
    # Create subplot; we want 5 in y and 1 in x, with incremented index - spn gives this
    spn = 510+n
    # Plot mode n
    ax.append(fig.add_subplot(spn,autoscale_on=False,xlim=(0,L),ylim=(-1.1,1.1)))
    # Plot at a number of times to show how the wave changes
    for t in range(4):
        ax[n-1].plot(x,sin(n*pi*x/L)*sin(omega*pi*t/2.0))

fig = figure(figsize=[12,7])
# Set length scale
L = 5.0
# We will put the subplots into an array, ax
ax = []
# Define x coordinates and omega
x = linspace(0,L,1000)
omega = 1.0
# Loop over first five modes, n
for n in range(1,6):
    # Create subplot; we want 5 in y and 1 in x, with incremented index - spn gives this
    spn = 510+n
    # Plot mode n
    ax.append(fig.add_subplot(spn,autoscale_on=False,xlim=(0,L),ylim=(-1.1,1.1)))
    # Plot at a number of times to show how the wave changes
    for t in range(4):
        ax[n-1].plot(x,cos(n*pi*x/L)*sin(omega*pi*t/2.0))