function T = transfer_matrix_1d(E_range, V0, V1, lambda_M, N_periods, mstar, hbar)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  TRANSFER_MATRIX_1D  1D transfer-matrix tunneling through moire barriers
%
%  Computes the transmission coefficient T(E) for a quasiparticle
%  tunneling through N_periods of the moire superlattice potential.
%
%  The potential along a 1D slice (through AA sites) is:
%    V(x) = V0 * cos(2*pi*x/lambda_M) + V1 * cos(4*pi*x/lambda_M)
%
%  Method:
%    Divide each period into N_slabs thin slabs of width dx.
%    For each slab, compute the local wavevector:
%      - Propagating: k_j = sqrt(2*m*(E - V_j)) / hbar
%      - Evanescent:  kappa_j = sqrt(2*m*(V_j - E)) / hbar
%    Build the 2x2 transfer matrix for each slab and multiply.
%    Transmission: T = 1 / |M_22|^2
%
%  Inputs:
%    E_range   - array of incident energies [J]
%    V0        - barrier height (first harmonic) [J]
%    V1        - second harmonic amplitude [J]
%    lambda_M  - moire period [m]
%    N_periods - number of moire periods to tunnel through
%    mstar     - effective mass [kg]
%    hbar      - reduced Planck constant [J.s]
%
%  Output:
%    T         - transmission coefficient array, same size as E_range
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N_slabs_per_period = 200;  % resolution per period
N_total = N_periods * N_slabs_per_period;
L_total = N_periods * lambda_M;
dx = L_total / N_total;

% Position array for the barrier region
x = linspace(0, L_total, N_total);

% Moire potential along the 1D slice
V_x = V0 * cos(2*pi*x/lambda_M) + V1 * cos(4*pi*x/lambda_M);
% Shift so minimum is at 0 (barrier measured from ground)
V_x = V_x - min(V_x);

T = zeros(size(E_range));

for iE = 1:length(E_range)
    E = E_range(iE);

    % Initialize total transfer matrix as identity
    M = eye(2);

    for j = 1:N_total
        dV = E - V_x(j);

        if dV > 0
            % Propagating region
            k_j = sqrt(2 * mstar * dV) / hbar;
            phi = k_j * dx;
            % Propagation matrix
            Mj = [cos(phi),         sin(phi)/k_j;
                  -k_j*sin(phi),    cos(phi)];
        else
            % Evanescent region (barrier)
            kappa_j = sqrt(2 * mstar * abs(dV)) / hbar;
            phi = kappa_j * dx;
            % Evanescent matrix
            Mj = [cosh(phi),           sinh(phi)/kappa_j;
                  kappa_j*sinh(phi),   cosh(phi)];
        end

        M = Mj * M;

        % Numerical stabilization: if matrix elements get too large,
        % we know T will be effectively zero
        if abs(M(2,2)) > 1e30
            break;
        end
    end

    % Transmission coefficient
    % T = |t|^2, and for symmetric potential: t = 1/M_22
    % More precisely, for asymmetric: T = 1/|M_22|^2
    % (assuming free-particle states on both sides with same mass)
    T(iE) = 1 / (abs(M(2,2))^2 + 1e-50);
    T(iE) = min(T(iE), 1.0);  % cap at unity
end

end