clc
clear
M=4;                                                                        % Size of The constellation
sayan=0;                                                                    % Counter
tic
for Power=0:2:12
    Power
    sayan=sayan+1;
    disp(sayan)
    P=10^(Power/10);                                                        % Total power of the source and relay nodes
    Pow=[1 1 1 1 1]*P/6;                                                    % Power of each realy node
    Ps=P/6;                                                                 % Power of the source node
    d=[0.0117, 0.1365,0.2844, 0.4692, 0.8938 ];                             % Location of the relay nodes   
    omsr=d.^(-3);                                                           % Variance of the channel coefficients
    omrd=(1-d).^(-3);
    omsd=1;
    L=length(omsr);
    beta(1)=1;                                                              % SEP of the source to relay nodes
    for k=1:L
    cp=Ps*omsr(k)*sin(pi/M)^2;
    beta(k+1)=1-sqrt(cp/(cp+1))*(atan(cot((M-1)/M*pi)*sqrt(cp/(cp+1)))/pi-.5) -(-atan(cot((M-1)/M*pi))/pi+.5);
    end
    mean=[1,omrd];
    Powclosed=[Ps Pow];
    %%%---- This part is for Calculating the closed form Expression of the
    %%%SEP in (24).
    S=0;
    for k=1:L+1
        S=S+closed_eval（M,Powclosed,mean,beta,k);
    end
    S=S-(atan(cot((M-1)/M*pi))/pi-.5);
    %%%---- This part is for evaluation of the SEP using the Monte-Carlo simulation
    iter=10^4;                      % Number of the independent runs
    counter=0;
    X=[1+j 1-j -1+j -1-j]/sqrt(2);
    for k=1:iter
        k
       symb=floor(rand(1,1)*4)+1;
        NR=(wgn(L,1,0)+j*wgn(L,1,0))/sqrt(2);
        ND=(wgn(L,1,0)+j*wgn(L,1,0))/sqrt(2);
        NSD=(wgn(1,1,0)+j*wgn(1,1,0))/sqrt(2);
        HSR=sqrt(omsr').*(wgn(L,1,0)+j*wgn(L,1,0))/sqrt(2);
        HRD=sqrt(omrd').*(wgn(L,1,0)+j*wgn(L,1,0))/sqrt(2);
        HSD=sqrt(omsd)*(wgn(1,1,0)+j*wgn(1,1,0))/sqrt(2);
        R_R=X(symb)*sqrt(Ps)*HSR+NR;
        YSD=X(symb)*sqrt(Ps)*HSD+NSD;
        %$$$$$$$$-------------- DF processing ----------------------------------
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        status=zeros(1,L);
        for k=1:L
            d=abs(R_R(k)*HSR(k)'-X);
            [val ind]=min(d);
            if ind == symb
                status(k)=1;
            end
        end
        Powused=Pow.*status;
        Y=X(symb)*sqrt(Powused').*HRD+ND;
        alpha= (sqrt(Powused').*HRD)'*Y+(sqrt(Ps).*HSD)'*YSD;
        d=abs(alpha-X);
        [val ind]=min(d);
        if  symb ~= ind
            counter=counter+1;
        end
    end
        SEPexact(sayan)=S;
        SEPsim(sayan)=counter/iter;
        Pdb(sayan)=Power;
end
semilogy(Pdb,(SEPexact));
hold on
semilogy(Pdb,(SEPsim),'ro');