R = 0.05;



delta = linspace(-4*pi,4*pi,3000);

F = 4*R/(1-R)^2;

I = F*sin(delta/2).^2./(1+F*sin(delta/2).^2);



plot(delta,I);