%change SNR_hbh and press F5


%%originally, this is for response for TWC paper. 
%Now, it determines the thresholds of MHD-AMC without power contral.
%This is for a L-hop link.
%Author: Chen Dong
clear
clc
tic
%%%%%%%%Input SNR of each hop in EbN0
SNR_hbh=[20 10];%EbN0, set average SNR, distance is useless.
hop_number=length(SNR_hbh);
SNR_hbh_dB=10*log10(SNR_hbh);%in dB

%%%%%%Input end-to-end targe BER
BER_target_e2e=0.03;
noise=1e-14;
pla=2;
kappa=((3e8/2.4e9)/4/pi)^pla;%all parameters are from Jing's paper.

%%%for compare with power control, now useless
%distance=[1 2];%also set the hop number
%hop_number=length(distance);
%distance=distance.*1000/sum(distance);%total distance is 1000m; useless in this code
%total_power=3e-04;
%power=total_power*distance.^pla/sum(distance.^pla);
%SNR_hbh=power./distance.^pla/noise*kappa;


%%%find the thresholds%each hop one TS, totally 2TS
BER_tar=ones(1,hop_number)*BER_target_e2e/hop_number;%original set
seg_basic=0.5;%channel SNR step length
M=[2 4 16 64];%modulation schemes
th=zeros(hop_number,4);%thresholds
restemp_count=1;
while 1
for hop_num_count=1:hop_number%find thresholds with a given tar_BER vector
for M_count=1:length(M)
    seg=seg_basic;
    i=0;
    while seg>0.0001
    while berfadingamc(SNR_hbh_dB(hop_num_count)-10*log10(log2(M(M_count))),M(M_count),1,i+seg,i+seg,1)>BER_tar(hop_num_count)&&i<10
       i=i+seg;
    end
    seg=seg/2;
    end
    th(hop_num_count,M_count)=i;
end
end
th%show current thresholds
thput=zeros(1,hop_number);
for i=1:hop_number
    thput(i)=exp(-th(i,1))+exp(-th(i,2))+2*exp(-th(i,3))+2*exp(-th(i,4));
end
thput%show current throughput

%%%%%%adjust target BER, 
%Increase the target BER of the hop with minimum throughput
%Decrease the target BER of other hops
min_thput_hop=find(thput==min(thput),1);
seg_BER=(max(thput)-min(thput));%BER adjustment step, decrease with running the code
seg_BER=min(seg_BER,0.2*(sum(BER_tar)-BER_tar(min_thput_hop))/BER_tar(min_thput_hop));%the maximum BER adjustment is 20% of each BER target value
BER_tar_new=zeros(1,hop_number);
BER_tar_new=BER_tar-BER_tar(min_thput_hop)*seg_BER*BER_tar/(sum(BER_tar)-BER_tar(min_thput_hop));
BER_tar_new(min_thput_hop)=BER_tar(min_thput_hop)*(1+seg_BER);

%%% check point when coding, now useless
%while sum(BER_tar_new<0)>0
%    seg_BER=seg_BER/2;
%    BER_tar_new=zeros(1,hop_number);
%BER_tar_new=BER_tar-BER_tar(min_thput_hop)*seg_BER*BER_tar/(sum(BER_tar)-BER_tar(min_thput_hop));
%BER_tar_new(min_thput_hop)=BER_tar(min_thput_hop)*(1+seg_BER);
%end
%if sum(BER_tar_new)>1.0001*BER_target_e2e
%    [BER_tar_new]
%    break
%end
BER_tar=BER_tar_new;
th_final=max(th);
result_temp(restemp_count)=exp(-th_final(1))+exp(-th_final(2))+2*exp(-th_final(3))+2*exp(-th_final(4));restemp_count=restemp_count+1;
if sum(result_temp>sum(max(result_temp)*0.995))>=5&&(max(thput)-min(thput))/min(thput)<0.005%stop condition:1) each hop has similar throughput 2)the throughput nearly constant
    break
end
%[result_temp(end) BER_tar seg_BER]
end
disp('The hop throughput are')
result_temp(end)%achievable throoughput
disp('The targe BER of each hop are')
BER_tar
disp('The thresholds are')
max(th)
toc