import numpy as np

# Define Pauli matrices and identity matrix
X = np.array([[0, 1], [1, 0]], dtype=np.complex128)
Y = np.array([[0, -1j], [1j, 0]], dtype=np.complex128)
Z = np.array([[1, 0], [0, -1]], dtype=np.complex128)
I = np.eye(2, dtype=np.complex128)

def tensor_product(ops):
    result = ops[0]
    for op in ops[1:]:
        result = np.kron(result, op)
    return result

def construct_transformed_hamiltonian_OBC(L, alpha, m0, eta):
    dim = 2 ** L  # Total dimension for a system of L qubits
    kinetic_term = np.zeros((dim, dim), dtype=np.complex128)
    mass_term = np.zeros((dim, dim), dtype=np.complex128)
    interaction_term = np.zeros((dim, dim), dtype=np.complex128)

    q = [(-1)**(k+1) for k in range(L)]  # Background charge distribution
    print(q)

    # Kinetic Term
    for k in range(1, L):
        Z_string_ops = [Z if i < k else I for i in range(L)]
        Z_string = tensor_product(Z_string_ops)
        kinetic_term += alpha / 2 * Z_string * q[k-1]  # q index adjusted to k-1

    # Mass Term
    # Z_{[1,L-1]} for the boundary condition term
    Z_boundary_ops = [Z if i < L-1 else I for i in range(L)]
    Z_boundary = tensor_product(Z_boundary_ops)
    mass_term += m0 / 2 * (-1)**L * q[L-1] * Z_boundary  # Apply (-1)^L and q_{[1,L]}

    # Sum term (-1)^k \tilde{Z}_k over k from 1 to L-1
    for k in range(1, L-1):
        Z_tilde_ops = [Z if i == k else I for i in range(L)]
        Z_tilde = tensor_product(Z_tilde_ops)
        mass_term -= m0 / 2 * (-1)**k * Z_tilde  # Summing the alternating sign mass terms


    # Interaction Term
    interaction_term = np.zeros((dim, dim), dtype=np.complex128)
    for k in range(1, L):
        X_tilde_ops = [X if i == k else I for i in range(1, L+1)]
        Y_tilde_ops = [Y if i == k else I for i in range(1, L+1)]

        if k < L:
            X_tilde_next_ops = [X if i == k + 1 else I for i in range(1, L+1)]
            Y_tilde_next_ops = [Y if i == k + 1 else I for i in range(1, L+1)]

            X_tilde_k = tensor_product(X_tilde_ops)
            Y_tilde_k = tensor_product(Y_tilde_ops)
            X_tilde_k_plus_1 = tensor_product(X_tilde_next_ops)
            Y_tilde_k_plus_1 = tensor_product(Y_tilde_next_ops)

            interaction_term += eta / 4 * (X_tilde_k.dot(X_tilde_k_plus_1) + Y_tilde_k.dot(Y_tilde_k_plus_1))


    # Combine all terms to form the full Hamiltonian
    H = kinetic_term + mass_term + interaction_term
    return H