% integrazione con il metodo del rigetto

fmax=4/pi;
nmax=10000;
x=rand(nmax,1);
y=rand(nmax,1)*fmax;
f=(4/pi)*sqrt(1-x.^2);
count=0;
ncount=0;
for k=1:nmax
    if y(k)<f(k)
       count=count+1;
       xacc(count)=x(k);
       yacc(count)=y(k);
   else
       ncount=ncount+1;
       xnacc(ncount)=x(k);
       ynacc(ncount)=y(k);

   end
end
int=fmax*count/nmax
xmin=0;
xmax=1;
dx=(1/nmax)*100;
xbin=xmin:dx:xmax;
hcont=histc(xacc,xbin)/(dx*nmax/fmax);
figure(1);
subplot(1,2,1);
f=(4/pi)*sqrt(1-(xbin).^2);
plot(xbin,f,'r')
hold on;
bar(xbin,hcont,'y');
axis([-0.1,1.1,0.,1.401]);
hold off;
subplot(1,2,2);
plot(xacc,yacc,'.')
hold on;
plot(xnacc,ynacc,'.r')
hold off;
