%%% This program produces makes a plot of C(f) on pp.475 of Proakis 2002
function EPh
clc;clear;close all
altarr = -0.05:-0.021:-3;ustarr = 0.05:0.021:3;tol = 1e-6;
resarr = [];
Eps = quade(@funEhs,-inf,inf,tol);
%Eps = quad(@funEhs,-10,10,tol);
for ilim = 1:length(altarr)
    alt = altarr(ilim);ust = ustarr(ilim);
rest = quade(@funEhs,alt,ust,tol);
%rest = quad(@funEhs,alt,ust,tol);
resarr = [resarr rest];
end
figure(1)
resarr = resarr / Eps;
plot(ustarr,resarr,'-k ','LineWidth',4);
set(gcf,'Color','w');set(gca,'FontSize',20)

xlabel('\bf\itW T','FontSize',24,'FontWeight','bold');%set(gca,'XTickLabel',Tarrstr); 
ylabel('\bf\it\epsilon_h / \epsilon_s','FontSize',30,'FontWeight','bold');
axis ([0 max(ustarr)*1.05 0 1.1]);
grid on;axis on
function y = funEhs(alfa)
y = (sin(pi*alfa)).^2./alfa.^2./(1 - alfa.^2).^2;