function f = MI_cone(N)
% matlab script to carry out the volume and moment-of-inertia calculation
% for a cone. We declare it a function so that we can call it with
% different number of points. Also, this allows us to put locally used
% functions at the bottom.



% first we dial up randoms to generate 3dim position vectors.
% N=10^3; % this is now provided by the function argument!
S1=random('unif',0,1,1,N); 
S2=random('unif',0,1,1,N); 
S3=random('unif',0,1,1,N); 
% to represent a cone we chose a reference volume that is a box of
% base area 2*r0 by 2*r0 by h, where h is the height of the cone.
% For this we remap the randoms as follows:

r0=1
h=3
RefVolume=(2*r0)^2*h;

X=(2*S1-1)*r0;
Y=(2*S2-1)*r0;
Z=S3*h;

% unlike in the maple implementation we have created a new
% set of arrays out of S1, S2, S3.

% We set up a function to determine the radius of the cone as a function
% of chosen height. We could use an inline function. We could define a
% new script M-file. But we can also insert the function that will only
% be used by the current M-file at the bottom of it. Look at the bottom!

% Now let us measure the volume:
for i0=1:10 
    Nhit(i0)=0; 
    for i=1:N/10 
        nu=(i0-1)*(N/10)+i;
        Rh=rad(Z(nu),r0,h);
        if X(nu)^2+Y(nu)^2 <= Rh^2 
            Nhit(i0)=Nhit(i0)+1; 
        end 
    end 
    Nhit(i0)=Nhit(i0)/(N/10); 
end 
Nhit_a=0; 
for i0=1:10 
    Nhit_a = Nhit_a + Nhit(i0); 
end 
Nhit_a=Nhit_a/10; 
%'average number:',Nhit_a
Nhit_d=0; 
for i0=1:10 
    Nhit_d = Nhit_d + (Nhit(i0)-Nhit_a)^2; 
end 
Nhit_d=sqrt(Nhit_d/(9*10)); 
%'with estimated deviation of: ',Nhit_d 
%'while the correct answer is: ',(1/3)*pi*r0^2*h 
disp('Exact answer for the cone volume:');
pi/3*r0^2*h
disp('MC estimate based on the number of points:');
N
disp('is given as:');
% multiply the MC average by the reference volume:
MC_Estimate=Nhit_a*RefVolume
disp('with deviation:');
% multiply the MC average by the reference volume:
Dev=Nhit_d*RefVolume
%f=(Nhit_a-Nhit_d Nhit_a Nhit_a+Nhit_d)*RefVolume;
R(1)=MC_Estimate-Dev;
R(2)=MC_Estimate;
R(3)=MC_Estimate+Dev;
R
% Now estimate the inertia by repeating the above loop, but
% summing the distance from the rotation axis instead of just
% updating the hit counter.

% Use an independent sequence of randoms:
S1=random('unif',0,1,1,N); 
S2=random('unif',0,1,1,N); 
S3=random('unif',0,1,1,N); 
X=(2*S1-1)*r0;
Y=(2*S2-1)*r0;
Z=S3*h;

% Now let us measure the inertia:
for i0=1:10 
    Nhit(i0)=0; 
    for i=1:N/10 
        nu=(i0-1)*(N/10)+i;
        Rh=rad(Z(nu),r0,h);
        if X(nu)^2+Y(nu)^2 <= Rh^2 
            Nhit(i0)=Nhit(i0)+X(nu)^2+Y(nu)^2; 
        end 
    end 
    Nhit(i0)=Nhit(i0)/(N/10); 
end 
Nhit_a=0; 
for i0=1:10 
    Nhit_a = Nhit_a + Nhit(i0); 
end 
Nhit_a=Nhit_a/10; 
%'average number:',Nhit_a
Nhit_d=0; 
for i0=1:10 
    Nhit_d = Nhit_d + (Nhit(i0)-Nhit_a)^2; 
end 
Nhit_d=sqrt(Nhit_d/(9*10)); 
%'with estimated deviation of: ',Nhit_d 
%'while the correct answer is: ',(1/3)*pi*r0^2*h 
disp('Exact answer for the cone inertia for a mass of 1:');
3/10*r0^2
disp('MC estimate based on the number of points:');
N
disp('is given as:');
% multiply the MC average by the reference volume, and multiply with the mass density of the cone:
MC_Estimate=(Nhit_a*RefVolume)*(1/(pi/3*r0^2*h))
disp('with deviation:');
% multiply the MC average by the reference volume:
Dev=Nhit_d*RefVolume*(1/(pi/3*r0^2*h))
%f=(Nhit_a-Nhit_d Nhit_a Nhit_a+Nhit_d)*RefVolume;
R(1)=MC_Estimate-Dev;
R(2)=MC_Estimate;
R(3)=MC_Estimate+Dev;
R


f=R;

function sf1 = rad(z,r0,h)
sf1 = r0-r0/h*z;