% spectral envelop
rhohalf = envlp(R,A,B,th);
rho = rhohalf.^2;
plot(th,rho);