function f = PE_Sphere(N)
% matlab script to carry out the graviatational potential energy calculation
% for a sphere. 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.
%
% this function is based on an adaptation of MI_Cone.m
% but the program is simplified to measure the deviation
% by sampling the square of the measurement function.

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

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

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


% Now let us measure the volume:
Nhit=0; 
for nu=1:N 
    if X(nu)^2+Y(nu)^2+Z(nu)^2 <= r0^2 
        Nhit=Nhit+1; 
    end 
end
Nhit=Nhit/N;
disp('Exact answer for the sphere volume:');
V_sphere=4*pi/3*r0^3
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*RefVolume
disp('with unknown deviation');
% Now estimate the potential energy by repeating the above loop, but
% summing the distance from the rotation axis instead of just
% updating the hit counter.
% We define a function at the bottom which yields the 
% potential energy between a probe charge on the x-axis
% at location xp, and a test mass inside the sphere of radius r0
% at location (x,y,z).
% Location of the probe mass:
xp=10;

Nhit=0; 
Nhit2=0;
for nu=1:N 
    if X(nu)^2+Y(nu)^2+Z(nu)^2 <= r0^2 
        Nhit=Nhit+PE(xp,X(nu),Y(nu),Z(nu)); 
        Nhit2=Nhit2+PE(xp,X(nu),Y(nu),Z(nu))^2; 
    end 
end
Nhit=Nhit/N;
Nhit2=Nhit2/N;
dev=sqrt((Nhit2-Nhit^2)/N);

disp('MC estimate of potential energy:');
MC_Estimate=Nhit*RefVolume*(1/V_sphere)
disp('with deviation');
Dev_for_Estimate=dev*RefVolume*(1/V_sphere)
R(1)=MC_Estimate-Dev_for_Estimate;
R(2)=MC_Estimate;
R(3)=MC_Estimate+Dev_for_Estimate;
R


f=R;

function sf1 = PE(xp,x,y,z)
sf1 = 1/sqrt((x-xp)^2+y^2+z^2);