from __future__ import division

t = linspace(-5,5,300) # redefine this here for convenience

fig,ax = subplots()

fs=5.0
ax.plot(t,sinc(fs * t))
ax.grid()
ax.annotate('This keeps going...',
            xy=(-4,0),
            xytext=(-5+.1,0.5),
            arrowprops={'facecolor':'green','shrink':0.05},fontsize=14)
ax.annotate('... and going...',
            xy=(4,0),
            xytext=(3+.1,0.5),
            arrowprops={'facecolor':'green','shrink':0.05},fontsize=14)

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)

def kernel(x,sigma=1):
    'convenient function to compute kernel of eigenvalue problem'
    x = np.asanyarray(x)
    y = pi*where(x == 0,1.0e-20, x)
    return sin(sigma/2*y)/y

nstep=100                # quick and dirty integral quantization
t = linspace(-1,1,nstep) # quantization of time
dt = diff(t)[0]          # differential step size
def eigv(sigma):
    return eigvalsh(kernel(t-t[:,None],sigma)).max() # compute max eigenvalue

sigma = linspace(0.01,4,15) # range of time-bandwidth products to consider

fig,ax = subplots()
ax.plot(sigma, dt*array([eigv(i) for i in sigma]),'-o')
ax.set_xlabel('time-bandwidth product $\sigma$',fontsize=14)
ax.set_ylabel('max eigenvalue',fontsize=14)
ax.axis(ymax=1.01)
ax.grid()

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)

sigma=3
w,v=eigh(kernel(t-t[:,None],sigma))
maxv=v[:, w.argmax()]
fig,ax=subplots()
ax.plot(t,maxv)
ax.set_xlabel('time',fontsize=14)
ax.set_title('Eigenvector corresponding to e-value=%2.4e;$\sigma$=%3.2f'%(w.max()*dt,sigma))

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)

import scipy.special

# quick definition to clean up the function signature and normalize
def pro(t,k=0,sigma=3): 
    'normalized prolate angular spherioidal wave function wrapper'
    tmp= scipy.special.pro_ang1(0,k,pi/2*sigma,t)[0] #kth prolate function
    den=linalg.norm(tmp[np.logical_not( np.isnan(tmp))])# drop those pesky NaNs at edges
    return tmp/den # since e-vectors are likewise normalized

fig,ax=subplots()

ax.plot(t,pro(t),'b.',t,maxv,'-g');
ax.set_xlabel('time')
ax.set_title('Comparison of $0^{th}$ prolate function')


# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)

W=1/2 # take bandwidth = 1 Hz. This makes the total signal power=1
# compute epsilon for this time-extent
epsilon=1 - 2/pi*scipy.special.sici( 2*pi )[0] # using sin-integral special function

fig,ax = subplots()

t = linspace(-3,3,100)
ax.plot(t,sinc(t)**2)
ax.hlines(0,-3,3)
tt = linspace(-1,1,20)
ax.fill_between(tt,0,sinc(tt)**2,facecolor='g',alpha=0.5)
ax.set_title('$x(t)^2$: $\epsilon=%2.4f$ over shaded region' % epsilon)
ax.set_xlabel('time')

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)

epsilon= 2/3. + 4/3/pi*scipy.special.sici(2*pi)[0]-8/3/pi*scipy.special.sici(4*pi)[0] # total energy is 2/3 in this case

fig,ax = subplots()
t = linspace(-3,3,100)
ax.plot(t,sinc(t)**4)
ax.hlines(0,-3,3)
tt = linspace(-1,1,20)
ax.fill_between(tt,0,sinc(tt)**4,facecolor='g',alpha=0.5)
ax.set_title('$x(t)^2$: $\epsilon=%2.4f$ over shaded region' % epsilon)
ax.set_xlabel('time')

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)


def kernel_tau(x,W=1):
    'convenient function to compute kernel of eigenvalue problem'
    x = np.asanyarray(x)
    y = pi*where(x == 0,1.0e-20, x)
    return sin(2*W*y)/y

nstep=300                # quick and dirty integral quantization
t = linspace(-1,1,nstep) # quantization of time
tt = linspace(-2,2,nstep)# extend interval

w,v=eig(kernel_tau(t-tt[:,None],5))
ii = argsort(w.real) 
maxv=v[:, w.real.argmax()].real
fig,ax = subplots()
ax.plot(tt,maxv)
##plot(tt,v[:,ii[-2]].real)
ax.set_xlabel('time',fontsize=14)
ax.set_title('$\phi_{max}(t),\sigma=10*2=20$',fontsize=16)

# fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)