function [L]=BCJR(r,v,PS,NS,OP,SIG0,SIGR,Lc)

N=length(PS);
AL(1,1)=0;
AL(2:N,1)=-Inf;

G0(1)=Gamma(0,SIG0(1,1),SIG0(1,2),PS,NS,OP,Lc,r);%*0.5*sum(r(1,:).*Codeword(PS,NS,OP,SIG0(1,1),SIG0(1,2)));%(+ ul*0.5*La(ul))
G0(2)=Gamma(0,SIG0(2,1),SIG0(2,2),PS,NS,OP,Lc,r);%Lc*0.5*sum(r(1,:).*Codeword(PS,NS,OP,SIG0(2,1),SIG0(2,2)));
s=size(r);
h=s(1)-v;

for i=1:h-1  %G(1:K-1,2*2^v) :- Rows 1:K-1=Time ;First 2^v columns=+1 input state pairs;Second 2^v columns=-1 input state pairs
    for j=1:N
    G(i,j)=Gamma(i,SIGR(i,j,1),SIGR(i,j,2),PS,NS,OP,Lc,r);%Lc*0.5*sum(r(i+1,:).*Codeword(PS,NS,OP,SIGR(i,j,1),SIGR(i,j,2)));
    G(i,N+j)=Gamma(i,SIGR(i,N+j,1),SIGR(i,N+j,2),PS,NS,OP,Lc,r);%Lc*0.5*sum(r(i+1,:).*Codeword(PS,NS,OP,SIGR(i,N+j,1),SIGR(i,N+j,2)));
    end
end

for i=h:h+v-1               %G(1:K-1,2*2^v) :- 1:K-1=Time ;First 2^v columns=+1 input ;Second 2^v columns=-1 input
    for j=1:N
    G(i,j)=Gamma(i,SIGR(i,j,1),SIGR(i,j,2),PS,NS,OP,Lc,r);%Lc*0.5*sum(r(i+1,:).*Codeword(PS,NS,OP,SIGR(i,j,1),SIGR(i,j,2)));
    G(i,N+j)=Gamma(i,SIGR(i,N+j,1),SIGR(i,N+j,2),PS,NS,OP,Lc,r);%Lc*0.5*sum(r(i+1,:).*Codeword(PS,NS,OP,SIGR(i,N+j,1),SIGR(i,N+j,2)));
    end
end

for i=1:N
    AL(i,2)=G0(N-i+1)+AL(1,1);
end

for i=3:h+v+1
    for j=1:N
        L=maxstar(-Inf,Gamma(i-2,1,j,PS,NS,OP,Lc,r)+AL(1,i-1));
        for k=2:N
            L=maxstar(L,Gamma(i-2,k,j,PS,NS,OP,Lc,r)+AL(k,i-1));
        end
        AL(j,i)=L;
    end
end

BE(:,h+v+1)=AL(:,h+v+1);

for i=h+v:-1:1
    for j=1:N
        L=maxstar(-Inf,Gamma(i-1,j,1,PS,NS,OP,Lc,r)+BE(1,i+1));
        for k=2:N
            L=maxstar(L,Gamma(i-1,j,k,PS,NS,OP,Lc,r)+BE(k,i+1));
        end
        BE(j,i)=L;
    end
end

for i=1:h+v-1
    L1=maxstar(-Inf,BE(SIGR(i,1,2),i+2)+G(i,1)+AL(SIGR(i,1,1),i+1));
    L2=maxstar(-Inf,BE(SIGR(i,N+1,2),i+2)+G(i,N+1)+AL(SIGR(i,N+1,1),i+1));
    for j=2:N
        L1=maxstar(L1,BE(SIGR(i,j,2),i+2)+G(i,j)+AL(SIGR(i,j,1),i+1));
        L2=maxstar(L2,BE(SIGR(i,N+j,2),i+2)+G(i,N+j)+AL(SIGR(i,N+j,1),i+1));
    end
    L(i+1)=L1-L2;
end

L(1)=(BE(SIG0(1,2),2)+G0(1)+AL(SIG0(1,1),1))-(BE(SIG0(2,2),2)+G0(2)+AL(SIG0(2,1),1));