%matplotlib inline
from nsm import *
plt.figure(figsize=(5,10))
k = 4
x = np.linspace(0,2*np.pi,500,endpoint=False)
w = np.exp(1j*x)
for j in range(-k,k+1):
    plt.subplot(2*k+1,1,j+k+1)
    plt.plot(np.real(w**j),'r')
    plt.plot(np.imag(w**j),'b')
    plt.title(f'j={j}')

%matplotlib notebook
ax = plt.figure().add_subplot(projection='3d')
k = 4
x = np.linspace(0,2*np.pi,500,endpoint=False)
w = np.exp(1j*x)
for j in [k]:
    plt.plot(x,np.real(w**j),np.imag(w**j))

%matplotlib inline
twopi = 2*np.pi
i = 1j
def f1(x): return np.exp(np.cos(x)+np.sin(x)/2)  # a smooth periodic function that is not a combination of trig functions
def f2(x): return x**2  # a function with discontinuous periodic extension
def f3(x): return x*(twopi-x)  # a function with continuous periodic extension but jump in derivative
def f4(x): x1 = (x-np.pi)/np.pi; return x1**4 - 2*x1**2 

f = f2

x = np.linspace(0,twopi,1000)
kstop = 6
plt.figure(figsize=(8,16))
for k in range(kstop):
    plt.subplot(kstop,1,k+1)
    plt.plot(x,f(x),'k',label='f');
    gk = np.zeros_like(x,dtype=complex)
    alpha = []
    for j in range(-k,k+1):
        wj = np.exp(i*j*x)
        def integrand(x): return f(x)*np.exp(-i*j*x)
        twopiaj,err = quadrature(integrand,0,twopi)   # gaussian quadrature
        alphaj = twopiaj/twopi
        alpha.append(alphaj)
        gk += alphaj*wj
    plt.plot(x,np.real(gk),label=f'$Re(g_{k})$',lw=3,alpha=0.5)
    plt.plot(x,np.imag(gk),label=f'$Im(g_{k})$',lw=3,alpha=0.5)
    plt.title(f'k = {k}')
    #print(np.round(a,4))
    plt.legend()

f=f4
x = np.linspace(0,twopi,1000)
kstop = 10
#plt.figure(figsize=(8,16))
kk = np.linspace(1,kstop,3)
for k in range(1,kstop):
    gk = np.zeros_like(x,dtype=complex)
    alpha = []
    for j in range(-k,k+1):
        wj = np.exp(i*j*x)
        def integrand(x): return f(x)*np.exp(-i*j*x)
        twopiaj,err = quadrature(integrand,0,twopi)   # gaussian quadrature
        alphaj = twopiaj/twopi
        alpha.append(alphaj)
        gk += alphaj*wj
    
    maxerror = np.abs(f(x)-gk).max()
    print(k,maxerror)
    if k==1:
        for p in range(1,5):
            plt.loglog(kk,kk**(-p)*maxerror,label='$k^{-'+str(p)+'}$',alpha=0.3)

    plt.loglog(k,maxerror,'ro')
plt.xlabel('k'); plt.ylabel('max error')
plt.legend();

from nsm import *
n = 8
x = np.linspace(0,2*np.pi,500,endpoint=False)
w = np.exp(1j*x)
for j in range(n):
    plt.subplot(n,1,j+1)
    plt.plot(np.real(w**j),'r')
    plt.plot(np.imag(w**j),'b')

np.set_printoptions(linewidth=400)
n = 8
w = np.exp(-1j*2*np.pi/n)
A = np.empty((n,n),dtype=complex)
v = [w**j for j in range(n)]
A[:,0] = 1
for j in range(1,n):
    A[:,j] = A[:,j-1]*v
A /= n
plt.imshow(np.real(A));

y = np.random.randn(n)
x = np.linspace(0,2*np.pi,n,endpoint=False)
plt.plot(x,y,'ko',alpha=0.5)
plt.plot(x,y*0,'ko',alpha=0.25)
nn = 1000
xx = np.linspace(0,2*np.pi,nn,endpoint=False)
ww = np.exp(1j*xx)
c = A@y
#print(c)
p = np.zeros(nn,dtype=complex)
for j in range(n):
    print(c[j])
    p += c[j]*ww**j
plt.plot(xx,np.real(p),'r')
plt.plot(xx,np.imag(p),'b')
plt.axhline(0)
plt.xlim(0,2*np.pi)

%matplotlib notebook
ax = plt.figure().add_subplot(projection='3d')
ax.plot( x,y,x*0,'o')
ax.plot(xx,np.real(p),np.imag(p));