#!/usr/bin/env python
# coding: utf-8

# # Chebyshev polynomials on the unit circle
# 
# This notebook illustrates the relation between Fourier modes on the unit circle and Chebyshev polynomials.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'notebook')


# In[2]:


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import *


# In[3]:


k = 15  # wave number.  Try 1, 2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
theta = linspace(0, 2*pi, 1000)
x = cos(theta)
y = sin(theta)
z = cos(k*theta)
ax.plot(x,y,z,'b')
ax.plot(x,y,0*z,'k')  # plot unit circle

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.draw()


# If you run this notebook you should be able to rotate the figure above.  If you view it in the $x$-$z$ plane, you get the figure below.
# 
# ## Chebyshev polynomial
# 
# Plotting $z = \cos(k \theta)$ as a function of $x = \cos(\theta)$ gives the Chebyshev polynomial $T_k(x)$:

# In[4]:


plt.figure(figsize=(10,5))
plt.plot(x,z)
plt.plot([-1,1],[0,0],'k')


# # Aliasing
# 
# Because of the relation between Chebyshev polynomials and Fourier modes on the unit circle, aliasing takes place as described in ATAP Theorem 4.1 and illustrated below.

# In[25]:


n = 7   # Use Chebyshev grid with n+1 points
         
k = 17  # Start with T_k
m = abs(mod(k+n-1,2*n) - (n-1))  # aliased to T_m

# Fine grid for plotting polynomials:
theta = linspace(0, pi, 500)
x = cos(theta)
Tk = cos(k*theta)
Tm = cos(m*theta)
         
# Grid points:
theta_j = linspace(0, pi, n+1)
xj = cos(theta_j)
Tk_j = cos(k*theta_j)

plt.figure()
plt.plot(x,Tk,'b')
plt.plot(xj, Tk_j, 'bo')
plt.plot(x,Tm,'r')

plt.title("The Chebyshev polynomial T_{%s} gets aliased to T_{%s}" % (k,m))


# ### On the unit circle
# 
# Note that each Chebyshev point in the interval (with the exception of the endpoints) maps to two points on the unit circle.  Below only one of these points is plotted.

# In[26]:


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
theta = linspace(0, 2*pi, 1000)
x = cos(theta)
y = sin(theta)
Tk = cos(k*theta)
Tm = cos(m*theta)
ax.plot(x,y,Tm,'b')
ax.plot(x,y,Tk,'r')
ax.plot(x,y,0*Tm,'k')  # plot unit circle

theta_j = linspace(0, pi, n+1)
xj = cos(theta_j)
yj = sin(theta_j)
Tk_j = cos(k*theta_j)
ax.plot(xj, yj, Tk_j, 'bo')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.draw()


# In[ ]: