%MHDOPP18.m is the main program for buffer-aided transmission with OR in a
%three node network
%Stage 1. find the trans order. For 3 node network, DN<RN<SN
%Stage 2. find the \bar{eng} with Effb=0;
%Stage 3. find the final result, the changes of Effb is based on systemerror
%this is only for 3 nodes. eng, TS/pac is right.
%Written by Chen Dong --Southampton University
clear
clc
tic
%parameters, maxpower is very large,e.g 1e-02
node(1).x=100;node(1).y=100;
node(2).x=1100;node(2).y=100;
node(3).x=400;node(3).y=300;%nearer to DN
total_node_num=length(node);
B=[16];%only for sim, update later
sim_number=100*B;%100B is OK, only for sim, update later
naka_m=1;
count_output=1;
clear output
output=[];
pla=2;kappa=(3e8/2.4e09/4/pi)^2;Noise=1e-14;th=10^(2/10);
maxpower=[30]*1e-05;%3e-04 
%k_power=3;%%%for thput eng relationship
%resultMHDOPP(1:2,:)= [2.7682e-02   8.0537e-07   1.4354e-05;  4.5726e-02   1.5803e-06   1.6758e-05];%%%for thput eng relationship
%maxpower=maxpower+(resultMHDOPP(k_power-1,3)-resultMHDOPP(k_power-2,3))/(resultMHDOPP(k_power-1,1)-resultMHDOPP(k_power-2,1))*0.02;%%%for thput eng relationship
%while maxpower<1e-02%%%for thput eng relationship
    clear the_out_temp the_out temppac tempTS temp1 temp2 temp systemerrortemp systemerror_pernode systemerror_old P E
for i=1:total_node_num
    node(i).energy=0;
    node(i).energy_avg=0;
    node(i).received=zeros(total_node_num,total_node_num);
    node(i).transmitted=zeros(total_node_num,total_node_num);
    node(i).b=0;
    node(i).B=B;
    node(i).outtime=0;
end
for i=1:total_node_num
    for j=1:total_node_num
        dis(i,j)=((node(i).x-node(j).x)^2+(node(i).y-node(j).y)^2)^0.5;
    end
end
C=th*Noise/kappa.*(dis.^pla);
th_out=C/maxpower;
node_order=zeros(1,total_node_num);
[temp node_order]=sort(dis(:,2));%this is not always right. It should be a interational code. first find the mininum distance, then each round the potential R is considered as SN, find the min-eng with all determined RN, DN, finally, all nodes are determined. All are the similar as OR(zuojing7.m) code.
Icood=[1 2;1 3;3 2];
%2 3 4 1
k=1;
for i=length(node_order):-1:2
    for j=1:i-1
        Icood(k,1)=node_order(i);
        Icood(k,2)=node_order(j);k=k+1;
    end
end
node_orderM=zeros(total_node_num,total_node_num);
for i=1:length(Icood)
    node_orderM(Icood(i,1),Icood(i,2))=1;
end


%%%%%%%%%theory, based on a roughly random node(i).energy_avg, everything
%%%%%%%%%is OK
P=zeros(total_node_num,total_node_num);P_old=P;
for ii=1:30
%%%%all P eng_avg written in integral
%%division points
for i=1:length(Icood)
    E_avg(i)=(node(1).energy_avg-node(Icood(i,1)).energy_avg)+(node(Icood(i,2)).energy_avg-node(2).energy_avg);
end


P=zeros(total_node_num,total_node_num);
E=zeros(total_node_num,total_node_num);
 
divp(1)=th_out(1,2);
for i=2:length(Icood)
if E_avg(i)>0
   divp(i)=C(1,2)/E_avg(i);%+-?upper
else
   divp(i)=10;
end
end
divp=min(divp,10);
if sum(divp<0)>0%only for debug
    disp('error divp, stage1, SD part')
    break
