function Throughput=CSI_Ray_SNR(K,MC_run_1)
%causal side information
tic
B_step=0.01;
loop=0;
itr_num=MC_run_1;
for SNR = 0:5:20 %SNR = 0:20 db
    loop = loop+1;
    SNR_rms = 10^((1/10)*SNR); %convert dB to rms
    %Fading
    r = SNR_rms*ones(1,K);
    %% K=1
    if K==1
        for itr_B1 = 1:itr_num
            flc = rand(1); %flip the coin
            if flc <=1/3
                B1 = 0.5;
            elseif flc<=2/3
                B1 = 1;
            elseif flc<=1
                B1 = 0;
            end
            E1_fun=@(x) exp(-x)./x;
            if B1==0
                J1_test(itr_B1)=0;
            else
                J1_test(itr_B1)=exp(1/(SNR_rms*B1))*integral(E1_fun,1/(SNR_rms*B1),inf)/log(2);
            end
        end
        J1(loop)=mean(J1_test);
    else
        %% K=2 or K=4
        for i = K:-1:1
            % first iteration : compute J_K
            if i == K
                for B_num=1:( floor(i/B_step)+1 )
                    B=(B_num-1)*B_step;
                    E1_fun=@(x) exp(-x)./x;
                    if B==0
                        J(i,B_num)=0;
                    else
                        J(i,B_num)=exp(1/(SNR_rms*B))*integral(E1_fun,1/(SNR_rms*B),inf)/log(2);
                    end
                end
                % iteration : compute J_k, k=K-1, K-2, ..., 2
            elseif i>1 & i<K
                for B_num=1:( floor(i/B_step)+1 )
                    B=(B_num-1)*B_step;
                    % compute average J_k+1
                    J_avg(i+1,B_num)=0;
                    for itr = 1:itr_num
                        flc = rand(1); %flip the coin
                        if flc <=1/3
                            H = 0.5;
                        elseif flc<=2/3
                            H = 1;
                        elseif flc<=1
                            H = 0;
                        end
                        B_temp=floor((B+H)/B_step)+1;
                        J_avg(i+1,B_num)=J_avg(i+1,B_num)+J(i+1,B_temp);
                    end
                    J_avg(i+1,B_num)=J_avg(i+1,B_num)/itr_num;
                    % compute J_k, using bisection search method
                    J_fun=@(x) exp(1/(SNR_rms*x+2e-3))*integral(E1_fun,1/(SNR_rms*x+2e-3),inf)/log(2)+J_avg(i+1,floor((B-x)/B_step)+1);
                    [J(i,B_num),T(i,B_num) ] = bisection(J_fun,0,B);
                end
                % iteration : compute J_1, which is maximum throughput
            elseif i==1
                % compute average J_2
                for B_num=1:(floor(1/B_step)+1)
                    B=(B_num-1)*B_step;
                    J_avg(i+1,B_num)=0;
                    for itr = 1:itr_num
                        flc = rand(1); %flip the coin
                        if flc <=1/3
                            H = 0.5;
                        elseif flc<=2/3
                            H = 1;
                        elseif flc<=1
                            H = 0;
                        end
                        B_temp=floor((B+H)/B_step)+1;
                        J_avg(i+1,B_num)=J_avg(i+1,B_num)+J(i+1,B_temp);
                    end
                    J_avg(i+1,B_num)=J_avg(i+1,B_num)/itr_num;
                end
                % random initial state B1
                for itr_B1 = 1:itr_num
                    flc = rand(1); %flip the coin
                    if flc <=1/3
                        B1 = 0.5;
                    elseif flc<=2/3
                        B1 = 1;
                    elseif flc<=1
                        B1 = 0;
                    end
                    J_fun=@(x) exp(1/(SNR_rms*x+2e-3))*integral(E1_fun,1/(SNR_rms*x+2e-3),inf)/log(2)+J_avg(i+1,floor((B1-x)/B_step)+1);
                    % compute J_1, using bisection search method
                    [J1_test(itr_B1),T1_test(itr_B1)] = bisection(J_fun,0,B1);
                end
                J1(loop)=mean(J1_test);
            end
            
        end
    end
end
%% output maximum throughput
Throughput=J1/K;
toc

