function [Average_AoI, Average_PeakAoI, Average_Energy_Consumption, Average_AoI_Analytical, Average_PeakAoI_Analytical, Average_Energy_Consumption_Analytical]=main(P,v,Omega,noi,p,L)
%Average_AoI Average_PeakAoI Average_Energy_Consumption are the simulation
%results, and Average_AoI_Analytical Average_PeakAoI_Analytical
%Average_Energy_Consumption_Analytical are the analytical results

% P is the transmit power in watts
% v = 2^R-1, and R is the transmission rate in bits/Hz, for R = 1, v = 1
% Omega is the channel power gain
% noi is the AWGN noise power at the receiver
% p is the generation probability of status updates
% L is the maximum retransmission times
q=1-exp(-noi.*v./P./Omega); % calculate outage probability of each transmission by Eq. (2) 

ran=10^8; %number of simulated time slots
retran_number=0; %retransmission number, 0 indicates that no stored status update at source
aoi=ones(1,ran+1); %AoI vector for ran time slots with additional initial time slot
sum_peakaoi=0;  % sum of peak AoI
count_peakaoi=0; % count of peak AoI used to calculate average peak AoI
E=0; %energy consumption
for i=1:ran
%Set up random values for the indicators of the generation of status updates
%,and outage events between source and destination
    indicator_generation=rand;    
    indicator_outage=rand;
%New packet generates at source and transmission fails 
    if (indicator_generation<=p)&&(indicator_outage<=q)
        aoi(i+1)=aoi(i)+1;
        retran_number=1;
        E=E+P;
%New packet generates at source and transmission succeeds
    elseif (indicator_generation<=p)&&(indicator_outage>q)
        aoi(i+1)=1;
        retran_number=1;
        E=E+P;
        %Note that destination may receive duplicated status updates which
        %does not count in peak AoIs
        if (aoi(i+1) < (aoi(i)+1))
        sum_peakaoi=sum_peakaoi+aoi(i);
        count_peakaoi=count_peakaoi+1;
        end
%No packet generates at source and no stored status update at source
    elseif (indicator_generation>p)&&(retran_number==0)
        aoi(i+1)=aoi(i)+1;
        retran_number=0;
%No packet generates at source and retransmission number reaches L
    elseif (indicator_generation>p)&&(retran_number==L)
        aoi(i+1)=aoi(i)+1;
        retran_number=0; %discard the stored status updates
%No packet generates at source, there exists a status update stored at the
%source with a retransmission number smaller than L, outage event happens
    elseif (indicator_generation>p)&&(retran_number~=0)&&(retran_number~=L)&&(indicator_outage<=q) 
        aoi(i+1)=aoi(i)+1;
        retran_number=retran_number+1;
        E=E+P;
%No packet generates at source, there exists a status update stored at the
%source with a retransmission number smaller than L, transmission succeeds
    else
        aoi(i+1)=retran_number+1;
        retran_number=retran_number+1;
        E=E+P;
        %Note that destination may receive duplicated status updates which
        %does not count in peak AoIs
        if (aoi(i+1) < (aoi(i)+1))
        sum_peakaoi=sum_peakaoi+aoi(i);
        count_peakaoi=count_peakaoi+1;
        end
    end
end

Average_AoI=sum(aoi)./ran;
Average_PeakAoI=sum_peakaoi/count_peakaoi;
Average_Energy_Consumption=E./ran;
Average_AoI_Analytical=(1-q+p*q)./(p-p*q)/(1-(q-p*q)^L); %Eq. (33)
Average_PeakAoI_Analytical=1/(1-q+p*q)-L*(q-p*q)^L/(1-(q-p*q)^L)+(1-q+p*q)/(p-p*q)/(1-(q-p*q)^L)-1; %Eq. (36)
Average_Energy_Consumption_Analytical=P*(1-(1-p)^L); %Eq. (45)