%% Calculates the error probability of 8PSK, where Es=1 (energy per symbol).
 
% This is done by integrating the pink area shown in the last slide of
% Lecture 3. Further, it is sufficient to only integrate what is above the
% x-axis and then multiply the result by 2.
 
% I did it in polar coordinates, so there is an extra "r" from the
% jacobian.
 
% Clearly, this is a problem in numerical analysis and there are many
% tricks that I did not use. As can be seen, the stepsize in the angular
% domain is very small.

Eb_N0dB = 10; % Eb/N0 in dB scale, %%%%%%%%%%%%%  PLAY WITH THIS VARIABLE

Eb_N0 = 10^(Eb_NodB/10);
 
N0=1/Eb_N0; % assume Eb = 1 without loss of generality

%%%%%%%%%%% PLAY WITH THESE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The bigger these numbers are, the more accurate the integral result
N_r = 100; % number of points in r that we are evaluating for the integral
max_r = 100*N0; % maximum r that we are considering (to get the actual result this should be inf)
N_phi = 69501; % number of points in phi that we are evaluating for the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=0;
r_range=linspace(0,max_r,N_r);
ste_r=r_range(2)-r_range(1);
phi_range=linspace(0,pi/8,N_phi);
phi_range=phi_range(1:end-1);
ste_phi=phi_range(2)-phi_range(1);
tic % Measure time for computing integral
for r=r_range
    for phi=phi_range
        
        %x=r*cos(phi)-1;
        %y=r*sin(phi);
        %I=I+2*r*ste_r*ste_phi*exp(-abs(x)^2/N0)*exp(-abs(y)^2/N0)/pi/N0;
        
        % The following line does the same as ther above 3. This follows by simple manipulations of the above 3 lines (easy)
        I=I+2*r*ste_r*ste_phi*exp(-1/N0)*exp(-r^2/N0)*exp(2*r*cos(phi)/N0)/pi/N0;
    end
end
toc %
Pe=1-I