import numpy as np from numpy.polynomial.legendre import leggauss exec(open('/workspace/rhpc/h2/dirichlet/opt.py').read().split('for A in')[0]) for R in [2.0, 2+50/1093, 2.1, 2.25, 2.5, 3.0, 4.0]: v,u,h=m(R,600); xg,w=leggauss(600); uu=xg*R/2; ww=w*R/2 cp=[v/abs((h*ww*np.exp(2j*np.pi*e*uu)).sum())**2 for e in (0.05,0.1)] print(f"R={R:.4f} m={v:.6f} simple>={1-v:.4%} close(0.05)={cp[0]:.4f} close(0.1)={cp[1]:.4f}")