n = input('Please input the number of points you would like to throw:')
disp('Now I am trying to do the integration...')

cmd = 'syms pi';
for i=1:n
    cmd = strcat(cmd, ' theta', int2str(i));
end
eval(cmd);

P1_1 = strcat('int(1, theta', int2str(n), ', theta', int2str(n-1), ', theta1+pi)');
for i=n-1:-1:2
     P1_1 = strcat('int(', P1_1, ', theta', int2str(i), ', theta', int2str(i-1), ', theta1+pi)');
end
P1_1 = strcat('int(', P1_1, ', theta1, -pi, 0)')
P1_1 = eval(P1_1)/(2*pi)^n

P1_2 = strcat('int(1, theta', int2str(n), ', theta', int2str(n-1), ', pi)');
for i=n-1:-1:2
    P1_2 = strcat('int(', P1_2, ', theta', int2str(i), ', theta', int2str(i-1), ', pi)');
end
P1_2 = strcat('int(', P1_2 , ', theta1, 0, pi)')
P1_2 = eval(P1_2)/(2*pi)^n

if(n==2)
    P2 = 'int(int(1, theta2, theta1+pi, pi), theta1, -pi, 0)'
else
    P2 = strcat('int(1, theta', int2str(n), ', theta', int2str(n-1), ', pi)');
    for i=n-1:-1:3
        P2 = strcat('int(', P2, ', theta', int2str(i), ', theta', int2str(i-1), ', pi)');
    end
    P2 = strcat('int(', P2, ', theta2, theta1+pi, pi)');
    P2 = strcat('int(', P2, ', theta1, -pi, 0)')    
end
P2 =eval(P2)/(2*pi)^n

P = (P1_1+P1_2+P2)*factorial(n)

clear pi;
disp('Done with the integration. Now I will do the simulation...')
num_sampling = input('Please input the number of sampling:')
num_samples = input('Please input the number of samples in every sampling:')
p = zeros(num_sampling,1);
for j=1:num_sampling
    hit = 0;
    for i=1:num_samples
        angles = rand(n,1)*(2*pi)-pi;
        angles = sort(angles);    
        if( (angles(n)-angles(1)<=pi) || (angles(2)-angles(1)>=pi) )
            hit = hit + 1;
        end
    end
    p(j) = p(j) + hit/num_samples;
end
P_sim = sum(p)/num_sampling