end%end of only for debug
[Y,I]=sort(divp);
Iflag=ones(1,length(Icood)+1);
I(length(Icood)+1)=length(Icood)+1;Y(length(Icood)+1)=10;
for i=1:length(Icood)
    Iflag(I(i))=0;
    P(Icood(I(1),1),Icood(I(1),2))=P(Icood(I(1),1),Icood(I(1),2))+quadgk(@(x)exp(-x).*...
        (1-exp(-C(Icood(I(3),1),Icood(I(3),2))./(C(1,2)./x-E_avg(I(3))))).^Iflag(I(3)).*...
        (1-exp(-C(Icood(I(2),1),Icood(I(2),2))./(C(1,2)./x-E_avg(I(2))))).^Iflag(I(2)),Y(i),Y(i+1));
    E(Icood(I(1),1),Icood(I(1),2))=E(Icood(I(1),1),Icood(I(1),2))+quadgk(@(x)C(1,2)./x.*exp(-x).*...
        (1-exp(-C(Icood(I(3),1),Icood(I(3),2))./(C(1,2)./x-E_avg(I(3))))).^Iflag(I(3)).*...
        (1-exp(-C(Icood(I(2),1),Icood(I(2),2))./(C(1,2)./x-E_avg(I(2))))).^Iflag(I(2)),Y(i),Y(i+1));
end


for Icood_count=2:length(Icood)
    if E_avg(Icood_count)<0
       divp(1)=-C(Icood(Icood_count,1),Icood(Icood_count,2))/E_avg(Icood_count);%upper
    else
       divp(1)=10;
    end
divp(Icood_count)=th_out(Icood(Icood_count,1),Icood(Icood_count,2));%
for i=2:length(Icood)
    if i~=Icood_count
       if E_avg(i)-E_avg(Icood_count)>0
          divp(i)=C(Icood(Icood_count,1),Icood(Icood_count,2))/(E_avg(i)-E_avg(Icood_count));%upper
       else
          divp(i)=10;
       end
    end
end
if sum(divp<0)>0%only for debug
    disp('error divp, stage1, non-SD part')
    break
