function I_wf_AWGN = FSI_TS(MC_run,AWGN_Fading_indicator)
tic
K = 8;
%full CSI
%H(1) = B_0, H(2) = H_1, H(3) = H_2, and etc.
%H_K is not used (casuality of harvasting energy.)
I_avg = zeros(MC_run,8);
SNR = 20;
for itr = 1:MC_run
    %AWGN
    tol = 10^(-4);
    S = [];
    loop = 0;
    for K_ind = 1:K %SNR = 0:20 db
        loop = loop+1;
        SNR_rms = 10^((1/10)*SNR); %convert dB to rms 
        K_num = 2^(K_ind-1);
        if AWGN_Fading_indicator == 1 %%AWGN
            r = SNR_rms*ones(1,K_num);
        elseif AWGN_Fading_indicator == 2 %%Fading
            r = exprnd(SNR_rms,[1,K_num]);
        end        
        H = zeros(1,K_num);
        for k_i = 1:K_num
            flc = rand(1); %flip the coin
            if flc <=1/3
                H(1,k_i) = 0.5;
            elseif flc<=2/3
                H(1,k_i) = 1;
            elseif flc<=1
                H(1,k_i) = 0;
            end
        end
        T_opt = zeros(1,K_num);
        for i = 1:K_num
            % Outer iteration: find t1*, t2* (TS)
            if i == 1
                t_in1 = 0; %First entering only
            end

            for k = K_num: -1: t_in1+1
                T = zeros(1,K_num);
                T(1:t_in1) = T_opt(1:t_in1);

                P_max = sum(H(t_in1+1:k)); %note that H_n = H(n+1), summing H from t_i-1 to k-1
                [T_out,lambda] = con_water_filling(r(t_in1+1:k), P_max,tol);
                T(t_in1+1:k) = T_out;
                [ind_constraint] = constraint_check(T(1:k),H(1:k));
                if ind_constraint == 1
                    S = [S,k];
                    t_in1 = k;
                    T_opt(1:t_in1) = T(1:t_in1);
                    break
                end
            end
            if t_in1 == K_num
                break
            end
        end
        I  = log2(1+T.*r);
        I_avg(itr,loop)= sum(I)/K_num;

    end
end
I_wf_AWGN = mean(I_avg,1);
toc


