Work on ps16-1

ps16-1.py

@@ -2,30 +2,35 @@ import numpy as np
 import pylab as pl
 
 a = 5
-F = 1
+F = 4000
 D = 1
 hbar = 1
 m = 1
-z = 6.933
+z = 2.851
 l = z / a
 V0 = 2
 
+
+def cot(x):
+    return 1 / np.tan(x)
+
+
 x = np.linspace(-10, 10, 1000)
 x_pos = x[x > 0]
 z0 = (a / hbar) * np.sqrt(2 * m * V0)
 print(f'z0 = {z0}')
-k = l * np.tan(l * a)
+k = -l * cot(l * a)
 
 psi_conditions = [
     (0 < x_pos) & (x_pos <= a),
     (x_pos > a)
 ]
 psi_funcs = [
-    lambda x: D * np.cos(l * x),
+    lambda x: D * np.sin(l * x),
     lambda x: F * np.exp(-k * x),
 ]
 psi_half = np.piecewise(x_pos, psi_conditions, psi_funcs)
-psi = np.concatenate((np.flip(psi_half), psi_half))
+psi = np.concatenate((-np.flip(psi_half), psi_half))
 
 V0_graph = np.piecewise(x, [(x < -a), (-a <= x) & (x <= a), (x > a)], [0, -V0, 0])
 

test.py

@@ -0,0 +1,66 @@
+import numpy as np
+import pylab as pl
+from scipy.optimize import fsolve
+
+
+def cot(x):
+    return 1 / np.tan(x)
+
+
+a = 5
+hbar = 1
+m = 1
+z = 2.851
+l = z / a
+V0 = 2
+D = 1
+
+x = np.linspace(-10, 10, 1000)
+x_pos = x[x > 0]
+z0 = (a / hbar) * np.sqrt(2 * m * V0)
+print(f'z0 = {z0}')
+k = -l * cot(l * a)
+
+
+# Define the functions for the wavefunction and its derivative
+def psi_inside(x):
+    return D * np.sin(l * x)
+
+
+def psi_outside(x, F):
+    return F * np.exp(-k * x)
+
+
+def dpsi_inside(x):
+    return D * l * np.cos(l * x)
+
+
+def dpsi_outside(x, F):
+    return -F * k * np.exp(-k * x)
+
+
+# Define the equations for the continuity conditions
+def equations(p):
+    return psi_inside(a) - psi_outside(a, p)
+
+
+# Solve for D and F
+F = fsolve(equations, np.ndarray(1))
+print(f'F = {F}')
+
+psi_conditions = [
+    (0 < x_pos) & (x_pos <= a),
+    (x_pos > a)
+]
+psi_funcs = [
+    lambda x: psi_inside(x),
+    lambda x: psi_outside(x, F),
+]
+psi_half = np.piecewise(x_pos, psi_conditions, psi_funcs)
+psi = np.concatenate((-np.flip(psi_half), psi_half))
+
+V0_graph = np.piecewise(x, [(x < -a), (-a <= x) & (x <= a), (x > a)], [0, -V0, 0])
+
+pl.plot(x, psi)
+pl.plot(x, V0_graph)
+pl.show()