%% Simplying the integral
% Let's be a bit more clever and see if we can integrate out one of the
% variables in closed form first so that the numerical integral is not so
% computationally tough.


syms N0 r phi; % Use matlab symbolic
f_rxr = r*exp(-r^2/N0)*exp(2*r*cos(phi)/N0); % now it's a symbolic function

int(f_rxr, phi, 0, pi/8) % doesn't help at all

int(f_rxr,r,0,Inf); % this gives a limit that we can evaluate
pretty(ans)

% after some manipulation by looking at the limits we can actually get the
% function that we have to integrate over phi as:
% int(f_rxr, r, 0, inf) = exp(-1/N0)*(pi^(1/2)*exp(cos(phi)^2/N0)*cos(phi)*(1/N0)^(1/2) + pi^(1/2)*exp(cos(phi)^2/N0)*erf(cos(phi)*(1/N0)^(1/2))*cos(phi)*(1/N0)^(1/2) + 1)/pi


% Now we can just do the numerical integral over phi of the previous
% expression
clear N0 r phi

%% Simplified integral
%%%%%%%%%%% PLAY WITH THESE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Eb_N0dB = 10; % Eb/N0 in dB scale
N_phi = 11960001; % number of points in phi that we are evaluating for the integral
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Eb_N0 = 10^(Eb_N0dB/10);
N0=1/Eb_N0;
I=0;
 
phi_range=linspace(0,pi/8,N_phi);
phi_range=phi_range(1:end-1);
ste_phi=phi_range(2)-phi_range(1);
I=0;
tic % Measure time for computing integral
for phi=phi_range
    I=I+ste_phi*exp(-1/N0)*(pi^(1/2)*exp(cos(phi)^2/N0)*cos(phi)*(1/N0)^(1/2) + pi^(1/2)*exp(cos(phi)^2/N0)*erf(cos(phi)*(1/N0)^(1/2))*cos(phi)*(1/N0)^(1/2) + 1)/pi;
end
toc %
Pe=1-I