#  Python 2.7 file for Numerical Methods II, Spring 2014
#  http://www.math.nyu.edu/faculty/goodman/teaching/NumericalMethodsII2014/index.html

#  File: NumericalPython.py                                              # bullet 1

#  Illustrate the basics of numerical python ...
#  ... by making a picture of aliasing

#  Type: python aliasing.py   ... and hope for the best

import numpy      as np                                                  # bullet 2
import matplotlib as mp
import pylab      as pl                                                  # bullet 3

print "Hello, aliasing"     # we will drop the "hello" soon

n  = 11           # number of mesh points                                # bullet 4
L  = 4.           # length of the interval                               # bullet 5
dx = L/n          # distance between mesh points

xm = np.zeros(n)  # locations of the mesh points                         # bullet 6
fm = np.zeros(n)  # function values at the mesh points
print "The type of xm is: " + str(type(xm))

k  = 3            # mode number                                          # bullet 7
ka = k + n        # aliased mode number
p  = 2*np.pi*k/L  # the wave number
pa = 2*np.pi*ka/L # the aliased wave number

for i in range(0,n):
  x     =  i*dx
  xm[i] = x
  fm[i] = np.cos( p *x)

m   = 1000        # number of plot points
dxp = L/m         # distance between x values for plotting

xa = np.zeros(m)  #  all x sample points
f1 = np.zeros(m)  #  corresponding values of f with k
f2 = np.zeros(m)  #  corresponding values of f with aliased k


for i in range(0,m):
  x  =  i*dxp
  xa[i] = x
  f1[i] = np.cos( p *x)
  f2[i] = np.cos( pa*x)

#    Plot

fig, ax = pl.subplots()                                                  # bullet 8
l1 = "k = " + str(k)                                                     # bullet 9
l2 = "k = " + str(ka)
ax.plot( xa, f1, 'b',   label = l1)                                      # bullet 10
ax.plot( xa, f2, 'k--', label = l2)                                      # bullet 11
ax.plot( xm, fm, 'ro',  label = 'grid vals')

# Now add the legend and labels
legend = ax.legend(loc='upper right')
ts = 'Two k values give the same grid function, n = ' + str(n) +' points'
pl.title(ts)
pl.xlabel('x')                                                           #bullet 12
pl.ylabel('cos(2*pi*k*x/L')

pl.savefig("aliasing.png")                                               #bullet 13

pl.show()                                                                #bullet 14