% function P = rll_coil_P(lmin, lmax, a)
% RLL_COIL_P Markov chain for RLL(lmin,lmax) with NRZI-like polarity flips on '1'.
% States: (r,+) and (r,-), r = 0..lmax.
% Column-stochastic: P(i,j)=Pr{next=i | prev=j}.
%
% Inputs
%   lmin, lmax : integers with 0 <= lmin <= lmax
%   a          : probabilities of writing '1' in the "free" region r=lmin..lmax-1.
%                - If scalar: same probability used for all eligible r.
%                - If vector: length must be (lmax - lmin), where
%                    a(1) corresponds to r=lmin,
%                    a(2) corresponds to r=lmin+1, ..., up to r=lmax-1.
%
% Example:
%   P = rll_coil_P(1,3,[a1 a2]);  % produces (1,3) code
function P = rll_coil_P(lmin, lmax, a)
    if nargin < 2
        error('Need at least lmin and lmax.');
    end
    if lmin < 0 || lmax < 0 || lmin > lmax || floor(lmin)~=lmin || floor(lmax)~=lmax
        error('Require integers 0 <= lmin <= lmax.');
    end

    % Number of "free" r values where a probability is needed:
    nfree = max(lmax - lmin, 0);  % r = lmin..lmax-1
    if nargin < 3 || isempty(a)
        if nfree > 0
            a = 0.5 * ones(nfree,1);
        else
            a = []; % no free region
        end
    end

    if nfree > 0
        if isscalar(a)
            a = repmat(a, nfree, 1);
        else
            a = a(:);
            if numel(a) ~= nfree
                error('a must be scalar or length (lmax-lmin) = %d.', nfree);
            end
        end
        if any(a < 0) || any(a > 1)
            error('Probabilities in a must be in [0,1].');
        end
    end

    % Total states = 2*(lmax+1)
    N = 2*(lmax+1);
    P = zeros(N,N);

    % index helper: r=0..lmax, pol=+1/-1
    % ordering: (0,+),(1,+),...,(lmax,+),(0,-),(1,-),...,(lmax,-)
    idx = @(r,pol) (pol==+1)*0 + (pol==-1)*(lmax+1) + (r+1);

    for pol = [+1, -1]
        for r = 0:lmax
            j = idx(r, pol);

            if r < lmin
                % must write 0: r -> r+1, pol unchanged
                i = idx(r+1, pol);
                P(i,j) = 1;

            elseif r == lmax
                % must write 1: r -> 0, pol flips
                i = idx(0, -pol);
                P(i,j) = 1;

            else
                % free region: lmin <= r < lmax
                p1 = a(r - lmin + 1);   % prob write 1
                p0 = 1 - p1;            % prob write 0

                % write 1: reset r=0, flip polarity
                i1 = idx(0, -pol);

                % write 0: r -> r+1, polarity unchanged
                i0 = idx(r+1, pol);

                P(i1,j) = p1;
                P(i0,j) = p0;
            end
        end
    end
end
