diff --git a/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/Haldane_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/Haldane_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
new file mode 100644
index 0000000..1273557
--- /dev/null
+++ b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/Haldane_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
@@ -0,0 +1,47 @@
+"""
+This code is supported by the website: https://www.guanjihuan.com
+The newest version of this code is on the web page: https://www.guanjihuan.com/archives/22691
+"""
+
+
+import guan
+import numpy as np
+import functools
+
+
+# 2D Haldane lattice
+k1_array = np.linspace(-2*np.pi, 2*np.pi, 100)
+k2_array = np.linspace(-2*np.pi, 2*np.pi, 100)
+eigenvalue_array = guan.calculate_eigenvalue_with_two_parameters(k1_array, k2_array, guan.hamiltonian_of_haldane_model)
+guan.plot_3d_surface(k1_array, k2_array, eigenvalue_array, xlabel='kx', ylabel='ky', zlabel='E')
+
+
+# 2D Haldane lattice for discrete ky array
+Ny = 5
+ky_array = np.linspace(-np.pi, np.pi, Ny*4) # important
+print(ky_array)
+
+kx_array = np.linspace(-np.pi, np.pi, 100)
+i0 = 0
+for ky in ky_array:
+    hamiltonian_function = functools.partial(guan.hamiltonian_of_haldane_model, k2=ky)
+    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(-kx_array, hamiltonian_function)
+    if i0 == 0:
+        eigenvalue_array_for_discrete_ky = eigenvalue_array
+    else:
+        eigenvalue_array_for_discrete_ky = np.append(eigenvalue_array_for_discrete_ky, eigenvalue_array, axis=1)
+    i0 += 1
+
+
+# 1D Haldane ribbon
+hamiltonian_function = functools.partial(guan.hamiltonian_of_haldane_model_in_quasi_one_dimension, N=Ny)
+eigenvalue_array_2 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# 1D Haldane ribbon with periodic boundary condition in y direction
+hamiltonian_function = functools.partial(guan.hamiltonian_of_haldane_model_in_quasi_one_dimension, N=Ny, period=1)
+eigenvalue_array_3 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# Plot figures
+guan.plot_three_array(kx_array, eigenvalue_array_for_discrete_ky, eigenvalue_array_2, eigenvalue_array_3, xlabel='kx', ylabel='E', style_1='-k', style_2='--r', style_3='.y', linewidth_1=3, markersize_2=3, markersize_3=3)
\ No newline at end of file
diff --git a/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/graphene_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/graphene_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
new file mode 100644
index 0000000..1f24a60
--- /dev/null
+++ b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/graphene_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
@@ -0,0 +1,46 @@
+"""
+This code is supported by the website: https://www.guanjihuan.com
+The newest version of this code is on the web page: https://www.guanjihuan.com/archives/22691
+"""
+
+import guan
+import numpy as np
+import functools
+
+
+# 2D graphene lattice
+k1_array = np.linspace(-2*np.pi, 2*np.pi, 300)
+k2_array = np.linspace(-2*np.pi, 2*np.pi, 300)
+eigenvalue_array = guan.calculate_eigenvalue_with_two_parameters(k1_array, k2_array, guan.hamiltonian_of_graphene)
+guan.plot_3d_surface(k1_array, k2_array, eigenvalue_array, xlabel='kx', ylabel='ky', zlabel='E', rcount=300, ccount=300)
+
+
+# 2D graphene lattice for discrete ky array
+Ny = 5
+ky_array = np.linspace(-np.pi, np.pi, Ny*4) # important
+print(ky_array)
+
+kx_array = np.linspace(-np.pi, np.pi, 100)
+i0 = 0
+for ky in ky_array:
+    hamiltonian_function = functools.partial(guan.hamiltonian_of_graphene, k2=ky)
+    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+    if i0 == 0:
+        eigenvalue_array_for_discrete_ky = eigenvalue_array
+    else:
+        eigenvalue_array_for_discrete_ky = np.append(eigenvalue_array_for_discrete_ky, eigenvalue_array, axis=1)
+    i0 += 1
+
+
+# 1D graphene ribbon
+hamiltonian_function = functools.partial(guan.hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension, N=Ny)
+eigenvalue_array_2 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# 1D graphene ribbon with periodic boundary condition in y direction
+hamiltonian_function = functools.partial(guan.hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension, N=Ny, period=1)
+eigenvalue_array_3 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# Plot figures
+guan.plot_three_array(kx_array, eigenvalue_array_for_discrete_ky, eigenvalue_array_2, eigenvalue_array_3, xlabel='kx', ylabel='E', style_1='-k', style_2='--r', style_3='.y', linewidth_1=3, markersize_2=3, markersize_3=3)
\ No newline at end of file
diff --git a/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/square_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/square_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
new file mode 100644
index 0000000..10105b3
--- /dev/null
+++ b/academic_codes/2022.07.14_relation_of_2D_discrete_ky_and_ribbon_with_Ny/square_relation_of_2D_discrete_ky_and_ribbon_with_Ny.py
@@ -0,0 +1,47 @@
+"""
+This code is supported by the website: https://www.guanjihuan.com
+The newest version of this code is on the web page: https://www.guanjihuan.com/archives/22691
+"""
+
+
+import guan
+import numpy as np
+import functools
+
+
+# 2D square lattice
+k1_array = np.linspace(-np.pi, np.pi, 100)
+k2_array = np.linspace(-np.pi, np.pi, 100)
+eigenvalue_array = guan.calculate_eigenvalue_with_two_parameters(k1_array, k2_array, guan.hamiltonian_of_square_lattice)
+guan.plot_3d_surface(k1_array, k2_array, eigenvalue_array, xlabel='kx', ylabel='ky', zlabel='E')
+
+
+# 2D square lattice for discrete ky array
+Ny = 10
+ky_array = np.linspace(-np.pi, np.pi, Ny) # important
+print(ky_array)
+
+kx_array = np.linspace(-np.pi, np.pi, 100)
+i0 = 0
+for ky in ky_array:
+    hamiltonian_function = functools.partial(guan.hamiltonian_of_square_lattice, k2=ky)
+    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+    if i0 == 0:
+        eigenvalue_array_for_discrete_ky = eigenvalue_array
+    else:
+        eigenvalue_array_for_discrete_ky = np.append(eigenvalue_array_for_discrete_ky, eigenvalue_array, axis=1)
+    i0 += 1
+
+
+# 1D square ribbon
+hamiltonian_function = functools.partial(guan.hamiltonian_of_square_lattice_in_quasi_one_dimension, N=Ny)
+eigenvalue_array_2 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# 1D square ribbon with periodic boundary condition in y direction
+hamiltonian_function = functools.partial(guan.hamiltonian_of_square_lattice_in_quasi_one_dimension, N=Ny, period=1)
+eigenvalue_array_3 = guan.calculate_eigenvalue_with_one_parameter(kx_array, hamiltonian_function)
+
+
+# Plot figures
+guan.plot_three_array(kx_array, eigenvalue_array_for_discrete_ky, eigenvalue_array_2, eigenvalue_array_3, xlabel='kx', ylabel='E', style_1='-k', style_2='--r', style_3='.y', linewidth_1=3, markersize_2=1, markersize_3=1)
\ No newline at end of file