diff --git a/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/2D_bilayer_model.py b/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/2D_bilayer_model.py
new file mode 100644
index 0000000..78a066b
--- /dev/null
+++ b/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/2D_bilayer_model.py	
@@ -0,0 +1,94 @@
+import numpy as np
+from math import *
+import cmath
+import functools
+import guan 
+
+# Installation: pip install --upgrade guan
+
+def main():
+    M = 0
+    t1 = 1
+    phi= pi/2
+    t2 = 0.03
+    tc = 0.31
+    bilayer_bands(M, t1, t2, phi, tc)
+    chern_number(M, t1, t2, phi)
+
+def hamiltonian_of_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3)): 
+    h0 = np.zeros((2, 2), dtype=complex)
+    h1 = np.zeros((2, 2), dtype=complex)
+    h2 = np.zeros((2, 2), dtype=complex)
+
+    h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
+    h1[0, 1] = h1[0, 1].conj()
+
+    h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
+    h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
+    matrix = h0 + h1 + h2 + h2.transpose().conj()
+    return matrix
+
+def hamiltonian_of_modified_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3)): 
+    h0 = np.zeros((2, 2), dtype=complex)
+    h1 = np.zeros((2, 2), dtype=complex)
+    h2 = np.zeros((2, 2), dtype=complex)
+
+    k1 = -k1  # kx  # Note that to get the unit cell the bilayer honeycomb lattice containing all possible layer hoppings, here one layer (HM layer or mHM layer） needs to be applied with its counterpart of mirror symmetry along x direction.
+
+    h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
+    h1[0, 1] = h1[1, 0].conj()
+
+    h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
+    h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
+
+    matrix = h0 + h1 + h2 + h2.transpose().conj()
+    return matrix
+
+def hamiltonian_of_bilayer(k1, k2, M, t1, t2, phi, tc):
+    hamiltonian1 = hamiltonian_of_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3))
+    hamiltonian2 = hamiltonian_of_modified_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3))
+    hamiltonian = np.zeros((4, 4), dtype=complex)
+    hamiltonian[0:2, 0:2] = hamiltonian1
+    hamiltonian[2:4, 2:4] = hamiltonian2
+
+    # AB stacking
+    hamiltonian[1, 2] = tc   # B on HM， A on mHM
+    hamiltonian[2, 1] = tc
+
+    # BA stacking
+    # hamiltonian[0, 3] = tc  # A on HM， B on mHM
+    # hamiltonian[3, 0] = tc
+
+    # AA stacking
+    # hamiltonian[0, 2] = tc
+    # hamiltonian[2, 0] = tc
+    # hamiltonian[1, 3] = tc
+    # hamiltonian[3, 1] = tc
+    return hamiltonian
+
+def bilayer_bands(M, t1, t2, phi, tc):
+    k1_array = np.linspace(0, 4*pi, 3000)
+    k2 = 0
+    hamiltonian = functools.partial(hamiltonian_of_bilayer, k2=k2, M=M, t1=t1, t2=t2, phi=phi, tc=tc)
+    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k1_array, hamiltonian)
+    plt, fig, ax = guan.import_plt_and_start_fig_ax(labelsize=25)
+    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 0],  linewidth=2)
+    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 1],  linewidth=2)
+    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 2],  linewidth=2, color='k')
+    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 3],  linewidth=2, color='k')
+    guan.plot_without_starting_fig_ax(plt, fig, ax, [], [], fontsize=30, y_max=1.5, y_min=-1.5, xlabel='$\mathrm{k_x}$', ylabel='$\mathrm{E}$')
+    plt.show()
+
+def chern_number(M, t1, t2, phi):
+    tc_array = np.arange(0.05, 0.8, 0.05)
+    chern_number_array = []
+    for tc in tc_array:
+        print(tc)
+        hamiltonian = functools.partial(hamiltonian_of_bilayer, M=M, t1=t1, t2=t2, phi=phi, tc=tc)
+        chern_number = guan.calculate_chern_number_for_honeycomb_lattice(hamiltonian, a=1/sqrt(3), precision=500, print_show=0)
+        print(chern_number, '\n')
+        chern_number_array.append(chern_number[0:2])
+    guan.plot(tc_array, np.real(chern_number_array), xlabel='$\mathrm{t_c}$', ylabel='$\mathrm{C}$', style='o-')
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file
diff --git a/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/bilayer_ribbon_model.py b/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/bilayer_ribbon_model.py
new file mode 100644
index 0000000..43f96a6
--- /dev/null
+++ b/2024.10.21_Topological hidden phase transition in honeycomb bilayers with a high Chern number/bilayer_ribbon_model.py	
@@ -0,0 +1,220 @@
+import numpy as np
+import functools
+from math import *
+import cmath
+import guan
+
+# Installation: pip install --upgrade guan
+
+def main():
+    N = 15
+    M = 0
+    t1 = 1
+    phi= pi/2
+    t2= 0.03
+    tc = 0.8
+    layer_num = 2
+    bands_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
+    conductance_of_bilayer_model_with_disorder_array(N, M, t1, t2, phi, tc, layer_num)
+
+# Haldane model
+
+def hamiltonian_of_Haldane_model(k, N, M, t1, t2, phi, sign):
+    h00 = h00_of_Haldane_model(N, M, t1, t2, phi, sign)
+    h01 = h01_of_Haldane_model(N, M, t1, t2, phi, sign)
+    hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)
+    return hamiltonian
+
+def h00_of_Haldane_model(N, M, t1, t2, phi, sign):
+    h00 = np.zeros((4*N, 4*N), dtype=complex)
+    for i0 in range(N):
+        h00[i0*4+0, i0*4+0] = M
+        h00[i0*4+1, i0*4+1] = -M
+        h00[i0*4+2, i0*4+2] = M
+        h00[i0*4+3, i0*4+3] = -M
+        h00[i0*4+0, i0*4+1] = t1
+        h00[i0*4+1, i0*4+0] = t1
+        h00[i0*4+1, i0*4+2] = t1
+        h00[i0*4+2, i0*4+1] = t1
+        h00[i0*4+2, i0*4+3] = t1
+        h00[i0*4+3, i0*4+2] = t1
+        h00[i0*4+0, i0*4+2] = t2*cmath.exp(1j*phi*sign)
+        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
+        h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
+        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
+    for i0 in range(N-1):
+        h00[i0*4+3, (i0+1)*4+0] = t1
+        h00[(i0+1)*4+0, i0*4+3] = t1
+        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
+        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
+        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
+        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
+    return h00
+
+def h01_of_Haldane_model(N, M, t1, t2, phi, sign):
+    h01 = np.zeros((4*N, 4*N), dtype=complex)
+    for i0 in range(N):
+        h01[i0*4+1, i0*4+0] = t1
+        h01[i0*4+2, i0*4+3] = t1
+        h01[i0*4+0, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
+        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+2, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
+        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
+        h01[i0*4+2, i0*4+0] = t2*cmath.exp(1j*phi*sign)
+        if i0 != 0:
+            h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
+    for i0 in range(N-1):
+        h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
+    return h01
+
+# modified_Haldane_model
+
+def hamiltonian_of_modified_Haldane_model(k, N, M, t1, t2, phi, sign):
+    h00 = h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign)
+    h01 = h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign)
+    hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)
+    return hamiltonian
+
+def h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
+    h00 = np.zeros((4*N, 4*N), dtype=complex)
+    for i0 in range(N):
+        h00[i0*4+0, i0*4+0] = M
+        h00[i0*4+1, i0*4+1] = -M
+        h00[i0*4+2, i0*4+2] = M
+        h00[i0*4+3, i0*4+3] = -M
+        h00[i0*4+0, i0*4+1] = t1
+        h00[i0*4+1, i0*4+0] = t1
+        h00[i0*4+1, i0*4+2] = t1
+        h00[i0*4+2, i0*4+1] = t1
+        h00[i0*4+2, i0*4+3] = t1
+        h00[i0*4+3, i0*4+2] = t1
+        h00[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
+        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
+        h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
+        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
+    for i0 in range(N-1):
+        h00[i0*4+3, (i0+1)*4+0] = t1
+        h00[(i0+1)*4+0, i0*4+3] = t1
+        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
+        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
+        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
+        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
+    return h00
+
+def h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
+    h01 = np.zeros((4*N, 4*N), dtype=complex)
+    for i0 in range(N):
+        h01[i0*4+1, i0*4+0] = t1
+        h01[i0*4+2, i0*4+3] = t1
+        h01[i0*4+0, i0*4+0] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+2, i0*4+2] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
+        h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
+        h01[i0*4+2, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
+        if i0 != 0:
+            h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
+    for i0 in range(N-1):
+        h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
+    return h01
+
+# bilayer model
+
+def hamiltonian_of_bilayer_model(k, N, M, t1, t2, phi, tc, layer_num):
+    modified_Haldane_model = hamiltonian_of_modified_Haldane_model(k=k, N=N, M=M, t1=t1, t2=t2, phi=phi, sign=1)
+    Haldane_model = hamiltonian_of_Haldane_model(k=k, N=N, M=M, t1=t1, t2=t2, phi=phi, sign=1)
+    hamiltonian = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex)
+    for layer in range(layer_num):
+        if np.mod(layer,2) == 0:
+            hamiltonian[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = modified_Haldane_model
+        if np.mod(layer,2) == 1:
+            hamiltonian[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = Haldane_model 
+    for layer in range(layer_num-1):
+        for i0 in range(N):
+            # AB stacking
+            hamiltonian[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+1] = tc
+            hamiltonian[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+0] = tc
+            hamiltonian[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+3] = tc
+            hamiltonian[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+2] = tc
+            
+            # # # BA stacking
+            # hamiltonian[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+0] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+1] = tc
+            # hamiltonian[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+2] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+3] = tc
+
+            # # AA stacking
+            # hamiltonian[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+0] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+0] = tc
+            # hamiltonian[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+1] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+1] = tc
+            # hamiltonian[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+2] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+2] = tc
+            # hamiltonian[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+3] = tc
+            # hamiltonian[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+3] = tc
+    return hamiltonian
+
+def h00_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
+    h00_modified = h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign=1)
+    h00_Haldane = h00_of_Haldane_model(N, M, t1, t2, phi, sign=1)
+    h00  = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex) 
+    for layer in range(layer_num):
+        if np.mod(layer,2) == 0:
+            h00[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h00_modified
+        if np.mod(layer,2) == 1:
+            h00[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h00_Haldane
+    for layer in range(layer_num-1):
+        for i0 in range(N):
+            # AB stacking
+            h00[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+1] = tc
+            h00[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+0] = tc
+            h00[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+3] = tc
+            h00[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+2] = tc
+
+            # # BA stacking
+            # h00[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+0] = tc
+            # h00[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+1] = tc
+            # h00[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+2] = tc
+            # h00[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+3] = tc
+
+            # # AA stacking
+            # h00[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+0] = tc
+            # h00[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+0] = tc
+            # h00[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+1] = tc
+            # h00[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+1] = tc
+            # h00[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+2] = tc
+            # h00[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+2] = tc
+            # h00[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+3] = tc
+            # h00[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+3] = tc
+    return h00
+
+def h01_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
+    h01_modified = h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign=1)
+    h01_Haldane = h01_of_Haldane_model(N, M, t1, t2, phi, sign=1)
+    h01  = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex)
+    for layer in range(layer_num):
+        if np.mod(layer,2) == 0:
+            h01[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h01_modified
+        if np.mod(layer,2) == 1:
+            h01[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h01_Haldane
+    return h01
+
+def bands_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
+    k_array = np.linspace(0, 2*pi, 300)
+    hamiltonian = functools.partial(hamiltonian_of_bilayer_model, N=N, M=M, t1=t1, t2=t2, phi=phi, tc=tc, layer_num=layer_num)
+    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k_array, hamiltonian, print_show=1)
+    plt, fig, ax = guan.import_plt_and_start_fig_ax(labelsize=25)
+    guan.plot_without_starting_fig_ax(plt, fig, ax, k_array, eigenvalue_array, xlabel='$\mathrm{k_x}$', ylabel='$\mathrm{E}$', style='k', fontsize=30,  y_max=0.2, y_min=-0.2,linewidth=None, markersize=None, color=None, fontfamily='Times New Roman')
+    plt.show()
+
+def conductance_of_bilayer_model_with_disorder_array(N, M, t1, t2, phi, tc, layer_num):
+    h00 = h00_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
+    h01 = h01_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
+    fermi_energy = 0
+    disorder_intensity_array = np.arange(0, 2.5, .05)
+    conductance_array = guan.calculate_conductance_with_disorder_intensity_array(fermi_energy, h00, h01, disorder_intensity_array, length=100, calculation_times=3, print_show=1) # length=2000, calculation_times=20
+    guan.plot(disorder_intensity_array, conductance_array, xlabel='$\mathrm{W_d}$', ylabel='$\mathrm{G(e^2/h)}$', style='o-')
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file
diff --git a/2024.10.21_python_format_output/python_format_output.py b/2024.10.21_python_format_output/python_format_output.py
new file mode 100644
index 0000000..5926515
--- /dev/null
+++ b/2024.10.21_python_format_output/python_format_output.py
@@ -0,0 +1,11 @@
+# 小数的格式化输出（浮点）
+value = 3.141592653589793
+print(f"pi={value:.2f}")
+print("pi={:.2f}".format(value))
+print("pi=%.2f" % value)
+
+# 整数的格式化输出
+value = 13
+print(f"a={value:5d}")
+print("a={:5d}".format(value))
+print("a=%5d" % value)
\ No newline at end of file