end%end of only for debug
divp=min(divp,10);
[Y,I]=sort(divp);
Iflag=ones(1,length(Icood)+1);
I(length(Icood)+1)=length(Icood)+1;Y(length(Icood)+1)=10;
for i=1:length(Icood)
Iflag(I(i))=0;
P(Icood(Icood_count,1),Icood(Icood_count,2))=P(Icood(Icood_count,1),Icood(Icood_count,2))+quadgk(@(x)exp(-x).*...
    (1-exp(-C(Icood(1,1),Icood(1,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(1)))).^Iflag(1).*...
    (1-exp(-C(Icood(2,1),Icood(2,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(2)))).^Iflag(2).*...
    (1-exp(-C(Icood(3,1),Icood(3,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(3)))).^Iflag(3),Y(i),Y(i+1));
E(Icood(Icood_count,1),Icood(Icood_count,2))=E(Icood(Icood_count,1),Icood(Icood_count,2))+quadgk(@(x)C(Icood(Icood_count,1),Icood(Icood_count,2))./x.*exp(-x).*...
    (1-exp(-C(Icood(1,1),Icood(1,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(1)))).^Iflag(1).*...
    (1-exp(-C(Icood(2,1),Icood(2,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(2)))).^Iflag(2).*...
    (1-exp(-C(Icood(3,1),Icood(3,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)-E_avg(3)))).^Iflag(3),Y(i),Y(i+1));
end
end

%P
systemerror=0;
for i=2:length(node_order)-1
    temp1=0;temp2=0;
    for j=length(node_order):-1:i+1
        temp1=temp1+P(node_order(j),node_order(i));
    end
    for j=1:i-1
        temp2=temp2+P(node_order(i),node_order(j));
    end
    systemerror=systemerror+abs(temp1-temp2);
end
for i=2:length(node_order)%2 3 4 1
    if i==2
       node(node_order(i)).energy_avg=E(node_order(i),2)/P(node_order(i),2);
    else
       E_temp=E(node_order(i),2);  
       P_temp=P(node_order(i),2);
       for j=2:i-1
           E_temp=E_temp+E(node_order(i),node_order(j))+P(node_order(i),node_order(j))*node(node_order(j)).energy_avg;
           P_temp=P_temp+P(node_order(i),node_order(j));
       end
       node(node_order(i)).energy_avg=E_temp/P_temp;
    end
end
%systemerror
%node(1).energy_avg
if sum(sum(abs(P-P_old)))<1e-06
   break
end
P_old=P;
end

Effb=node(1).energy_avg/10;%base step of Effb
aim_err=0.001;pb_counter=1;
signs=zeros(1,length(Icood));  
tiny=1e-08;
while systemerror>aim_err*3
for i=1:length(Icood)
    E_avg(i)=(node(1).energy_avg-node(Icood(i,1)).energy_avg)+(node(Icood(i,2)).energy_avg-node(2).energy_avg);
end
 
P=zeros(total_node_num,total_node_num);
E=zeros(total_node_num,total_node_num);
   %%list all Effbchange_temp
   systemerror=0;
   for i=2:length(node_order)-1
       temp1=0;temp2=0;
       for j=length(node_order):-1:i+1
           temp1=temp1+P_old(node_order(j),node_order(i));
       end
       for j=1:i-1
           temp2=temp2+P_old(node_order(i),node_order(j));
       end
       systemerror=systemerror+abs(temp1-temp2);
   end
   Effbchange_temp=zeros(length(Icood),length(Icood));
   Effbchange_temp_factor=zeros(length(Icood),length(Icood));
   for i=1:length(Icood)
   for j=1:length(Icood)
       P_old_temp=P_old;
       P_old_temp(Icood(i,1),Icood(i,2))=P_old_temp(Icood(i,1),Icood(i,2))-tiny;
       P_old_temp(Icood(j,1),Icood(j,2))=P_old_temp(Icood(j,1),Icood(j,2))+tiny;
       systemerrortemp=0;
       for i1=2:length(node_order)-1
           temp1=0;temp2=0;
           for j1=length(node_order):-1:i1+1
               temp1=temp1+P_old_temp(node_order(j1),node_order(i1));
           end
           for j1=1:i1-1
               temp2=temp2+P_old_temp(node_order(i1),node_order(j1));
           end
           systemerrortemp=systemerrortemp+abs(temp1-temp2);
       end
       Effbchange_temp_factor(i,j)=(systemerrortemp-systemerror)/tiny;
   end
   end
   Effbchange_temp_factor=round(Effbchange_temp_factor);
   Effbchange=zeros(1,length(Icood));
   for i=2:length(Icood)
       if sum(Effbchange_temp_factor(i,:)<0)==0
          Effbchange(i)=systemerror;
       end
   end
   for i=2:length(Icood)
       if sum(Effbchange_temp_factor(i,:)>0)==0
          Effbchange(i)=-systemerror;
       end
   end
     
   signs=signs-Effbchange;

   for availablechannels_num=1:2^length(Icood)-1
       availablechannels_num_temp=dec2bin(availablechannels_num,length(Icood));
       for j=1:length(Icood)
           availablechannels(j)=bin2dec(availablechannels_num_temp(j));%bound all channels available
       end
   factor_temp=1;
   for i=1:length(Icood)
       factor_temp=factor_temp*exp(-th_out(Icood(i,1),Icood(i,2)))^0*(1-exp(-th_out(Icood(i,1),Icood(i,2))))^(1-availablechannels(i));
   end
   [P_part,E_part]=MHDOPP18PE(total_node_num,C,Icood,signs,Effb,th_out,E_avg,availablechannels);P=P+factor_temp*P_part;E=E+factor_temp*E_part;
   end
   
   
% divp(1)=th_out(1,2);
% for i=2:length(Icood)
% if E_avg(i)+signs(i)*Effb>0
%    divp(i)=C(1,2)/(E_avg(i)+signs(i)*Effb);%+-?upper
% else
%    divp(i)=10;
% end
% end
% divp=min(divp,10);
% if sum(divp<0)>0%only for debug
%     disp('error divp, stage2, SD part')
%     break
% end%end of only for debug
% [Y,I]=sort(divp);
% Iflag=ones(1,length(Icood)+1);
% I(length(Icood)+1)=length(Icood)+1;Y(length(Icood)+1)=10;
% for i=1:length(Icood)
% Iflag(I(i))=0;
% P(Icood(I(1),1),Icood(I(1),2))=P(Icood(I(1),1),Icood(I(1),2))+quadgk(@(x)exp(-x).*...
%     (1-exp(-C(Icood(I(3),1),Icood(I(3),2))./(C(1,2)./x-E_avg(I(3))-signs(I(3))*Effb))).^Iflag(I(3)).*...
%     (1-exp(-C(Icood(I(2),1),Icood(I(2),2))./(C(1,2)./x-E_avg(I(2))-signs(I(2))*Effb))).^Iflag(I(2)),Y(i),Y(i+1),'AbsTol',1e-06,'RelTol',0);
% E(Icood(I(1),1),Icood(I(1),2))=E(Icood(I(1),1),Icood(I(1),2))+quadgk(@(x)C(1,2)./x.*exp(-x).*...
%     (1-exp(-C(Icood(I(3),1),Icood(I(3),2))./(C(1,2)./x-E_avg(I(3))-signs(I(3))*Effb))).^Iflag(I(3)).*...
%     (1-exp(-C(Icood(I(2),1),Icood(I(2),2))./(C(1,2)./x-E_avg(I(2))-signs(I(2))*Effb))).^Iflag(I(2)),Y(i),Y(i+1),'AbsTol',1e-10,'RelTol',0);
% end
% 
% for Icood_count=2:length(Icood)
%     if (E_avg(Icood_count)+signs(Icood_count)*Effb)<0
%        divp(1)=-C(Icood(Icood_count,1),Icood(Icood_count,2))/(E_avg(Icood_count)+signs(Icood_count)*Effb);%upper
%     else
%        divp(1)=10;
%     end
% divp(Icood_count)=th_out(Icood(Icood_count,1),Icood(Icood_count,2));%+-?
% for i=2:length(Icood)
%     if i~=Icood_count
%        if E_avg(i)+signs(i)*Effb-E_avg(Icood_count)-signs(Icood_count)*Effb>0
%           divp(i)=C(Icood(Icood_count,1),Icood(Icood_count,2))/(E_avg(i)+signs(i)*Effb-E_avg(Icood_count)-signs(Icood_count)*Effb);%upper
%        else
%           divp(i)=10;
%        end
%     end
% end
% divp=min(divp,10);
% if sum(divp<0)>0%only for debug
%     disp('error divp, stage2, other part')
%     break
% end%end of only for debug
% [Y,I]=sort(divp);
% Iflag=ones(1,length(Icood)+1);
% I(length(Icood)+1)=length(Icood)+1;Y(length(Icood)+1)=10;
% for i=1:length(Icood)
% Iflag(I(i))=0;
% P(Icood(Icood_count,1),Icood(Icood_count,2))=P(Icood(Icood_count,1),Icood(Icood_count,2))+quadgk(@(x)exp(-x).*...
%     (1-exp(-C(Icood(1,1),Icood(1,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(1)-signs(1)*Effb))).^Iflag(1).*...
%     (1-exp(-C(Icood(2,1),Icood(2,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(2)-signs(2)*Effb))).^Iflag(2).*...
%     (1-exp(-C(Icood(3,1),Icood(3,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(3)-signs(3)*Effb))).^Iflag(3),Y(i),Y(i+1),'AbsTol',1e-06,'RelTol',0);
% E(Icood(Icood_count,1),Icood(Icood_count,2))=E(Icood(Icood_count,1),Icood(Icood_count,2))+quadgk(@(x)C(Icood(Icood_count,1),Icood(Icood_count,2))./x.*exp(-x).*...
%     (1-exp(-C(Icood(1,1),Icood(1,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(1)-signs(1)*Effb))).^Iflag(1).*...
%     (1-exp(-C(Icood(2,1),Icood(2,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(2)-signs(2)*Effb))).^Iflag(2).*...
%     (1-exp(-C(Icood(3,1),Icood(3,2))./(C(Icood(Icood_count,1),Icood(Icood_count,2))./x+E_avg(Icood_count)+signs(Icood_count)*Effb-E_avg(3)-signs(3)*Effb))).^Iflag(3),Y(i),Y(i+1),'AbsTol',1e-06,'RelTol',0);
% end
% end

P
systemerror_old=systemerror;
systemerror=0;
for i=2:length(node_order)-1
    temp1=0;temp2=0;
    for j=length(node_order):-1:i+1
        temp1=temp1+P(node_order(j),node_order(i));
    end
    for j=1:i-1
        temp2=temp2+P(node_order(i),node_order(j));
    end
    systemerror_pernode(node_order(i))=abs(temp1-temp2);
    systemerror=systemerror+abs(temp1-temp2);
end
for i=2:length(node_order)%2 3 4 1
    if i==2
       node(node_order(i)).energy_avg=E(node_order(i),2)/P(node_order(i),2);
    else
       E_temp=E(node_order(i),2);  
       P_temp=P(node_order(i),2);
       for j=2:i-1
           E_temp=E_temp+E(node_order(i),node_order(j))+P(node_order(i),node_order(j))*node(node_order(j)).energy_avg;
           P_temp=P_temp+P(node_order(i),node_order(j));
       end
       node(node_order(i)).energy_avg=E_temp/P_temp;
    end
end
systemerror
node(1).energy_avg
if sum(sum(abs(P-P_old)))<1e-06
   break
end
P_old=P;
if systemerror<0.05%stat: pridict_bound
   for i=length(node_order):-1:2
       pridict(pb_counter,(length(node_order)-i)*2+2)=node(node_order(i)).energy_avg;
       pridict(pb_counter,(length(node_order)-i)*2+1)=systemerror_pernode(node_order(i));
   end
   pridict(pb_counter,1)=systemerror;
   pb_counter=pb_counter+1;
end
end

%%%pridict bound
if exist('pridict')~=0%for large probability of outage, only for 3nodes
for i=length(node_order):-1:2
temp=polyfit(pridict(:,2*(length(node_order)-i)+1),pridict(:,2*(length(node_order)-i)+2),1);%linear
node(node_order(i)).energy_avg=temp(2);
end
end
%%%end of pridict
%P
[node(1).energy_avg node(3).energy_avg]
%systemerror
the_out_temp=1-exp(-th_out);the_out=1;
for i=1:length(Icood)
    the_out=the_out*the_out_temp(Icood(i,1),Icood(i,2));
end
disp('out only based on channel')
sum(the_out)

for i=2:length(node_order)%2 3 4 1
    if i==2
       node(node_order(i)).R2Dthput=1;%1TS/1pac
    else 
       clear tempTS temppac
       tempTS=P(node_order(i),2);
       temppac=P(node_order(i),2);
       for j=2:i-1
           tempTS=tempTS+P(node_order(i),node_order(j))*(1+node(node_order(j)).R2Dthput);
           temppac=temppac+P(node_order(i),node_order(j));
       end
       node(node_order(i)).R2Dthput=tempTS/temppac;
    end
end
disp('the thput')
node(1).R2Dthput=node(1).R2Dthput/(1-sum(the_out));%thput when trans/out
node(1).R2Dthput
disp('time for theory')
toc


%resultMHDOPP(k_power,1)=P(1,2)+P(1,3);%%%for thput eng relationship
%resultMHDOPP(k_power,2)=sum(sum(E));%%%for thput eng relationship
%resultMHDOPP(k_power,3)=maxpower;%%%for thput eng relationship
%k_power=k_power+1;%%%for thput eng relationship
%resultMHDOPP
%maxpower=maxpower+(resultMHDOPP(k_power-1,3)-resultMHDOPP(k_power-2,3))/(resultMHDOPP(k_power-1,1)-resultMHDOPP(k_power-2,1))*0.02;%%%for thput eng relationship
%end%enf of while %%%for thput eng relationship
%resultMHDOPP



%%%%%%%%%%%%%simulation
output
for B=[8];%only for sim
     if B>64
     equal=1;%find the exact solution when B<=64
     elseif B<=64
     equal=0;
     end
     [the_out_exact,P_exact,E_exact,tao_the]=MHDOPP12outage(total_node_num,node,B,C,Icood,signs,Effb,th_out,E_avg,equal);%theory
     node(1).energy_avg_exact=sum(sum(E_exact))/(P_exact(1,2)+(sum(P_exact(1,3:end))+sum(P_exact(3:end,2)))/2);
    for sim_round=1:1
    for i=1:total_node_num
       node(i).B=B;
    end  
    sim_number=max(400*B,32*500);%100B is OK, only for sim%minimum 300*
    TS=1;
    [node,TS]=MHDOPP15(node,C,maxpower,total_node_num,naka_m,sim_number,Icood,signs,Effb,node_orderM,th_out);
    P_sim=zeros(total_node_num,total_node_num);
    node(1).energy_avg_sim=0;
    for i=1:total_node_num
       P_sim=P_sim+node(i).transmitted/TS;
       node(1).energy_avg_sim=node(1).energy_avg_sim+node(i).energy_real;
    end
    P_sim
    disp('simulation:sim/the')%the principle is that the PMF is even at all RNs. only consider 1 channel is unavailable at one time. the energy consumption rises to node(1).energy_avg*all_channel_num/(all_channel_num-1).
    output(count_output,:)=[maxpower B node(1).energy_avg_sim/sim_number node(1).energy_avg_exact node(1).energy_avg  sum(node(1).out_history(50*B+1:end-50*B)==1)/(TS-100*B) the_out_exact sum(the_out) TS/sim_number 1/(P_exact(1,2)+(sum(P_exact(1,3:end))+sum(P_exact(3:end,2)))/2) node(1).R2Dthput 1/(exp(-th_out(1,2))/2+1/2) signs*Effb]%B, sim_eng the_eng_exact the_eng_bound sim_out the_out_exact the_out_bound sim_thput the_thput_exact the_thput_bound upper_thput_bound effb*signs
    count_output=count_output+1;
    disp('time for theory+sim')
    toc   
    end%sim_round
end
tao=0;
for i=1:length(node(1).delayPMF)
tao=tao+i*node(1).delayPMF(i);
end
tao=tao/sim_number;
%[tao tao_the]
disp('maxpower Buffer_size, sim_PED theoretical_exact_PED theoretical_PED_bound sim_outage theoretical_exact_out theoretical_out_bound sim_throughput theoretical_exact_throughput theoretical_throughput_bound')
output(:,1:11)

%B
% 
% %%%%only for thput-eng figure
% %%%direct transmission
% k=1;
% for thput=0.02:0.02:1
%     temp=reallog(1/thput);
%     resultdir(k,1)=thput;
%     resultdir(k,2)=C(1,2)*expint(temp);
%     k=k+1;
% end
% %resultdir
% 
%%%%%conventional 2hop
% k=1;
% %power allocation
% PAfactor=dis(1,2)^2/(dis(1,3)^2+dis(3,2)^2);
% for thput=0.02:0.02:0.48
%     temp=reallog(1/(2*thput));
%     result2hop(k,1)=thput;
%     result2hop(k,2)=(C(1,2)+C(3,2))*expint(temp)/PAfactor;
%     k=k+1;
% end
% result2hop
% 
% %%%%%MHD, only add a diversity
% k=1;
% %power allocation
% PAfactor=dis(1,2)^2/(dis(1,3)^2+dis(3,2)^2);
% for thput=0.02:0.02:0.48
%     temp=reallog(1/(1-(1-2*thput)^(1/2)));
%     result2hopMHD(k,1)=thput;
%     result2hopMHD(k,2)=(C(1,2)+C(3,2))*2*(expint(temp)-expint(temp*2))/PAfactor;
%     k=k+1;
% end
% result2hopMHD


%%%%%conventional OR

% k=1;
% maxpower=2e-05;
% while maxpower<1e-02%%%for thput eng relationship
%     maxpower=maxpower*1.15;
% th_out=C/maxpower;
% %%energy distance of R
% node(3).energy_avg=expint(th_out(3,2))*C(3,2);
% %%probability of SR when S select
% tempE(1,3)=quadgk(@(r)exp(-r).*(1-exp(-C(1,2)./(C(1,3)./r+node(3).energy_avg)))./r,max(th_out(1,3),C(1,3)/(C(1,2)/th_out(1,2)-node(3).energy_avg)),10)*C(1,3);
% tempP(1,3)=quadgk(@(r)exp(-r).*(1-exp(-C(1,2)./(C(1,3)./r+node(3).energy_avg))),max(th_out(1,3),C(1,3)/(C(1,2)/th_out(1,2)-node(3).energy_avg)),10);
% 
% tempE(1,2)=quadgk(@(r)exp(-r).*(1-exp(-C(1,3)./(C(1,2)./r-node(3).energy_avg)))./r,max(th_out(1,2),C(1,2)/(C(1,3)/th_out(1,3)+node(3).energy_avg)),C(1,2)/node(3).energy_avg)*C(1,2)+...
%            expint(C(1,2)/node(3).energy_avg)*C(1,2);
% tempP(1,2)=quadgk(@(r)exp(-r).*(1-exp(-C(1,3)./(C(1,2)./r-node(3).energy_avg))),max(th_out(1,2),C(1,2)/(C(1,3)/th_out(1,3)+node(3).energy_avg)),C(1,2)/node(3).energy_avg)+...
%            exp(-C(1,2)/node(3).energy_avg);
% 
% node(1).energy_avg=(tempE(1,2)+tempE(1,3)+tempP(1,3)*node(3).energy_avg)/(tempP(1,2)+tempP(1,3));%energy per packet
% thput=(tempP(2)+tempP(3))/(1+tempP(3)/exp(-th_out(3,2)));
% resultOR(k,1)=thput;
% resultOR(k,2)=(tempE(1,2)+tempE(1,3)+tempP(1,3)*node(3).energy_avg)/(1+tempP(3)/exp(-th_out(3,2)));
% k=k+1;
% end
% resultOR