function [V_moire, V_type] = compute_moire_potential(X, Y, lambda_M, V0, V1, theta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  COMPUTE_MOIRE_POTENTIAL  Moire superlattice potential for twisted bilayer
%
%  The moire potential arises from the spatial modulation of interlayer
%  stacking. In a twisted TMD bilayer, the three high-symmetry points
%  (AA, AB, BA) create a triangular superlattice with period lambda_M.
%
%  V_moire(r) = V0 * sum_{j=1}^{3} cos(G_j . r)
%             + V1 * sum_{j=1}^{3} cos(2 * G_j . r)
%
%  where G_j are the three moire reciprocal lattice vectors at 60 deg
%  intervals, and V0, V1 are the first and second harmonic amplitudes.
%
%  Inputs:
%    X, Y      - meshgrid coordinate arrays [m]
%    lambda_M  - moire period [m]
%    V0        - first harmonic amplitude [J]
%    V1        - second harmonic amplitude [J]
%    theta     - twist angle [radians]
%
%  Outputs:
%    V_moire   - potential energy landscape [J] (same size as X)
%    V_type    - stacking type indicator: 0=AA, 0.5=SP, 1=AB/BA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Moire reciprocal lattice vectors
% G_M = 4*pi / (sqrt(3) * lambda_M)
G_M = 4*pi / (sqrt(3) * lambda_M);

% Three G vectors at 0, 60, 120 degrees
G = zeros(3, 2);
for j = 1:3
    ang = (j-1) * 2*pi/3;
    G(j, :) = G_M * [cos(ang + pi/6), sin(ang + pi/6)];
end

% First harmonic: honeycomb potential
V_moire = zeros(size(X));
for j = 1:3
    phase = G(j,1)*X + G(j,2)*Y;
    V_moire = V_moire + cos(phase);
end
V_moire = V0 * V_moire / 3;  % normalize so max ~ V0

% Second harmonic: sharper features at AA sites
if V1 ~= 0
    V2 = zeros(size(X));
    for j = 1:3
        phase2 = 2*(G(j,1)*X + G(j,2)*Y);
        V2 = V2 + cos(phase2);
    end
    V_moire = V_moire + V1 * V2 / 3;
end

% Stacking type indicator (0 = AA-aligned, 1 = AB/BA)
% At AA sites, all three cosines = +1, giving V ~ +V0 (maximum)
% The potential minimum (AB/BA) is the tunneling-favorable site
% We invert so AA nodes (Josephson nodes) are valleys
V_moire = -V_moire;

% Compute stacking-type map
% AA sites: where |sum cos(G.r)| ~ 3 (all aligned)
stacking_sum = zeros(size(X));
for j = 1:3
    stacking_sum = stacking_sum + cos(G(j,1)*X + G(j,2)*Y);
end
V_type = 1 - (stacking_sum + 3) / 6;  % 0 at AA, 1 at AB/BA

end