# just some standard imports
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.special import jv as jn
import matplotlib as mpl
mpl.cm.register_cmap(name='cubehelix2', data=mpl._cm.cubehelix(gamma=1.0, s=0.3, r=-0.5, h=1.5))
%matplotlib inline

N = 36
NP = 2*N
k = (np.arange(N)+1.0)
alpha = 2/3.0
t1 = 1.5 * np.pi * (np.arange(NP)+0.5)/NP
r1 = 1./np.maximum( np.abs(np.sin(t1)), np.abs(np.cos(t1)))
t2 = 1.5 * np.pi*sp.rand(NP)
r2 = sp.rand(NP)/np.maximum( np.abs(np.sin(t2)), np.abs(np.cos(t2)) )
t = np.r_[t1,t2]
r = np.r_[r1,r2]

xs = r*np.cos(t)
ys = r*np.sin(t)
plt.scatter(xs,ys);
plt.axis('equal');
plt.axis('off');

def mat(lam):
    return np.sin(alpha*t[:,None]*k[None,:])*jn(alpha*k[None,:], np.sqrt(lam)*r[:,None])
    
def sigma(lam):
    A = mat(lam)
    Q,R = sp.linalg.qr(A, mode='economic')
    s = sp.linalg.svdvals(Q[:NP]).min()
    return s

lamvec = np.r_[0.2:40:0.2]
plt.plot(lamvec, map(sigma, lamvec));
plt.xlabel(r"$\lambda$")
plt.ylabel(r"$\sigma(\lambda)$");

from scipy.optimize import minimize_scalar
lams = [  minimize_scalar(sigma,(z-0.1,z+0.1), tol=1e-16).x for z in [9.6, 15.2, 19.7, 29.5, 31.9, 41.5] ]
print lams

x, y = np.ogrid[-1:1:200j, -1:1:200j]
rr = np.sqrt( x**2 + y**2 )
tt = (sp.arctan2(y,x)+2*np.pi)%(2*np.pi)

fig,axs = plt.subplots(2,3, figsize=(10,5))

for pk,ax in enumerate(axs.flat):
    Q,R = sp.linalg.qr(mat(lams[pk]), mode='economic')
    u,s,vt = sp.linalg.svd( Q[:NP,:], full_matrices=False )
    c = sp.linalg.solve(R, vt[-1])
    ans = (np.sin(alpha*tt[:,:,None]*k[None,None,:])*jn(alpha*k[None,None,:], np.sqrt(lams[pk])*rr[:,:,None])).dot(
         c)
    ans[100:,:100] = np.nan
    #ax.imshow(  ans[:,::-1],cmap='cubehelix2', extent=[-1,1,-1,1]);
    ax.contourf( ans[:,::-1], 10, cmap='cubehelix2')
    ax.set_title(r"$\lambda_{} = {}$".format(pk, lams[pk]))
    ax.axis('off');
    ax.axis('equal');
plt.tight_layout()