function mainExit(snrCh,SNR,numOfMemElem1,g1_1,g2_1,numOfMemElem2,g1_2,g2_2,stim_em,Improved_Results,cont,code1,code2,filename)
fh = figure(1);
snrshift = 0;
%Improved_Results = imp                   ;   %I/P
%if Improved_Results
 %   Numbits   = 20000                    ; 
%else
    Numbits   = 20000                     ; 
%end
if code1 == 0
numOfMemElem1 = 1;
g1_1 = 3;
g2_1 = 2;

elseif code1 == 1
numOfMemElem1 = 2;

g1_1 = 7;
g2_1 = 5;

elseif code1 == 2
numOfMemElem1 = 3;

g1_1 = 13;
g2_1 = 15;

elseif code1 == 3
numOfMemElem1 = 4;

g1_1 = 23;
g2_1 = 37;

elseif code1 == 4
numOfMemElem1 = 5;

g1_1 = 67;
g2_1 = 45;

elseif code1 == 5
numOfMemElem1 = 6;

g1_1 = 147;
g2_1 = 117;

elseif code1 == 6
numOfMemElem1 = 8;
g1_1 = 435;
g2_1 = 777;
snrshift = 1;

elseif code1 == 7
numOfMemElem1 = 4;
g1_1 = 37;
g2_1 = 21;

end

if code2 == 0
numOfMemElem2 = 1;

g1_2 = 3;
g2_2 = 2;

elseif code2 == 1
numOfMemElem2 = 2;

g1_2 = 7;
g2_2 = 5;

elseif code2 == 2
numOfMemElem2 = 3;

g1_2 = 13;
g2_2 = 15;

elseif code2 == 3
numOfMemElem2 = 4;

g1_2 = 23;
g2_2 = 37;

elseif code2 == 4
numOfMemElem2 = 5;

g1_2 = 67;
g2_2 = 45;

elseif code2 == 5
numOfMemElem2 = 6;

g1_2 = 147;
g2_2 = 117;

elseif code2 == 6
numOfMemElem2 = 8;
g1_2 = 435;
g2_2 = 777;
snrshift = 1;

elseif code2 == 7
numOfMemElem2 = 4;
g1_2 = 37;
g2_2 = 21;
end
if snrCh && snrshift
SNR = [1.1 1.5 2 2.5];
elseif snrCh
SNR = 1:0.5:2.5;
end

if (numOfMemElem1 == 1 || numOfMemElem2 == 1) && snrCh
SNR = [2.6 3 3.5 4];
end
 
global NumberOfbits 
global Ia           
global BitsUncoded  
global histLimit    
if ( (length(SNR) == 1) && (SNR<1 && SNR>-1)) || ( (length(SNR)>1) && min(SNR)<1 && any(ceil(SNR(SNR<1))~=SNR(SNR<1)))
Ia      =  0:0.05:1                     ;
stim_em =  1                            ;

elseif Improved_Results
	 Ia = 0:0.05:1								 ;
	 
	else
    Ia  = 0:0.1:1                       ; 
end
NumberOfbits = Numbits                  ;
BitsUncoded = randi(2,1,NumberOfbits)-1 ;                 
histLimit   = 32                        ;                 


% Generate Trellis
if numOfMemElem1 == numOfMemElem2 && g1_1 == g1_2 && g2_1 == g2_2
    Identical_decoders = 1;    
else
    Identical_decoders = 0;
end
if Identical_decoders == 1
    [trellis, template, ~, notPuncInd, ~, BitsBPSK]     = Construct(numOfMemElem1, g1_1, g2_1, 1) ;
else
    [trellis1, template1, ~, notPuncInd1, ~, BitsBPSK1] = Construct(numOfMemElem1, g1_1, g2_1, 1) ;
    [trellis2, template2, ~, notPuncInd2, ~, BitsBPSK2] = Construct(numOfMemElem2, g1_2, g2_2, 1) ;
end
    %figure;
    if Identical_decoders == 1
        drawExitChart(g1_1,g2_1,trellis,template, BitsBPSK, SNR, stim_em, Identical_decoders, notPuncInd,cont)  ;
    else
        drawExitChart({g1_1,g1_2},{g2_1,g2_2},{trellis1,trellis2},{template1,template2},...
            {BitsBPSK1,BitsBPSK2}, SNR, stim_em, Identical_decoders, {notPuncInd1,notPuncInd2},cont)            ;
    end
    pbaspect([1 1 1]);
    %saveas(fh,filename);
 % Create an image file. This png image is displayed on the webpage. Its size 
