#!/usr/bin/env python3

import matplotlib.pyplot as plt
import scipy.linalg as la
import numpy as np

def Main():

    h = 1.015

    x = 2**np.arange(4, 11)

    energy_gap = []
    S = []

    print('L\tD')
    for L in x:

        E, ... = ModelSolution(L, h)

        energy_gap.append( ... )
        print(f'{L}\t{energy_gap[-1]}')

    plt.loglog(1/x, energy_gap, marker='.', linestyle='')
    plt.xlabel('$1/L$')
    plt.ylabel('$\\delta$')
    plt.show()
    plt.close()

    print('l\tS')
    l = x[:-1] # np.arange(1, L//2+1)
    for u in l:
        S.append( Entropy(CorrMatrix(..., u)) )
        print(f'{u}\t{S[-1]}')

    X = np.log(...)/6

    plt.plot(l, S, marker='.', linestyle='')
    plt.xlabel('$l$')
    plt.ylabel('$S$')
    plt.show()
    plt.close()

    return

def ModelSolution(L, h):
    A = np.zeros((L, L))
    B = np.zeros((L, L))

    for l in range(...):
        B[...] = 
        A[...] = 
        ...

    M = np.zeros((2*L, 2*L))

    M[:L,:L] = A
    M[L:,L:] = -A
    M[:L,L:] = B
    M[L:,:L] = -B

    E, U = la.eigh(M, subset_by_index=[L, 2*L-1])
    # E is the vector of positive eigenvalues from smallest to largest
    # U is the matrix such that U^T * M * U^T = diag(E)

    U1 = U[:L,:L]
    U2 = U[:L,L:]
    U3 = U[L:,:L]
    U4 = U[L:,L:]

    return E, ...

def CorrMatrix(U1, ..., u):

    ### Truncated matrices ###
    T1 = U1[...]
    ...

    M = np.zeros((2*T1.shape[0], 2*T1.shape[0]))
    M[:u,:u] = ...
    M[u:,u:] = ...
    M[:u,u:] = ...
    M[u:,:u] = ...
    ### np.dot(a, b) for vector/matrix products

    return M

def Entropy(C):
    Lambda = la.eigh(C, eigvals_only=True, subset_by_value=[1e-4, 1-1e-4])
    return -np.sum( Lambda*np.log(Lambda) )

Main()