function [P,E]=MHDOPP18PE(total_node_num,C,Icood,signs,Effb,th_out,E_avg,availablechannels)
%part of MHDOPP12.m
%calculate P and E based on existing Effb,E_avg,availablechannels
P=zeros(total_node_num,total_node_num);
E=zeros(total_node_num,total_node_num);

for Icood_count=1:length(Icood)
    if availablechannels(Icood_count)==1%if this channel is available
divp(Icood_count)=th_out(Icood(Icood_count,1),Icood(Icood_count,2));%+-?
for i=1: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
Iflag=availablechannels;%put available hop here
divp=divp+(Iflag==0)*20;
[Y,I]=sort(divp);
Y=min(10,Y);
I(length(Icood)+1)=length(Icood)+1;Y(length(Icood)+1)=10;
Iflag(I(1:find(I==Icood_count)-1))=zeros;%set all divp<th_out(1,2) as 0. which means regardless the channel quality of that channel
for i=find(I==Icood_count):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).*...
    (exp(-th_out(Icood(1,1),Icood(1,2)))-exp(-max(th_out(Icood(1,1),Icood(1,2)),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).*...
    (exp(-th_out(Icood(2,1),Icood(2,2)))-exp(-max(th_out(Icood(2,1),Icood(2,2)),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).*...
    (exp(-th_out(Icood(3,1),Icood(3,2)))-exp(-max(th_out(Icood(3,1),Icood(3,2)),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).*...
    (exp(-th_out(Icood(1,1),Icood(1,2)))-exp(-max(th_out(Icood(1,1),Icood(1,2)),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).*...
    (exp(-th_out(Icood(2,1),Icood(2,2)))-exp(-max(th_out(Icood(2,1),Icood(2,2)),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).*...
    (exp(-th_out(Icood(3,1),Icood(3,2)))-exp(-max(th_out(Icood(3,1),Icood(3,2)),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
end

end%function