% must be 640x480.
 %print(filename, '-dgif','-S640,480', '-r0');
  print(fh, '-dpng','-r72', filename); 
end
%% Functions used

    function [trellis, template, Encoded, notPuncInd, toBeSent, BitsBPSK] = Construct(numOfMemElem, g1, g2, punc)
    global BitsUncoded;
    trellis        = poly2trellis( numOfMemElem + 1, [g1 g2] )             ;     %Encoding
    trellis        = rcs( trellis )                                        ;                   
    template       = constructTemplate( trellis )                          ;     %Constructing template
    Encoded        = convenc( BitsUncoded, trellis )                       ;
    if punc == 1
        [tBS,nP]   = puncture( Encoded )                                   ;     %Puncturing
    else
        tBS        = Encoded                                               ;
        nP         = 1:length( Encoded )                                   ;
    end
    BitsBPSK       = 2*tBS -1                                              ;      %BPSK Modulation
    notPuncInd     = nP                                                    ;
    toBeSent       = tBS                                                   ; 
    end
    
    function [trell]    = rcs(t)
    x                         = de2bi( t.outputs,log2(max(max(t.outputs))+1));
    y                         = reshape( x(:, 2), t.numStates, 2 )           ;
    for i = 1:length( t.nextStates(:, 1) ) %nextStates
        if y(i,1) == 1
            temp              = t.nextStates(i,1)                            ;
            t.nextStates(i,1) = t.nextStates(i,2)                            ;
            t.nextStates(i,2) = temp                                         ;
            temp              = t.outputs(i,1)                               ;
            t.outputs(i,1)    = t.outputs(i,2)                               ;
            t.outputs(i,2)    = temp                                         ;
        end
    end
            trell             = t                                            ;
    end

    function [ output]  = LOGMAP( template, trellis, BitsBPSK, Ia, noise, variance, L_a_input,  notPuncInd)
    
    global NumberOfbits
    global histLimit    
    L_c                    = 2 / variance                                                                  ;    
    received               = BitsBPSK + noise                                                              ; 
    BitsReceived           = zeros( 1, 2*NumberOfbits )                                                    ;
    nP                     = notPuncInd                                                                    ;
    BitsReceived(nP)       = received                                                                      ;
    L_a                    = L_a_input                                                                     ;
    gammas                 = constructGammas( template, trellis, L_a,BitsReceived, variance )              ;
    [Alphas,Betas]         = constructAlphasBetas( trellis, gammas )                                       ;
    Betas_LLR              = Betas  ( reshape( trellis.nextStates, 2*trellis.numStates, 1 ) + 1, :)        ;
    Alphas_LLR             = Alphas ( ( 1:size(Alphas,1))' * ones(1, 2), :)                                ;
    Gammas_LLR             = [ gammas( 1:2:length(gammas(:,1)), :); gammas( 2:2:length(gammas(:, 1)), :) ] ;
    LLR_pre                = Betas_LLR + Alphas_LLR + Gammas_LLR                                           ;
    x1                     = LLR_pre( trellis.numStates+1:2*trellis.numStates, : )                         ;
    x0                     = LLR_pre( 1:trellis.numStates, : )                                             ;
    LLR_log1               = max(x1) + log( sum( exp(x1 - repmat( max(x1), trellis.numStates, 1 ) ), 1 ) ) ;
    LLR_log0               = max(x0) + log( sum( exp(x0 - repmat( max(x0), trellis.numStates, 1 ) ), 1 ) ) ;
    LLR_log                = LLR_log1 - LLR_log0                                                           ;
    L_e                    = LLR_log - L_a  -L_c *( BitsReceived( 1:2:end ) )                              ;
    L_e(L_e >=  histLimit) = histLimit                                                                     ;
    L_e(L_e <= -histLimit) = -histLimit                                                                    ;
    Ie                     = measure_mutual_information_histogram(L_e,Ia)                                  ;
    output                 = Ie                                                                            ;    
    end
    
    function [Template] = constructTemplate(trell)
        global NumberOfbits ;   
        % Constructing Output matrix
        output_Matrix         = nan(trell.numStates, 2*trell.numStates)                                       ;
        Input_Matrix          = nan(trell.numStates, trell.numStates)                                         ;
        for i=1:trell.numStates 
            for j=1:2
                output_Matrix( i, 2*(trell.nextStates(i,j))+1:2*(trell.nextStates(i,j))+2 )...
                              = 2*de2bi(trell.outputs(i,j)', 2, 'left-msb') - 1                               ;
                Input_Matrix ( i, trell.nextStates(i,j)+1) ...
                              = find( trell.nextStates(i,:) == trell.nextStates(i,j) ) - 1                    ;
            end
        end
        trellisBranchOutputs  = 2*de2bi( trell.outputs', log2( max(max(trell.outputs)) + 1 ), 'left-msb') - 1 ;
        %Constructin previous states matrix
        prevStates            = [ ]                                                                           ;
        for i=1:trell.numStates
            prevStates        = [ prevStates; ( find( sum( (trell.nextStates == i - 1)' )' == 1 ) - 1 )' ]    ;
        end
        previous              = prevStates                                                                    ;
        for i = 2:trell.numStates-1
            temp              = previous(i, 1)                                                                ;
            previous(i,1)     = previous(i, 2)                                                                ;
            previous(i,2)     = temp                                                                          ;
        end
        [a_prev, b_prev]      = sort( reshape( previous', 2*trell.numStates, 1 ) + 1, 'ascend' )              ;
        for i = 3:2:(2*trell.numStates) - 2
            temp              = b_prev(i  , 1)                                                                ;
            b_prev(i, 1)      = b_prev(i+1, 1)                                                                ;
            b_prev(i+1, 1)    = temp                                                                          ;
        end
        % Entrance stage for the trellis
        activeStates          = 0                                                                             ;
        count                 = 0                                                                             ;
        output_Conc           = [ ]                                                                           ;
        while length( activeStates ) < trell.numStates
            output            = [ ]                                                                           ;
            tN                = trell.nextStates( activeStates+1, : )                                         ;
            for i = 1:trell.numStates
                if ismember( i, activeStates + 1 )
                    output    = [ output; 2*de2bi( trell.outputs(i, :)', log2( max( max(trell.outputs))+1)...
                                 , 'left-msb') - 1 ]                                                          ; 
                else
                    output    = [ output; nan nan; nan nan ]                                                  ;
                end
            end
            output_Conc       = [ output_Conc output]                                                         ;
            [ newStates, not]    = unique( reshape( tN, 1, numel(tN) ) )                                         ;
            for i = 1:length(newStates)
              if ~ismember( newStates(i), activeStates )
                 activeStates = [ activeStates newStates(i) ]                                                 ;
              end
            end
            count             = count + 1                                                                     ;
        end
        %Exit stage for the trellis
        activeStates          = 0                                                                             ;
        countPrev             = 0                                                                             ;
        output_Conc_Prev      = [ ]                                                                           ;
        while length(activeStates) < trell.numStates
            output            = [ ]                                                                           ;
            output_temp       = nan( 2*trell.numStates, 2 )                                                   ;
            for i = 1:trell.numStates
                if ismember( i, activeStates + 1 )
                    xx        = previous(i, :)                                                                ;
                    for k = 1:2
                        offset                           = Input_Matrix (xx(k)+1, i        )                  ;
                        output_temp(2*xx(k)+1+offset,:)  = output_Matrix(xx(k)+1, 2*i-1:2*i)                  ;
                    end
                end
            end
            output            = [ output_temp output ]                                                        ;
            output_Conc_Prev  =[output output_Conc_Prev ]                                                     ;
            tP                = prevStates(activeStates + 1 , :)                                              ;
            [newPrevStates, not] = unique( reshape(tP, 1, numel(tP) ) )                                          ;
            for i = 1:length(newPrevStates)
                if ~ismember( newPrevStates(i), activeStates)
                    activeStates = [ activeStates newPrevStates(i) ]                                          ;
                end
            end
            countPrev            = countPrev + 1                                                              ;
        end
        steady                   = repmat( trellisBranchOutputs, 1, NumberOfbits - count - countPrev )        ;    
        Template                 = [ output_Conc steady output_Conc_Prev ]                                    ;
    end

    function [gammas]   = constructGammas(template, trellis, L_a,BitsReceived,variance)
     inf                   = 1000                                                                                 ;
     L_a_Matrix            = repmat( L_a, 2*trellis.numStates, 1 )                                                ;         
     BitsReceived_Matrix   = repmat( BitsReceived, 2*trellis.numStates, 1 )                                       ;
     gammas                = 0.5*(((2/variance)*BitsReceived_Matrix(:,1:2:end)+L_a_Matrix).*template(:,1:2:end)...
                           +       (2/variance)*BitsReceived_Matrix(:,2:2:end)            .*template(:,2:2:end))  ;               
     gammas(isnan(gammas)) = -inf                                                                                 ;
    end

    function [Alphas_log,Betas_log] = constructAlphasBetas(trellis, gammas)
    global NumberOfbits ;   
    numStates         = trellis.numStates;
    inf               = 1000;
    % gammas_ratio = exp(gammas);
    [a,b]             = sort( reshape( trellis.nextStates', 2*trellis.numStates, 1 ) + 1, 'ascend') ;
    bBeta             = reshape( trellis.nextStates', 2*trellis.numStates, 1 ) +1                   ;
    gammas1           = gammas(b, :)                                                                ;
    %Alpha elements
    alphas            = zeros( 2*trellis.numStates,1)                                               ;
    alphas_final      = zeros( trellis.numStates,NumberOfbits)                                      ;
    alphas_final(:,1) = [ 0; -1000*ones( trellis.numStates-1, 1 ) ]                                 ;
    %Beta elements
    iBeta             = NumberOfbits - 1                                                            ;
    Betas             = zeros( 2*trellis.numStates, 1 )                                             ;
    Betas_final       = zeros( trellis.numStates, NumberOfbits)                                     ;
    nn                = NumberOfbits                                                                ;
    Betas_final(:,nn) = [0; -1000*ones( trellis.numStates-1, 1 ) ]                                  ;
    
    for i=2:NumberOfbits
        xx            = alphas_final(:, i-1)'                                                       ;
        alphas_temp   = reshape( xx([1:size(xx,1)]' * ones(1,2), :),2*numel(alphas_final(:,i-1)),1) ;
        alphas_temp   = alphas_temp (b,:)                                                           ;
        alphas        = (alphas_temp + gammas1(:,i-1))                                              ;
        x             = Betas_final(:,iBeta+1)                                                      ;
        Betas_temp    = x(bBeta,:)                                                                  ;
        Betas         = (Betas_temp + gammas(:,iBeta+1))                                            ;
        ii            = 1                                                                           ;
    for j = 1:trellis.numStates
             alphas_final(j,i)    = max( alphas(ii), alphas(ii+1) ) ...
                                  + log( 1 + exp( -abs( alphas(ii) - alphas(ii+1) ) ) )             ;
             Betas_final(j,iBeta) = max( Betas (ii), Betas (ii+1) ) ...
                                  + log( 1 + exp( -abs( Betas (ii) - Betas (ii+1) ) ) )             ;
             ii                   = ii + 2                                                          ;
        end
        
        if ( i <= (log2(numStates)) || i >= NumberOfbits - (log2( numStates ) ) + 1 )
            alphas_final(isnan(alphas_final))   = -inf                                              ;
            Betas_final(isnan(Betas_final))     = -inf                                              ;
        end
        iBeta                                   = iBeta - 1                                         ;
    end     
    Betas_log   = Betas_final                                                                       ;
    Alphas_log  = alphas_final                                                                      ;
    end

    function [ JI ] = JInverseFunction( I )
        I_tp     =  0.3646                                                               ;
        a_alpha1 =  1.09542                                                              ;
        b_alpha1 =  0.214217                                                             ;
        c_alpha1 =  2.33727                                                              ;
        a_alpha2 =  0.706692                                                             ;
        b_alpha2 =  0.386013                                                             ;
        c_alpha2 = -1.75017                                                              ;
        JI_temp  = zeros(1,length(I))                                                    ;
        JI_temp(I >= 0   & I <= I_tp)  = (a_alpha1*(I(I >= 0 & I <= I_tp).^2)) ...
                                       + (b_alpha1*(I(I >= 0 & I <= I_tp)))    ...
                                       + (c_alpha1*(I(I >= 0 & I <= I_tp).^(1/2)))       ;
        JI_temp(I > I_tp & I <= 1   )  = (-a_alpha2*log(-b_alpha2*((I(I > I_tp ...
                                       & I <= 1))-1)))-(c_alpha2*(I(I > I_tp & ...
                                       I <= 1 ) ) )                                      ;
        JI      = JI_temp                                                                ;
    end 
    
    function [punc, notPuncInd] = puncture(enc) 
    puncInd            = 2:4:length(enc)                   ;
    notPuncInd         = setdiff( 1:length(enc), puncInd ) ;
    punc_temp          = enc                               ;
    punc_temp(puncInd) = [ ]                               ;
    punc               = punc_temp                         ;
    end

    function [LLR] = generate_La(Ia)
        global BitsUncoded;
        global histLimit;
        sigma_a = JInverseFunction( Ia )                         ;
        mu_a    = ( sigma_a^2 ) / 2                              ;
        u       = 2*BitsUncoded - 1                              ;
        if     Ia==1  
        LLR     = u*histLimit                                    ;
        elseif Ia==0                                             ;
        LLR     =  zeros(1, length(u) )                          ;
        else
        LLR     =  normrnd( u*(mu_a), sigma_a, 1, length(u) )    ;
        end
    
    end 

    function [mutual_information] = measure_mutual_information_histogram(llrs,Ia_sent)
        global BitsUncoded;
        global histLimit;
        indices1           = BitsUncoded == 1                                                                       ;            
        indices0           = BitsUncoded == 0                                                                       ;
        Num_HistBins       = 400                                                                                    ;                       
        xAxis              = linspace( -histLimit,histLimit, Num_HistBins )                                         ;   
        hist1              = hist( llrs(indices1), xAxis )                                                          ;
        hist0              = hist( llrs(indices0), xAxis )                                                          ;
        count1             = sum( BitsUncoded )                                                                     ;
        count0             = length(BitsUncoded) - count1                                                           ;
        if count1~=0;        hist1=hist1/(count1); end                                                              ;
        if count0~=0;        hist0=hist0/(count0); end                                                              ;
        t0                 = sum( hist0(hist0~=0).*log2( (2*hist0(hist0~=0))./(hist1(hist0~=0)+hist0(hist0~=0))) )  ;   
        t1                 = sum( hist1(hist1~=0).*log2( (2*hist1(hist1~=0))./(hist1(hist1~=0)+hist0(hist1~=0))) )  ;
        mutual_information = 0.5 * ( t0 + t1 )                                                                      ;
        if Ia_sent == 1; mutual_information = 1; end                                                                ;
  end 

    function drawExitChart(g1,g2,trellis, template, BitsBPSK, SNR,stim_em,ident, notPuncInd,cont)
    global NumberOfbits ;   
    global Ia;
    %if length(SNR)>1 ; ident =1; end
    if ident ~=1
        trellis1     = trellis{1,1}                                        ;
        trellis2     = trellis{1,2}                                        ;
        template1    = template{1,1}                                       ;
        template2    = template{1,2}                                       ;
        BitsBPSK1    = BitsBPSK{1,1}                                       ;
        BitsBPSK2    = BitsBPSK{1,2}                                       ;
        notPuncInd1  = notPuncInd{1,1}                                     ;
        notPuncInd2  = notPuncInd{1,2}                                     ;
        noEnc        = 2                                                   ;
        g1_1         = g1{1,1}                                             ;
        g2_1         = g2{1,1}                                             ;
        g1_2         = g1{1,2}                                             ;
        g2_2         = g2{1,2}                                             ;    
    else
         trellis1    = trellis                                             ;
         BitsBPSK1   = BitsBPSK                                            ;
         template1   = template                                            ;
         notPuncInd1 = notPuncInd                                          ;
         trellis2    = trellis                                             ;
         BitsBPSK2   = BitsBPSK                                            ;
         template2   = template                                            ;
         notPuncInd2 = notPuncInd                                          ;
         noEnc       = 1                                                   ;
    end
        first        = 0                                                   ;
        L_a_toIa     = zeros( length(Ia), NumberOfbits )                   ;
        cstring      = 'gbrkcmbygc'; % color string
        icolor       = 1;           % color counter
        legendString = { } ;
        leg          = { } ;
        disp('Working on Exit chart')
        SNR          = sort (SNR, 'ascend');
        
        if length(SNR) == 1 % Exit chart for a single SNR
           Ie        = zeros(1,length(Ia))                                 ;
           Ie2       = zeros(1,length(Ia))                                 ;
           SNR_ratio = 10^(SNR / 10)                                       ;
           Eb        = 1                                                   ;
           N0        = Eb / SNR_ratio                                      ;
           variance  = N0                                                  ;
           noise     = normrnd( 0, sqrt(variance), 1, length(BitsBPSK1) )  ;
           disp ( ['SNR =', num2str(SNR)] )                                ;
           for iia = 1:length(Ia)
               L_a_toIa(iia,:) = generate_La(Ia(iia))                      ;
               Ie(iia)         = LOGMAP(template1,trellis1, BitsBPSK1,...
                   Ia(iia), noise, variance, L_a_toIa(iia,:), notPuncInd1) ;
           end
           if ident == 1
                Ie2 = smooth(Ie)                                           ;
                Ia2 = Ia                                                   ;
           else
              for iia = 1:length(Ia)
                 Ie2(iia)     = LOGMAP(template2,trellis2, BitsBPSK2,...
                    Ia(iia), noise, variance,L_a_toIa(iia,:), notPuncInd2) ;
              end
                Ie2  = smooth(Ie2)                                         ;
                Ia2  = Ia                                                  ;
            end
        if ~isrow( Ie2 ) ; Ie2 = Ie2' ; end
            
              plot(Ia,smooth(Ie),cstring(icolor),'linewidth',2)                    ;
              hold on
              leg = [ leg, [ '1^s^t decoder' ] ]                   ; 
              plot( Ie2, Ia2, [ cstring(icolor+1), '-x' ],'markersize', 1,'linewidth',2 )                ; 
              hold on;
              leg = [ leg, [ '2^n^d decoder']]            
              if ( trapz( Ia, smooth(Ie)' ) <  trapz( Ie2, Ia2 ) ) ...
                          || ~all( smooth(Ie)' >= Ia )                        %not Converging
                          contf =0;
                   disp( ['not converging for SNR = ', num2str(SNR),'dB']) ;
                   drawStairSteps_emulated( noEnc, Ia, {Ie,Ie2}, SNR)      ;
                   leg = [ leg, 'not converging' ]                         ;

                  else % converging
                  contf = 1;
                  ebno = SNR_ratio;
                     if stim_em == 1
                        it_sim = drawStairSteps_simulated( trellis,     ...
                            template, BitsBPSK, variance, Ia(1)   ,     ...
                            noise, notPuncInd )                            ;
                        leg    = [leg, ['simulated in ', num2str(it_sim)...
                            , ' iterations' ] ]                            ;
                    elseif stim_em == 0
                        it_em  = drawStairSteps_emulated ( noEnc, Ia,   ...
                            { Ie, Ie2 }, SNR )                             ;
                        leg    = [leg, [ 'emulated in ', num2str(it_em) ...
                            ' iterations' ] ]                              ;
                    elseif stim_em == 2
                        it_sim = drawStairSteps_simulated( trellis,     ...
                            template, BitsBPSK, variance, Ia(1)   ,     ...
                            noise, notPuncInd)                             ;
                        hold on
                        it_em  = drawStairSteps_emulated ( noEnc, Ia,   ...
                            { Ie, Ie2 }, SNR )                             ;
                        leg    ={ '1^s^t decoder', '2^n^d decoder',     ...
                            [ 'simulated in ', num2str( it_sim )  ,     ...
                            ' iterations' ], [ 'emulated in '     ,     ...
                            num2str(it_em) ' iterations' ] }               ;
                    else
                        error('Only 1,0 and 2 are the only valid entries') ;
                    end
                  end
              legendString   = leg                                         ;
        else %more than 1 SNR 
        for  iSNR = 1:length(SNR)
            leg       = { }                                                ;
            Ie        = zeros(1,length(Ia))                                ;
            Ie2       = zeros(1,length(Ia))                                ;
            SNR_ratio = 10^(SNR(iSNR) / 10)                                ;
            Eb        = 1                                                  ;
            N0        = Eb / SNR_ratio                                     ;
            variance  = N0                                                 ;
            noise     = normrnd( 0, sqrt(variance), 1, length(BitsBPSK1) ) ;
            disp( [ 'SNR =', num2str( SNR(iSNR) ) ] )                      ;
                
            for iia = 1:length(Ia)
                if iSNR == 1
                L_a_toIa(iia,:) = generate_La( Ia(iia) )                   ;
                end
                Ie(iia)         = LOGMAP(template1,trellis1, BitsBPSK1, ...
                    Ia(iia), noise, variance,L_a_toIa(iia,:), notPuncInd1) ;
            end
                if ident
                Ie2             = smooth(Ie)                               ;
                Ia2             = Ia                                       ;
            else
            	 for iia = 1:length(Ia)
                 Ie2(iia)     = LOGMAP(template2,trellis2, BitsBPSK2,...
                    Ia(iia), noise, variance,L_a_toIa(iia,:), notPuncInd2) ;
              end
                Ie2  = smooth(Ie2)                                         ;
                Ia2  = Ia                                                  ;
            end               
            if ~isrow( Ie2 ) ; Ie2 = Ie2' ; end
            
               
            if ( trapz( Ia, smooth(Ie)' ) <  trapz(Ie2, Ia2)) ...
                    || ~all( smooth(Ie)' >= Ia )
                disp( [ 'not converging for SNR = ', num2str(SNR(iSNR)) ...
                    ' dB' ] )                                              ;
                leg            = { '' }                                    ;
                contf = 0;
             else %converging          
                plot (Ia,smooth(Ie),cstring(icolor),'linewidth',2)
                leg            = {['SNR = ' num2str(SNR(iSNR)) ' dB' ]};
                hold on
                if first == 0
                contf = 1;
                ebno = SNR_ratio;
                    first      = 1                                         ;
                    plot( Ie2, Ia2, [ cstring(icolor), '--' ],'linewidth',2 )            ;
                    leg        = [ leg, [ '2^n^d (Iterative) decoder for'...
                        , ' ', num2str( SNR(iSNR) ), ' dB' ] ]             ;
                    hold on
                    drawStairSteps_emulated(noEnc, Ia, {Ie,Ie2}, SNR(iSNR));
                    leg        = [ leg, [ 'iterations for ', num2str(   ...
                        SNR(iSNR)), ' dB' ] ]                             ;
                    hold on
                end
                icolor         = icolor + 1                                ; 
            end
            
                legendString   = [ legendString, leg{ : } ]                ;
        end
        
       

        end
         ylabel ('I_e')
         xlabel ('I_a')
         axis   ([0 1 0 1]);
         if ident ==1
               title ('Exit chart for the given decoders' )      ;
         end
         legend (legendString,'Location','SouthEast');
         if cont && contf
               ie=0:0.005:1;
               
               index = 1;
               for ia = 0:0.005:1
                 pb(:,index) = 0.5*erfc((sqrt((4*ebno)+(JInverseFunction(ie).^2)+(JInverseFunction(ia).^2)))/(2*sqrt(2)));
                 index = index + 1;
                    end
               hold on;
               [C,h] = contour(0:0.005:1,0:0.005:1,pb,[0.1 0.08 0.06 0.04 0.03 0.02 0.01 0.005 0.001]);
            clabel(C,h,'FontSize',10,'Color','k','Rotation',0,'LabelSpacing',700);
               end
    end
    
    function [iteration] = drawStairSteps_simulated(trellis, template, BitsBPSK, variance,Ia_1,noise, notPuncInd)
        Ia1         = Ia_1            ;
        Ie1         = 1               ;
        Ie2         = 0               ;
        xaxis_stair = [ ]             ;
        yaxis_stair = [ ]             ;
        recordx     = [ ]             ;
        recordy     = [ ]             ;
        iteration   = 0               ;
        failed      = 0               ;
        if iscell(trellis)
        trellis1    = trellis{1,1}    ;
        trellis2    = trellis{1,2}    ;
        template1   = template{1,1}   ;
        template2   = template{1,2}   ;
        BitsBPSK1   = BitsBPSK{1,1}   ;
        BitsBPSK2   = BitsBPSK{1,2}   ;
        notPuncInd1 = notPuncInd{1,1} ;
        notPuncInd2 = notPuncInd{1,2} ;
        else
        trellis1    = trellis         ;
        trellis2    = trellis         ;
        template1   = template        ;
        template2   = template        ;
        BitsBPSK1   = BitsBPSK        ;
        BitsBPSK2   = BitsBPSK        ;
        notPuncInd1 = notPuncInd      ;
        notPuncInd2 = notPuncInd      ;
        end
            
        while~( abs( Ie1 - Ie2 ) < 0.001 && ( Ie1 > 0.99 || Ie2 > 0.99 ) )
            Ie1 = LOGMAP(template1, trellis1, BitsBPSK1, Ia1,noise, variance, generate_La(Ia1), notPuncInd1)    ;
            Ia2 = Ie1                                                                                           ;
            Ie2 = LOGMAP(template2, trellis2, BitsBPSK2, Ia2,noise, variance, generate_La(Ia2), notPuncInd2)    ;
            if Ie2 > Ia1
            iteration   = iteration+1
            xaxis_stair = [ xaxis_stair Ia1 Ie2]    ;
            yaxis_stair = [ yaxis_stair Ia2 Ie1]    ;
            Ia1         = Ie2                       ;
            else
                failed  = failed + 1                ;
                recordx = [recordx Ia1 Ie2]         ;
                recordy = [recordy Ia2 Ie1]         ;
                if failed > 5
                 break                              ;
                end
            end
            
        end
        if ( Ie1~=1 || Ie2~=1 ) && ( abs(Ie1-Ie2)<0.001 && ( Ie1>0.99 || Ie2>0.99 ) )
                xaxis_stair = [ xaxis_stair xaxis_stair(end) 1]    ;
                yaxis_stair = [ yaxis_stair 1 1]                   ;
                iteration   = iteration+1
        end
        plot([0 xaxis_stair],[0 yaxis_stair])                      ;
       
    end

    function [iteration] = drawStairSteps_emulated (noOfEnc,Ia,Ie,SNR)
    if noOfEnc == 2
        IeEnc1 = Ie{1,1}     ;
        IeEnc2 = Ie{1,2}     ;
    else
        IeEnc1 = Ie{1,1}	 ;
        IeEnc2 = IeEnc1      ;
    end
    Ia_interp  = 0:0.001:1   ;
    
        Ie_interp                = interp1(Ia,IeEnc1,Ia_interp,'spline')    ;
        Ie2_interp               = interp1(Ia,IeEnc2,Ia_interp,'spline')    ;
        Ie_interp (Ie_interp>1 ) = 1                                        ;
        Ie2_interp(Ie2_interp>1) = 1                                        ;
        
        iteration                = 0                                        ;
        step                     = Ia_interp(2) - Ia_interp(1)              ;
        Ie1                      = Ie_interp(1)                             ;
        Ie2                      = 0                                        ; %initially random
        Ia1                      = Ia_interp(1)                             ;
        staircase_x              = []                                       ;
        staircase_y              = []                                       ;
        failed                   = 0                                        ;
        while ~(Ie1==1 && Ie2~=1)
            if (failed>5 & (Ie1>=0.99 || Ie2>=0.99))
                iteration        = iteration+1
                staircase_x      = [ staircase_x staircase_x(end) 1]        ;
                staircase_y      = [ staircase_y 1 1]                       ;
                break                                                       ;
            elseif failed>5
                break                                                       ;
                
            end
            if  Ie1 == round(Ie1/step)*step
                if Ie1 + step >= 1
                    Ie1 = 1                                      ;
                else
                    Ie1 = mod((Ie1 + 10*step),1)                 ;
                end
            end
            Ie1         = round(Ie1/step)*step                   ;
            Ia2         = Ie1                                    ; 
            Ie2         = Ie2_interp(abs(Ia2-Ia_interp)<1e-10)   ;
            if  Ie2 == round(Ie2/step)*step
                if Ie2 + step >= 1
                    Ie2 = 1                                      ;
                else
                    Ie2 = mod((Ie2 +step),1)                     ;
                end
            end
            Ie2 = round(Ie2/step)*step;
            if Ie2 > Ia1
                staircase_x = [staircase_x Ia1 Ie2]                   ;
                staircase_y = [staircase_y Ie1 Ia2]                   ;
                Ia1         = Ie2                                     ;
                Ie1         = Ie_interp(abs(Ia_interp-Ia1)<1e-10)     ;
                if Ie1 == 1 & staircase_y(end)~=1
                    Ie1   = Ie1 - step                                ;
                end
                iteration = iteration+1
            else
                failed    = failed+1                                  ;
            end
        end
    plot([0 staircase_x],[0 staircase_y],'r')
    
    end