Add plots

Reviewed-on: #2

data/probability_deviation.txt

@@ -1,320 +0,0 @@
-0,    1.0000000000e+00,    1.0000000000e+00
-1,    1.0000000000e+00,    1.0000000000e+00
-2,    1.0000000000e+00,    1.0000000000e+00
-3,    1.0000000000e+00,    1.0000000000e+00
-4,    1.0000000000e+00,    1.0000000000e+00
-5,    1.0000000000e+00,    1.0000000000e+00
-6,    1.0000000000e+00,    1.0000000000e+00
-7,    1.0000000000e+00,    1.0000000000e+00
-8,    1.0000000000e+00,    1.0000000000e+00
-9,    1.0000000000e+00,    1.0000000000e+00
-10,    1.0000000000e+00,    1.0000000000e+00
-11,    1.0000000000e+00,    1.0000000000e+00
-12,    1.0000000000e+00,    1.0000000000e+00
-13,    1.0000000000e+00,    1.0000000000e+00
-14,    1.0000000000e+00,    1.0000000000e+00
-15,    1.0000000000e+00,    1.0000000000e+00
-16,    1.0000000000e+00,    1.0000000000e+00
-17,    1.0000000000e+00,    1.0000000000e+00
-18,    1.0000000000e+00,    1.0000000000e+00
-19,    1.0000000000e+00,    1.0000000000e+00
-20,    1.0000000000e+00,    1.0000000000e+00
-21,    1.0000000000e+00,    1.0000000000e+00
-22,    1.0000000000e+00,    1.0000000000e+00
-23,    1.0000000000e+00,    1.0000000000e+00
-24,    1.0000000000e+00,    1.0000000000e+00
-25,    1.0000000000e+00,    1.0000000000e+00
-26,    1.0000000000e+00,    1.0000000000e+00
-27,    1.0000000000e+00,    1.0000000000e+00
-28,    1.0000000000e+00,    1.0000000000e+00
-29,    1.0000000000e+00,    1.0000000000e+00
-30,    1.0000000000e+00,    1.0000000000e+00
-31,    1.0000000000e+00,    1.0000000000e+00
-32,    1.0000000000e+00,    1.0000000000e+00
-33,    1.0000000000e+00,    1.0000000000e+00
-34,    1.0000000000e+00,    1.0000000000e+00
-35,    1.0000000000e+00,    1.0000000000e+00
-36,    1.0000000000e+00,    1.0000000000e+00
-37,    1.0000000000e+00,    1.0000000000e+00
-38,    1.0000000000e+00,    1.0000000000e+00
-39,    1.0000000000e+00,    1.0000000000e+00
-40,    1.0000000000e+00,    1.0000000000e+00
-41,    1.0000000000e+00,    1.0000000000e+00
-42,    1.0000000000e+00,    1.0000000000e+00
-43,    1.0000000000e+00,    1.0000000000e+00
-44,    1.0000000000e+00,    1.0000000000e+00
-45,    1.0000000000e+00,    1.0000000000e+00
-46,    1.0000000000e+00,    1.0000000000e+00
-47,    1.0000000000e+00,    1.0000000000e+00
-48,    1.0000000000e+00,    1.0000000000e+00
-49,    1.0000000000e+00,    1.0000000000e+00
-50,    1.0000000000e+00,    1.0000000000e+00
-51,    1.0000000000e+00,    1.0000000000e+00
-52,    1.0000000000e+00,    1.0000000000e+00
-53,    1.0000000000e+00,    1.0000000000e+00
-54,    1.0000000000e+00,    1.0000000000e+00
-55,    1.0000000000e+00,    1.0000000000e+00
-56,    1.0000000000e+00,    1.0000000000e+00
-57,    1.0000000000e+00,    1.0000000000e+00
-58,    1.0000000000e+00,    1.0000000000e+00
-59,    1.0000000000e+00,    1.0000000000e+00
-60,    1.0000000000e+00,    1.0000000000e+00
-61,    1.0000000000e+00,    1.0000000000e+00
-62,    1.0000000000e+00,    1.0000000000e+00
-63,    1.0000000000e+00,    1.0000000000e+00
-64,    1.0000000000e+00,    1.0000000000e+00
-65,    1.0000000000e+00,    1.0000000000e+00
-66,    1.0000000000e+00,    1.0000000000e+00
-67,    1.0000000000e+00,    1.0000000000e+00
-68,    1.0000000000e+00,    1.0000000000e+00
-69,    1.0000000000e+00,    1.0000000000e+00
-70,    1.0000000000e+00,    1.0000000000e+00
-71,    1.0000000000e+00,    1.0000000000e+00
-72,    1.0000000000e+00,    1.0000000000e+00
-73,    1.0000000000e+00,    1.0000000000e+00
-74,    1.0000000000e+00,    1.0000000000e+00
-75,    1.0000000000e+00,    1.0000000000e+00
-76,    1.0000000000e+00,    1.0000000000e+00
-77,    1.0000000000e+00,    1.0000000000e+00
-78,    1.0000000000e+00,    1.0000000000e+00
-79,    1.0000000000e+00,    1.0000000000e+00
-80,    1.0000000000e+00,    1.0000000000e+00
-81,    1.0000000000e+00,    1.0000000000e+00
-82,    1.0000000000e+00,    1.0000000000e+00
-83,    1.0000000000e+00,    1.0000000000e+00
-84,    1.0000000000e+00,    1.0000000000e+00
-85,    1.0000000000e+00,    1.0000000000e+00
-86,    1.0000000000e+00,    1.0000000000e+00
-87,    1.0000000000e+00,    1.0000000000e+00
-88,    1.0000000000e+00,    1.0000000000e+00
-89,    1.0000000000e+00,    1.0000000000e+00
-90,    1.0000000000e+00,    1.0000000000e+00
-91,    1.0000000000e+00,    1.0000000000e+00
-92,    1.0000000000e+00,    1.0000000000e+00
-93,    1.0000000000e+00,    1.0000000000e+00
-94,    1.0000000000e+00,    1.0000000000e+00
-95,    1.0000000000e+00,    1.0000000000e+00
-96,    1.0000000000e+00,    1.0000000000e+00
-97,    1.0000000000e+00,    1.0000000000e+00
-98,    1.0000000000e+00,    1.0000000000e+00
-99,    1.0000000000e+00,    1.0000000000e+00
-100,    1.0000000000e+00,    1.0000000000e+00
-101,    1.0000000000e+00,    1.0000000000e+00
-102,    1.0000000000e+00,    1.0000000000e+00
-103,    1.0000000000e+00,    1.0000000000e+00
-104,    1.0000000000e+00,    1.0000000000e+00
-105,    1.0000000000e+00,    1.0000000000e+00
-106,    1.0000000000e+00,    1.0000000000e+00
-107,    1.0000000000e+00,    1.0000000000e+00
-108,    1.0000000000e+00,    1.0000000000e+00
-109,    1.0000000000e+00,    1.0000000000e+00
-110,    1.0000000000e+00,    1.0000000000e+00
-111,    1.0000000000e+00,    1.0000000000e+00
-112,    1.0000000000e+00,    1.0000000000e+00
-113,    1.0000000000e+00,    1.0000000000e+00
-114,    1.0000000000e+00,    1.0000000000e+00
-115,    1.0000000000e+00,    1.0000000000e+00
-116,    1.0000000000e+00,    1.0000000000e+00
-117,    1.0000000000e+00,    1.0000000000e+00
-118,    1.0000000000e+00,    1.0000000000e+00
-119,    1.0000000000e+00,    1.0000000000e+00
-120,    1.0000000000e+00,    1.0000000000e+00
-121,    1.0000000000e+00,    1.0000000000e+00
-122,    1.0000000000e+00,    1.0000000000e+00
-123,    1.0000000000e+00,    1.0000000000e+00
-124,    1.0000000000e+00,    1.0000000000e+00
-125,    1.0000000000e+00,    1.0000000000e+00
-126,    1.0000000000e+00,    1.0000000000e+00
-127,    1.0000000000e+00,    1.0000000000e+00
-128,    1.0000000000e+00,    1.0000000000e+00
-129,    1.0000000000e+00,    1.0000000000e+00
-130,    1.0000000000e+00,    1.0000000000e+00
-131,    1.0000000000e+00,    1.0000000000e+00
-132,    1.0000000000e+00,    1.0000000000e+00
-133,    1.0000000000e+00,    1.0000000000e+00
-134,    1.0000000000e+00,    1.0000000000e+00
-135,    1.0000000000e+00,    1.0000000000e+00
-136,    1.0000000000e+00,    1.0000000000e+00
-137,    1.0000000000e+00,    1.0000000000e+00
-138,    1.0000000000e+00,    1.0000000000e+00
-139,    1.0000000000e+00,    1.0000000000e+00
-140,    1.0000000000e+00,    1.0000000000e+00
-141,    1.0000000000e+00,    1.0000000000e+00
-142,    1.0000000000e+00,    1.0000000000e+00
-143,    1.0000000000e+00,    1.0000000000e+00
-144,    1.0000000000e+00,    1.0000000000e+00
-145,    1.0000000000e+00,    1.0000000000e+00
-146,    1.0000000000e+00,    1.0000000000e+00
-147,    1.0000000000e+00,    1.0000000000e+00
-148,    1.0000000000e+00,    1.0000000000e+00
-149,    1.0000000000e+00,    1.0000000000e+00
-150,    1.0000000000e+00,    1.0000000000e+00
-151,    1.0000000000e+00,    1.0000000000e+00
-152,    1.0000000000e+00,    1.0000000000e+00
-153,    1.0000000000e+00,    1.0000000000e+00
-154,    1.0000000000e+00,    1.0000000000e+00
-155,    1.0000000000e+00,    1.0000000000e+00
-156,    1.0000000000e+00,    1.0000000000e+00
-157,    1.0000000000e+00,    1.0000000000e+00
-158,    1.0000000000e+00,    1.0000000000e+00
-159,    1.0000000000e+00,    1.0000000000e+00
-160,    1.0000000000e+00,    1.0000000000e+00
-161,    1.0000000000e+00,    1.0000000000e+00
-162,    1.0000000000e+00,    1.0000000000e+00
-163,    1.0000000000e+00,    1.0000000000e+00
-164,    1.0000000000e+00,    1.0000000000e+00
-165,    1.0000000000e+00,    1.0000000000e+00
-166,    1.0000000000e+00,    1.0000000000e+00
-167,    1.0000000000e+00,    1.0000000000e+00
-168,    1.0000000000e+00,    1.0000000000e+00
-169,    1.0000000000e+00,    1.0000000000e+00
-170,    1.0000000000e+00,    1.0000000000e+00
-171,    1.0000000000e+00,    1.0000000000e+00
-172,    1.0000000000e+00,    1.0000000000e+00
-173,    1.0000000000e+00,    1.0000000000e+00
-174,    1.0000000000e+00,    1.0000000000e+00
-175,    1.0000000000e+00,    1.0000000000e+00
-176,    1.0000000000e+00,    1.0000000000e+00
-177,    1.0000000000e+00,    1.0000000000e+00
-178,    1.0000000000e+00,    1.0000000000e+00
-179,    1.0000000000e+00,    1.0000000000e+00
-180,    1.0000000000e+00,    1.0000000000e+00
-181,    1.0000000000e+00,    1.0000000000e+00
-182,    1.0000000000e+00,    1.0000000000e+00
-183,    1.0000000000e+00,    1.0000000000e+00
-184,    1.0000000000e+00,    1.0000000000e+00
-185,    1.0000000000e+00,    1.0000000000e+00
-186,    1.0000000000e+00,    1.0000000000e+00
-187,    1.0000000000e+00,    1.0000000000e+00
-188,    1.0000000000e+00,    1.0000000000e+00
-189,    1.0000000000e+00,    1.0000000000e+00
-190,    1.0000000000e+00,    1.0000000000e+00
-191,    1.0000000000e+00,    1.0000000000e+00
-192,    1.0000000000e+00,    1.0000000000e+00
-193,    1.0000000000e+00,    1.0000000000e+00
-194,    1.0000000000e+00,    1.0000000000e+00
-195,    1.0000000000e+00,    1.0000000000e+00
-196,    1.0000000000e+00,    1.0000000000e+00
-197,    1.0000000000e+00,    1.0000000000e+00
-198,    1.0000000000e+00,    1.0000000000e+00
-199,    1.0000000000e+00,    1.0000000000e+00
-200,    1.0000000000e+00,    1.0000000000e+00
-201,    1.0000000000e+00,    1.0000000000e+00
-202,    1.0000000000e+00,    1.0000000000e+00
-203,    1.0000000000e+00,    1.0000000000e+00
-204,    1.0000000000e+00,    1.0000000000e+00
-205,    1.0000000000e+00,    1.0000000000e+00
-206,    1.0000000000e+00,    1.0000000000e+00
-207,    1.0000000000e+00,    1.0000000000e+00
-208,    1.0000000000e+00,    1.0000000000e+00
-209,    1.0000000000e+00,    1.0000000000e+00
-210,    1.0000000000e+00,    1.0000000000e+00
-211,    1.0000000000e+00,    1.0000000000e+00
-212,    1.0000000000e+00,    1.0000000000e+00
-213,    1.0000000000e+00,    1.0000000000e+00
-214,    1.0000000000e+00,    1.0000000000e+00
-215,    1.0000000000e+00,    1.0000000000e+00
-216,    1.0000000000e+00,    1.0000000000e+00
-217,    1.0000000000e+00,    1.0000000000e+00
-218,    1.0000000000e+00,    1.0000000000e+00
-219,    1.0000000000e+00,    1.0000000000e+00
-220,    1.0000000000e+00,    1.0000000000e+00
-221,    1.0000000000e+00,    1.0000000000e+00
-222,    1.0000000000e+00,    1.0000000000e+00
-223,    1.0000000000e+00,    1.0000000000e+00
-224,    1.0000000000e+00,    1.0000000000e+00
-225,    1.0000000000e+00,    1.0000000000e+00
-226,    1.0000000000e+00,    1.0000000000e+00
-227,    1.0000000000e+00,    1.0000000000e+00
-228,    1.0000000000e+00,    1.0000000000e+00
-229,    1.0000000000e+00,    1.0000000000e+00
-230,    1.0000000000e+00,    1.0000000000e+00
-231,    1.0000000000e+00,    1.0000000000e+00
-232,    1.0000000000e+00,    1.0000000000e+00
-233,    1.0000000000e+00,    1.0000000000e+00
-234,    1.0000000000e+00,    1.0000000000e+00
-235,    1.0000000000e+00,    1.0000000000e+00
-236,    1.0000000000e+00,    1.0000000000e+00
-237,    1.0000000000e+00,    1.0000000000e+00
-238,    1.0000000000e+00,    1.0000000000e+00
-239,    1.0000000000e+00,    1.0000000000e+00
-240,    1.0000000000e+00,    1.0000000000e+00
-241,    1.0000000000e+00,    1.0000000000e+00
-242,    1.0000000000e+00,    1.0000000000e+00
-243,    1.0000000000e+00,    1.0000000000e+00
-244,    1.0000000000e+00,    1.0000000000e+00
-245,    1.0000000000e+00,    1.0000000000e+00
-246,    1.0000000000e+00,    1.0000000000e+00
-247,    1.0000000000e+00,    1.0000000000e+00
-248,    1.0000000000e+00,    1.0000000000e+00
-249,    1.0000000000e+00,    1.0000000000e+00
-250,    1.0000000000e+00,    1.0000000000e+00
-251,    1.0000000000e+00,    1.0000000000e+00
-252,    1.0000000000e+00,    1.0000000000e+00
-253,    1.0000000000e+00,    1.0000000000e+00
-254,    1.0000000000e+00,    1.0000000000e+00
-255,    1.0000000000e+00,    1.0000000000e+00
-256,    1.0000000000e+00,    1.0000000000e+00
-257,    1.0000000000e+00,    1.0000000000e+00
-258,    1.0000000000e+00,    1.0000000000e+00
-259,    1.0000000000e+00,    1.0000000000e+00
-260,    1.0000000000e+00,    1.0000000000e+00
-261,    1.0000000000e+00,    1.0000000000e+00
-262,    1.0000000000e+00,    1.0000000000e+00
-263,    1.0000000000e+00,    1.0000000000e+00
-264,    1.0000000000e+00,    1.0000000000e+00
-265,    1.0000000000e+00,    1.0000000000e+00
-266,    1.0000000000e+00,    1.0000000000e+00
-267,    1.0000000000e+00,    1.0000000000e+00
-268,    1.0000000000e+00,    1.0000000000e+00
-269,    1.0000000000e+00,    1.0000000000e+00
-270,    1.0000000000e+00,    1.0000000000e+00
-271,    1.0000000000e+00,    1.0000000000e+00
-272,    1.0000000000e+00,    1.0000000000e+00
-273,    1.0000000000e+00,    1.0000000000e+00
-274,    1.0000000000e+00,    1.0000000000e+00
-275,    1.0000000000e+00,    1.0000000000e+00
-276,    1.0000000000e+00,    1.0000000000e+00
-277,    1.0000000000e+00,    1.0000000000e+00
-278,    1.0000000000e+00,    1.0000000000e+00
-279,    1.0000000000e+00,    1.0000000000e+00
-280,    1.0000000000e+00,    1.0000000000e+00
-281,    1.0000000000e+00,    1.0000000000e+00
-282,    1.0000000000e+00,    1.0000000000e+00
-283,    1.0000000000e+00,    1.0000000000e+00
-284,    1.0000000000e+00,    1.0000000000e+00
-285,    1.0000000000e+00,    1.0000000000e+00
-286,    1.0000000000e+00,    1.0000000000e+00
-287,    1.0000000000e+00,    1.0000000000e+00
-288,    1.0000000000e+00,    1.0000000000e+00
-289,    1.0000000000e+00,    1.0000000000e+00
-290,    1.0000000000e+00,    1.0000000000e+00
-291,    1.0000000000e+00,    1.0000000000e+00
-292,    1.0000000000e+00,    1.0000000000e+00
-293,    1.0000000000e+00,    1.0000000000e+00
-294,    1.0000000000e+00,    1.0000000000e+00
-295,    1.0000000000e+00,    1.0000000000e+00
-296,    1.0000000000e+00,    1.0000000000e+00
-297,    1.0000000000e+00,    1.0000000000e+00
-298,    1.0000000000e+00,    1.0000000000e+00
-299,    1.0000000000e+00,    1.0000000000e+00
-300,    1.0000000000e+00,    1.0000000000e+00
-301,    1.0000000000e+00,    1.0000000000e+00
-302,    1.0000000000e+00,    1.0000000000e+00
-303,    1.0000000000e+00,    1.0000000000e+00
-304,    1.0000000000e+00,    1.0000000000e+00
-305,    1.0000000000e+00,    1.0000000000e+00
-306,    1.0000000000e+00,    1.0000000000e+00
-307,    1.0000000000e+00,    1.0000000000e+00
-308,    1.0000000000e+00,    1.0000000000e+00
-309,    1.0000000000e+00,    1.0000000000e+00
-310,    1.0000000000e+00,    1.0000000000e+00
-311,    1.0000000000e+00,    1.0000000000e+00
-312,    1.0000000000e+00,    1.0000000000e+00
-313,    1.0000000000e+00,    1.0000000000e+00
-314,    1.0000000000e+00,    1.0000000000e+00
-315,    1.0000000000e+00,    1.0000000000e+00
-316,    1.0000000000e+00,    1.0000000000e+00
-317,    1.0000000000e+00,    1.0000000000e+00
-318,    1.0000000000e+00,    1.0000000000e+00
-319,    1.0000000000e+00,    1.0000000000e+00

data/probability_deviation_no_slits.txt

@@ -0,0 +1,321 @@
+0	-1.33227e-14
+2.5e-05	-1.37668e-14
+5e-05	-1.24345e-14
+7.5e-05	-1.15463e-14
+0.0001	-1.17684e-14
+0.000125	-1.06581e-14
+0.00015	-1.02141e-14
+0.000175	-1.02141e-14
+0.0002	-1.02141e-14
+0.000225	-9.76996e-15
+0.00025	-9.76996e-15
+0.000275	-1.15463e-14
+0.0003	-7.99361e-15
+0.000325	-1.08802e-14
+0.00035	-9.10383e-15
+0.000375	-9.32587e-15
+0.0004	-1.06581e-14
+0.000425	-9.32587e-15
+0.00045	-1.26565e-14
+0.000475	-1.19904e-14
+0.0005	-1.15463e-14
+0.000525	-9.54792e-15
+0.00055	-1.37668e-14
+0.000575	-1.26565e-14
+0.0006	-9.76996e-15
+0.000625	-7.99361e-15
+0.00065	-9.76996e-15
+0.000675	-1.02141e-14
+0.0007	-1.24345e-14
+0.000725	-1.02141e-14
+0.00075	-8.88178e-15
+0.000775	-1.15463e-14
+0.0008	-7.99361e-15
+0.000825	-8.21565e-15
+0.00085	-8.43769e-15
+0.000875	-7.10543e-15
+0.0009	-1.15463e-14
+0.000925	-1.19904e-14
+0.00095	-7.77156e-15
+0.000975	-6.43929e-15
+0.001	-7.99361e-15
+0.001025	-8.21565e-15
+0.00105	-9.76996e-15
+0.001075	-8.88178e-15
+0.0011	-1.06581e-14
+0.001125	-9.32587e-15
+0.00115	-1.19904e-14
+0.001175	-9.54792e-15
+0.0012	-7.54952e-15
+0.001225	-9.54792e-15
+0.00125	-1.06581e-14
+0.001275	-8.43769e-15
+0.0013	-8.88178e-15
+0.001325	-7.99361e-15
+0.00135	-7.10543e-15
+0.001375	-1.08802e-14
+0.0014	-8.88178e-15
+0.001425	-7.10543e-15
+0.00145	-9.32587e-15
+0.001475	-7.54952e-15
+0.0015	-1.11022e-14
+0.001525	-1.19904e-14
+0.00155	-1.19904e-14
+0.001575	-9.76996e-15
+0.0016	-1.26565e-14
+0.001625	-9.76996e-15
+0.00165	-1.19904e-14
+0.001675	-9.32587e-15
+0.0017	-1.46549e-14
+0.001725	-1.46549e-14
+0.00175	-1.11022e-14
+0.001775	-1.04361e-14
+0.0018	-1.5099e-14
+0.001825	-9.99201e-15
+0.00185	-1.19904e-14
+0.001875	-9.99201e-15
+0.0019	-1.46549e-14
+0.001925	-1.19904e-14
+0.00195	-1.33227e-14
+0.001975	-1.37668e-14
+0.002	-1.28786e-14
+0.002025	-1.28786e-14
+0.00205	-1.28786e-14
+0.002075	-8.43769e-15
+0.0021	-9.32587e-15
+0.002125	-1.24345e-14
+0.00215	-1.33227e-14
+0.002175	-1.31006e-14
+0.0022	-1.37668e-14
+0.002225	-1.59872e-14
+0.00225	-1.37668e-14
+0.002275	-1.4877e-14
+0.0023	-1.55431e-14
+0.002325	-1.37668e-14
+0.00235	-1.55431e-14
+0.002375	-1.35447e-14
+0.0024	-1.46549e-14
+0.002425	-1.55431e-14
+0.00245	-1.55431e-14
+0.002475	-1.42109e-14
+0.0025	-1.64313e-14
+0.002525	-1.46549e-14
+0.00255	-1.37668e-14
+0.002575	-1.39888e-14
+0.0026	-1.59872e-14
+0.002625	-1.39888e-14
+0.00265	-1.55431e-14
+0.002675	-1.44329e-14
+0.0027	-1.5099e-14
+0.002725	-1.42109e-14
+0.00275	-1.39888e-14
+0.002775	-1.37668e-14
+0.0028	-1.37668e-14
+0.002825	-1.5099e-14
+0.00285	-1.55431e-14
+0.002875	-1.35447e-14
+0.0029	-1.55431e-14
+0.002925	-1.5099e-14
+0.00295	-1.46549e-14
+0.002975	-1.42109e-14
+0.003	-1.55431e-14
+0.003025	-1.42109e-14
+0.00305	-1.62093e-14
+0.003075	-1.70974e-14
+0.0031	-1.33227e-14
+0.003125	-1.64313e-14
+0.00315	-1.68754e-14
+0.003175	-1.68754e-14
+0.0032	-1.19904e-14
+0.003225	-1.37668e-14
+0.00325	-1.5099e-14
+0.003275	-1.37668e-14
+0.0033	-1.5099e-14
+0.003325	-1.55431e-14
+0.00335	-1.19904e-14
+0.003375	-1.5099e-14
+0.0034	-1.79856e-14
+0.003425	-1.39888e-14
+0.00345	-1.37668e-14
+0.003475	-1.46549e-14
+0.0035	-1.46549e-14
+0.003525	-1.5099e-14
+0.00355	-1.46549e-14
+0.003575	-1.55431e-14
+0.0036	-1.35447e-14
+0.003625	-1.53211e-14
+0.00365	-1.33227e-14
+0.003675	-1.37668e-14
+0.0037	-1.33227e-14
+0.003725	-1.37668e-14
+0.00375	-1.46549e-14
+0.003775	-1.33227e-14
+0.0038	-1.73195e-14
+0.003825	-1.73195e-14
+0.00385	-1.35447e-14
+0.003875	-1.68754e-14
+0.0039	-1.64313e-14
+0.003925	-1.42109e-14
+0.00395	-1.82077e-14
+0.003975	-1.26565e-14
+0.004	-1.86517e-14
+0.004025	-1.46549e-14
+0.00405	-1.5099e-14
+0.004075	-1.68754e-14
+0.0041	-1.46549e-14
+0.004125	-1.64313e-14
+0.00415	-1.19904e-14
+0.004175	-1.17684e-14
+0.0042	-1.73195e-14
+0.004225	-1.28786e-14
+0.00425	-1.28786e-14
+0.004275	-1.4877e-14
+0.0043	-1.24345e-14
+0.004325	-1.19904e-14
+0.00435	-7.54952e-15
+0.004375	-9.10383e-15
+0.0044	-1.46549e-14
+0.004425	-8.21565e-15
+0.00445	-7.10543e-15
+0.004475	-4.21885e-15
+0.0045	-1.15463e-14
+0.004525	-6.21725e-15
+0.00455	-3.9968e-15
+0.004575	-7.54952e-15
+0.0046	-1.77636e-15
+0.004625	-3.55271e-15
+0.00465	-4.44089e-15
+0.004675	-4.88498e-15
+0.0047	-1.33227e-15
+0.004725	-6.21725e-15
+0.00475	-5.32907e-15
+0.004775	-4.88498e-15
+0.0048	0
+0.004825	2.66454e-15
+0.00485	-2.66454e-15
+0.004875	3.10862e-15
+0.0049	-3.77476e-15
+0.004925	1.9984e-15
+0.00495	-2.66454e-15
+0.004975	9.99201e-16
+0.005	1.0103e-14
+0.005025	-4.44089e-16
+0.00505	0
+0.005075	3.55271e-15
+0.0051	1.88738e-15
+0.005125	1.55431e-15
+0.00515	8.88178e-16
+0.005175	3.9968e-15
+0.0052	6.21725e-15
+0.005225	0
+0.00525	7.10543e-15
+0.005275	1.77636e-15
+0.0053	1.36557e-14
+0.005325	1.42109e-14
+0.00535	3.77476e-15
+0.005375	7.32747e-15
+0.0054	1.27676e-14
+0.005425	1.4877e-14
+0.00545	7.77156e-15
+0.005475	9.54792e-15
+0.0055	1.5099e-14
+0.005525	9.10383e-15
+0.00555	1.79856e-14
+0.005575	1.84297e-14
+0.0056	1.62093e-14
+0.005625	9.10383e-15
+0.00565	1.59872e-14
+0.005675	1.33227e-14
+0.0057	9.54792e-15
+0.005725	1.42109e-14
+0.00575	9.32587e-15
+0.005775	1.53211e-14
+0.0058	1.82077e-14
+0.005825	1.59872e-14
+0.00585	2.10942e-14
+0.005875	1.22125e-14
+0.0059	2.07612e-14
+0.005925	1.95399e-14
+0.00595	1.69864e-14
+0.005975	2.02061e-14
+0.006	2.35367e-14
+0.006025	1.9984e-14
+0.00605	1.60982e-14
+0.006075	1.77636e-14
+0.0061	1.67644e-14
+0.006125	2.53131e-14
+0.00615	2.5091e-14
+0.006175	2.35367e-14
+0.0062	2.68674e-14
+0.006225	3.17524e-14
+0.00625	2.75335e-14
+0.006275	1.9984e-14
+0.0063	2.4869e-14
+0.006325	2.57572e-14
+0.00635	2.62013e-14
+0.006375	2.40918e-14
+0.0064	2.44249e-14
+0.006425	2.62013e-14
+0.00645	2.75335e-14
+0.006475	3.44169e-14
+0.0065	2.08722e-14
+0.006525	2.62013e-14
+0.00655	3.15303e-14
+0.006575	3.04201e-14
+0.0066	2.55351e-14
+0.006625	3.53051e-14
+0.00665	3.19744e-14
+0.006675	2.78666e-14
+0.0067	3.61933e-14
+0.006725	3.03091e-14
+0.00675	3.21965e-14
+0.006775	3.17524e-14
+0.0068	3.88578e-14
+0.006825	3.64153e-14
+0.00685	2.77556e-14
+0.006875	2.52021e-14
+0.0069	2.83107e-14
+0.006925	2.17604e-14
+0.00695	2.22045e-14
+0.006975	2.38698e-14
+0.007	1.9762e-14
+0.007025	1.84297e-14
+0.00705	1.74305e-14
+0.007075	1.14353e-14
+0.0071	1.60982e-14
+0.007125	9.99201e-15
+0.00715	1.39888e-14
+0.007175	1.31006e-14
+0.0072	1.36557e-14
+0.007225	1.23235e-14
+0.00725	1.11022e-14
+0.007275	1.07692e-14
+0.0073	1.0103e-14
+0.007325	1.34337e-14
+0.00735	9.65894e-15
+0.007375	8.88178e-15
+0.0074	5.77316e-15
+0.007425	5.55112e-15
+0.00745	5.9952e-15
+0.007475	6.88338e-15
+0.0075	2.22045e-15
+0.007525	2.22045e-15
+0.00755	-6.66134e-16
+0.007575	1.77636e-15
+0.0076	4.44089e-15
+0.007625	1.33227e-15
+0.00765	4.32987e-15
+0.007675	4.88498e-15
+0.0077	1.22125e-15
+0.007725	4.44089e-16
+0.00775	-2.22045e-15
+0.007775	-8.88178e-16
+0.0078	-3.55271e-15
+0.007825	0
+0.00785	-2.22045e-15
+0.007875	-3.77476e-15
+0.0079	-6.88338e-15
+0.007925	-3.55271e-15
+0.00795	-4.44089e-15
+0.007975	-2.88658e-15
+0.008	-7.54952e-15

data/probability_deviation_slits.txt

@@ -0,0 +1,321 @@
+0	-4.88498e-15
+2.5e-05	-8.88178e-15
+5e-05	-7.54952e-15
+7.5e-05	-9.10383e-15
+0.0001	-5.9952e-15
+0.000125	-9.76996e-15
+0.00015	-7.10543e-15
+0.000175	-6.66134e-15
+0.0002	-8.65974e-15
+0.000225	-5.32907e-15
+0.00025	-6.88338e-15
+0.000275	-6.21725e-15
+0.0003	-7.10543e-15
+0.000325	-5.77316e-15
+0.00035	-6.21725e-15
+0.000375	-6.88338e-15
+0.0004	-7.99361e-15
+0.000425	-8.88178e-15
+0.00045	-6.21725e-15
+0.000475	-8.21565e-15
+0.0005	-1.19904e-14
+0.000525	-7.54952e-15
+0.00055	-6.88338e-15
+0.000575	-4.88498e-15
+0.0006	-1.02141e-14
+0.000625	-6.21725e-15
+0.00065	-6.21725e-15
+0.000675	-7.32747e-15
+0.0007	-4.66294e-15
+0.000725	-5.10703e-15
+0.00075	-8.21565e-15
+0.000775	-5.32907e-15
+0.0008	-8.21565e-15
+0.000825	-7.10543e-15
+0.00085	-7.54952e-15
+0.000875	-6.21725e-15
+0.0009	-3.55271e-15
+0.000925	-9.99201e-15
+0.00095	-8.21565e-15
+0.000975	-5.9952e-15
+0.001	-5.10703e-15
+0.001025	-6.88338e-15
+0.00105	-6.21725e-15
+0.001075	-3.55271e-15
+0.0011	-7.32747e-15
+0.001125	-8.21565e-15
+0.00115	-4.88498e-15
+0.001175	-2.44249e-15
+0.0012	2.22045e-16
+0.001225	-6.66134e-15
+0.00125	-3.10862e-15
+0.001275	-5.32907e-15
+0.0013	-5.9952e-15
+0.001325	-5.77316e-15
+0.00135	1.66533e-15
+0.001375	-1.55431e-15
+0.0014	1.11022e-15
+0.001425	-8.43769e-15
+0.00145	-8.88178e-16
+0.001475	-5.77316e-15
+0.0015	-1.55431e-15
+0.001525	-2.66454e-15
+0.00155	-7.77156e-15
+0.001575	-7.99361e-15
+0.0016	-5.77316e-15
+0.001625	-5.10703e-15
+0.00165	-1.17684e-14
+0.001675	-9.10383e-15
+0.0017	-1.15463e-14
+0.001725	-1.22125e-14
+0.00175	-1.02141e-14
+0.001775	-7.77156e-15
+0.0018	-1.06581e-14
+0.001825	-1.02141e-14
+0.00185	-7.10543e-15
+0.001875	-7.10543e-15
+0.0019	-9.54792e-15
+0.001925	-1.08802e-14
+0.00195	-1.9762e-14
+0.001975	-1.19904e-14
+0.002	-1.42109e-14
+0.002025	-9.76996e-15
+0.00205	-1.35447e-14
+0.002075	-1.17684e-14
+0.0021	-8.43769e-15
+0.002125	-1.24345e-14
+0.00215	-1.70974e-14
+0.002175	-1.17684e-14
+0.0022	-7.10543e-15
+0.002225	-1.26565e-14
+0.00225	-1.73195e-14
+0.002275	-1.82077e-14
+0.0023	-1.13243e-14
+0.002325	-1.88738e-14
+0.00235	-1.37668e-14
+0.002375	-1.28786e-14
+0.0024	-1.90958e-14
+0.002425	-1.5099e-14
+0.00245	-1.37668e-14
+0.002475	-2.04281e-14
+0.0025	-1.31006e-14
+0.002525	-1.24345e-14
+0.00255	-1.39888e-14
+0.002575	-1.26565e-14
+0.0026	-1.11022e-14
+0.002625	-1.35447e-14
+0.00265	-1.77636e-14
+0.002675	-1.88738e-14
+0.0027	-1.17684e-14
+0.002725	-1.22125e-14
+0.00275	-1.22125e-14
+0.002775	-1.55431e-14
+0.0028	-1.33227e-14
+0.002825	-9.99201e-15
+0.00285	-1.59872e-14
+0.002875	-1.59872e-14
+0.0029	-1.57652e-14
+0.002925	-1.19904e-14
+0.00295	-1.39888e-14
+0.002975	-9.99201e-15
+0.003	-1.31006e-14
+0.003025	-1.62093e-14
+0.00305	-1.37668e-14
+0.003075	-1.79856e-14
+0.0031	-1.35447e-14
+0.003125	-1.59872e-14
+0.00315	-1.35447e-14
+0.003175	-1.35447e-14
+0.0032	-1.77636e-14
+0.003225	-1.31006e-14
+0.00325	-1.02141e-14
+0.003275	-1.11022e-14
+0.0033	-1.37668e-14
+0.003325	-1.44329e-14
+0.00335	-1.06581e-14
+0.003375	-1.39888e-14
+0.0034	-1.55431e-14
+0.003425	-9.32587e-15
+0.00345	-1.31006e-14
+0.003475	-1.02141e-14
+0.0035	-1.17684e-14
+0.003525	-1.68754e-14
+0.00355	-1.64313e-14
+0.003575	-1.39888e-14
+0.0036	-1.62093e-14
+0.003625	-1.15463e-14
+0.00365	-1.28786e-14
+0.003675	-2.06501e-14
+0.0037	-1.62093e-14
+0.003725	-1.33227e-14
+0.00375	-1.55431e-14
+0.003775	-1.64313e-14
+0.0038	-1.4877e-14
+0.003825	-1.02141e-14
+0.00385	-1.28786e-14
+0.003875	-1.68754e-14
+0.0039	-1.31006e-14
+0.003925	-1.15463e-14
+0.00395	-1.68754e-14
+0.003975	-1.44329e-14
+0.004	-1.44329e-14
+0.004025	-1.84297e-14
+0.00405	-1.17684e-14
+0.004075	-1.44329e-14
+0.0041	-1.59872e-14
+0.004125	-1.39888e-14
+0.00415	-1.28786e-14
+0.004175	-1.35447e-14
+0.0042	-1.04361e-14
+0.004225	-1.15463e-14
+0.00425	-1.31006e-14
+0.004275	-8.88178e-15
+0.0043	-1.46549e-14
+0.004325	-1.9762e-14
+0.00435	-1.4877e-14
+0.004375	-1.06581e-14
+0.0044	-1.19904e-14
+0.004425	-1.19904e-14
+0.00445	-1.28786e-14
+0.004475	-1.93179e-14
+0.0045	-1.02141e-14
+0.004525	-1.37668e-14
+0.00455	-1.55431e-14
+0.004575	-1.22125e-14
+0.0046	-1.68754e-14
+0.004625	-1.66533e-14
+0.00465	-1.19904e-14
+0.004675	-1.5099e-14
+0.0047	-1.66533e-14
+0.004725	-1.86517e-14
+0.00475	-1.44329e-14
+0.004775	-1.04361e-14
+0.0048	-1.28786e-14
+0.004825	-1.53211e-14
+0.00485	-1.79856e-14
+0.004875	-1.37668e-14
+0.0049	-1.57652e-14
+0.004925	-1.82077e-14
+0.00495	-1.86517e-14
+0.004975	-1.13243e-14
+0.005	-1.42109e-14
+0.005025	-9.76996e-15
+0.00505	-9.76996e-15
+0.005075	-9.76996e-15
+0.0051	-1.11022e-14
+0.005125	-1.11022e-14
+0.00515	-1.62093e-14
+0.005175	-1.02141e-14
+0.0052	-1.55431e-14
+0.005225	-1.37668e-14
+0.00525	-1.46549e-14
+0.005275	-1.24345e-14
+0.0053	-1.77636e-14
+0.005325	-7.99361e-15
+0.00535	-1.02141e-14
+0.005375	-1.53211e-14
+0.0054	-1.11022e-14
+0.005425	-1.64313e-14
+0.00545	-1.66533e-14
+0.005475	-1.55431e-14
+0.0055	-1.22125e-14
+0.005525	-1.59872e-14
+0.00555	-1.42109e-14
+0.005575	-1.37668e-14
+0.0056	-1.33227e-14
+0.005625	-1.24345e-14
+0.00565	-1.5099e-14
+0.005675	-9.54792e-15
+0.0057	-1.68754e-14
+0.005725	-1.42109e-14
+0.00575	-1.68754e-14
+0.005775	-1.73195e-14
+0.0058	-1.35447e-14
+0.005825	-1.26565e-14
+0.00585	-1.11022e-14
+0.005875	-1.08802e-14
+0.0059	-1.86517e-14
+0.005925	-1.11022e-14
+0.00595	-1.39888e-14
+0.005975	-1.75415e-14
+0.006	-1.5099e-14
+0.006025	-1.33227e-14
+0.00605	-1.19904e-14
+0.006075	-1.37668e-14
+0.0061	-1.57652e-14
+0.006125	-1.77636e-14
+0.00615	-1.11022e-14
+0.006175	-1.75415e-14
+0.0062	-1.68754e-14
+0.006225	-1.79856e-14
+0.00625	-1.44329e-14
+0.006275	-1.19904e-14
+0.0063	-1.06581e-14
+0.006325	-1.22125e-14
+0.00635	-1.70974e-14
+0.006375	-1.70974e-14
+0.0064	-1.55431e-14
+0.006425	-1.68754e-14
+0.00645	-1.5099e-14
+0.006475	-1.39888e-14
+0.0065	-1.64313e-14
+0.006525	-1.46549e-14
+0.00655	-1.93179e-14
+0.006575	-1.24345e-14
+0.0066	-1.31006e-14
+0.006625	-1.26565e-14
+0.00665	-1.37668e-14
+0.006675	-1.39888e-14
+0.0067	-1.37668e-14
+0.006725	-1.59872e-14
+0.00675	-1.64313e-14
+0.006775	-1.77636e-14
+0.0068	-1.42109e-14
+0.006825	-1.59872e-14
+0.00685	-1.37668e-14
+0.006875	-1.57652e-14
+0.0069	-1.59872e-14
+0.006925	-1.22125e-14
+0.00695	-1.19904e-14
+0.006975	-1.02141e-14
+0.007	-1.33227e-14
+0.007025	-1.28786e-14
+0.00705	-1.82077e-14
+0.007075	-1.79856e-14
+0.0071	-1.42109e-14
+0.007125	-1.73195e-14
+0.00715	-1.64313e-14
+0.007175	-1.13243e-14
+0.0072	-1.5099e-14
+0.007225	-1.42109e-14
+0.00725	-1.24345e-14
+0.007275	-1.08802e-14
+0.0073	-1.70974e-14
+0.007325	-1.44329e-14
+0.00735	-1.55431e-14
+0.007375	-1.5099e-14
+0.0074	-1.42109e-14
+0.007425	-1.24345e-14
+0.00745	-1.19904e-14
+0.007475	-1.59872e-14
+0.0075	-1.37668e-14
+0.007525	-1.39888e-14
+0.00755	-1.75415e-14
+0.007575	-1.62093e-14
+0.0076	-1.90958e-14
+0.007625	-1.33227e-14
+0.00765	-1.19904e-14
+0.007675	-1.44329e-14
+0.0077	-1.64313e-14
+0.007725	-2.26485e-14
+0.00775	-1.55431e-14
+0.007775	-1.68754e-14
+0.0078	-1.59872e-14
+0.007825	-1.42109e-14
+0.00785	-1.59872e-14
+0.007875	-1.46549e-14
+0.0079	-1.5099e-14
+0.007925	-1.15463e-14
+0.00795	-1.64313e-14
+0.007975	-1.35447e-14
+0.008	-1.37668e-14

include/WaveSimulation.hpp

@@ -21,31 +21,96 @@
 class WaveSimulation {
 protected:
     uint32_t M;
+    int32_t N;
+    arma::cx_mat V;
+    arma::cx_mat U;
     arma::sp_cx_mat B;
+    arma::sp_cx_mat A;
     double h;
     double dt;
     double T;
 
-    void build_A();
-    void build_B();
+    /*  @brief Initialize the U matrix using an unormalized Gaussian wave
+     * packet.
+     *
+     *  @param x_c The center of the packet in the x direction.
+     *  @param y_c The center of the packet in the y direction.
+     *  @param sigma_x The The initial width in the x direction.
+     *  @param sigma_y The The initial width in the y direction.
+     *  @param p_x The wave packet momentum in the x direction.
+     *  @param p_y The wave packet momentum in the y direction.
+     * **/
     void initialize_U(double x_c, double y_c, double sigma_x, double sigma_y,
                       double p_x, double p_y);
-    void build_V(double thickness, double pos_x, double aperture_separation,
-                 double aperture, uint32_t slits);
+
+    /*  @brief Initialize the V matrix.
+     *
+     *  @param thickness The thickness of the wall in the x direction.
+     *  @param pos_x The center of the wall in the x direction.
+     *  @param ap_sep The separation between each aperture.
+     *  @param ap The aperture width.
+     *  @param slits The number of slits.
+     * **/
+    void initialize_V(double thickness, double pos_x,
+                      double aperture_separation, double aperture,
+                      uint32_t slits);
+
+    /*  @brief Initialize the V matrix with no wall.
+     * **/
+    void initialize_V();
+
+    /*  @brief Initialize the A matrix according to the Crank-Nicolson method
+     * **/
+    void initialize_A();
+
+    /*  @brief Initialize the B matrix according to the Crank-Nicolson method
+     * **/
+    void initialize_B();
+
 public:
-    int32_t N;
-    arma::cx_mat V;
-    arma::cx_mat U;
-    arma::sp_cx_mat A;
+    /*  @brief Constructor for the WaveSimulation class.
+     *
+     *  @param h The step size in the x and y direction.
+     *  @param dt The step size in the temporal dimension.
+     *  @param T The total time to simulate.
+     *  @param x_c The center of the packet in the x direction.
+     *  @param y_c The center of the packet in the y direction.
+     *  @param sigma_x The The initial width in the x direction.
+     *  @param sigma_y The The initial width in the y direction.
+     *  @param p_x The wave packet momentum in the x direction.
+     *  @param p_y The wave packet momentum in the y direction.
+     *  @param thickness The thickness of the wall in the x direction.
+     *  @param pos_x The center of the wall in the x direction.
+     *  @param ap_sep The separation between each aperture.
+     *  @param ap The aperture width.
+     *  @param slits The number of slits.
+     * **/
     WaveSimulation(double h, double dt, double T, double x_c, double y_c,
                    double sigma_x, double sigma_y, double p_x, double p_y,
                    double thickness, double pos_x, double ap_sep, double ap,
                    uint32_t slits);
+
+    /*  @brief Constructor for the WaveSimulation class with no wall.
+     *
+     *  @param h The step size in the x and y direction.
+     *  @param dt The step size in the temporal dimension.
+     *  @param T The total time to simulate.
+     *  @param x_c The center of the packet in the x direction.
+     *  @param y_c The center of the packet in the y direction.
+     *  @param sigma_x The The initial width in the x direction.
+     *  @param sigma_y The The initial width in the y direction.
+     *  @param p_x The wave packet momentum in the x direction.
+     *  @param p_y The wave packet momentum in the y direction.
+     * **/
     WaveSimulation(double h, double dt, double T, double x_c, double y_c,
                    double sigma_x, double sigma_y, double p_x, double p_y);
-    virtual void solve(std::ofstream &ofile);
-    void write_U(std::ofstream &ofile);
+
     void step();
+    void solve(std::string outfile, bool write_each_step = false);
+    void solve(std::string outfile, std::vector<double> &steps);
+    void probability_deviation(std::string outfile,
+                               bool write_each_step = false);
+    void write_U(std::ofstream &ofile);
 };
 
 #endif

latex/Makefile

@@ -1,4 +1,4 @@
-SRC=ising_model
+SRC=schrodinger_simulation
 PDFLATEX=pdflatex
 BIBTEX=bibtex
 

latex/draft.tex

@@ -0,0 +1,440 @@
+\documentclass[./schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\section{Theoretical background}\label{sec:theory}
+Scientists have made use of light in devices leading to...
+
+The use of light have facilitated many discoveries through history. The invention  
+of the compound microscope, which used lenses to focus light, led to the first 
+observation of animal cells. And the telescope, which led to a greater understanding 
+of the universe. 
+
+The nature of light has long been a subject of interest and discussion. Around year
+1600 the first microscope was invented, which made use of light through lenses.  
+The observation made using the microscope led to a greater understanding of the 
+world at a microscopic scale. However, the study of light itself have led to (something).
+
+Through the 1600s, the view shifted from particle to wave. 
+% Introduction?
+The nature of light has long been a subject of interest and discussion, from the 1500s 
+and the invention of microscopes, through the 1700s where both a particle theory 
+and a wave theory. In the late 1600s, Christiaan Huygens proposed the wave theory of light, 
+which was challenged by Isaac Newton's particle theory. The particle theory was 
+the leading theory in the beginning of 1800s, when Thomas Young demonstrated 
+the interference of light, through his double-slit experiment. The wave theory found 
+new hold. 
+
+In the 1800s, the study of ideal black bodies done by Gustav R. Kirchhoff, lead to a 
+better understanding of heat radiation. Wilhelm Wien started working on determining the 
+spectral energy distribution, and Wien's law. The law did make sense for high frequencies, 
+however, as there were inconsistencies when the frequency were lower than a threshold value. 
+Wiens law led to an exponential curve, which disagree with the law of conservation.
+
+Max Planck guessed a solution to the problem, which led to Planck's radiation law.
+He later derived the radiation law using Boltzmanns statistical interpretation of 
+the second law of thermodynamics. 
+Plancks findings gave rise to the quantum hypothesis, and later Einsteins wave-particle 
+duality \cite{britannica:1998:planck}.  
+
+For small atoms classical mechanics are not sufficient in describing the position of a particle,
+and Heisenberg uncertainty priciple suggest that the particles have a wavelike behavior.
+The wave-particle duality was later proposed to apply to particles by Louis de Broglie,
+which inspired Erwin Schrödinger who proposed a wave function to describe the quantum 
+state of a particle, resulting in the wave equation.
+
+% Methods?
+% Double-slit experiment
+Thomas Young first performed the double-slit experiment in 1801, to demonstrate 
+the principle of interference of light \cite{britannica:2023:young}, postulating 
+light as waves. In the study of blackbodies, scientists were not able to describe 
+the radiated intensity of increased frequencies using classical mechanichs, as they 
+contradicted the principle of conservation of energy \cite{britannica:1998:planck}. 
+Max Planck assumed that the radiated energy consist of discrete values, or quanta, 
+to describe the peak in radiated energy. 
+
+% Reference: https://tex.stackexchange.com/questions/469109/how-to-change-arrowheads-to-lie-on-a-plane
+\begin{tikzpicture}[scale=1.25,every node/.append style={transform shape}]
+    \foreach \x in {-1,-0.75,...,0} {
+        \draw (\x,-1) -- (\x,1);
+    }
+    \draw[fill=black!10] (0.5,-2,-1) -- (0.5,-2,1) -- (0.5,2,1) -- (0.5,2,-1) -- (0.5,-2,-1);
+    \fill (0.5,0,0) circle (0.05);
+    \foreach \r in {0.25,0.5,...,1.75} {
+        \draw (0.5,0) ++(-60:\r) arc (-60:60:\r);
+    }
+    \draw[fill=black!10] (2,-2,-1) -- (2,-2,1) -- (2,2,1) -- (2,2,-1) -- (2,-2,-1);
+    %\fill (2,0.5) circle (0.05) (2,-0.5) circle (0.05);
+    \foreach \r in {0.25,0.5,...,2} {
+        \draw (2,0.5) ++(-60:\r) arc (-60:60:\r);
+        \draw (2,-0.5) ++(-60:\r) arc (-60:60:\r);
+    }
+    \draw[fill=black!10] (4,-2,-1) -- (4,-2,1) -- (4,2,1) -- (4,2,-1) -- (4,-2,-1);
+%       LABELLING
+        \begin{scope}[canvas is yz plane at x=2,rotate=-90]
+            \node[circle,inner sep=0.5mm,fill,label=above:{S${}_1$}] at (0,0.5){};
+            \node[circle,inner sep=0.5mm,fill,label=below:{S${}_2$}] at (0,-0.5) {};
+            \draw[dash pattern=on 1.5pt off 1pt,thin] (-0.5,0.5) -- (0,0.5)
+            (-0.5,-0.5) -- (0,-0.5);
+            \pgflowlevelsynccm          
+            \draw[|<->|] (-0.5,0.5) -- (-0.5,-0.5) node[midway,right=-0.1cm] {d};
+        \end{scope}
+        \begin{scope}[canvas is yz plane at x=0.5,rotate=-90]
+            \node[below left=-0.1cm] at (0,0) {S${}_0$};
+        \end{scope}
+    \begin{scope}[xshift=4cm,yshift=2cm,rotate=-90,canvas is xy plane at z=0]
+        \fill[white] (0,0) rectangle (4,4);
+        \begin{axis}[
+            width=5.575cm,
+            xmin=-0.5,
+            xmax=0.5,
+            ticks=none
+        ]
+            \addplot [samples=1000,blue
+            ]
+            {(cos(deg(5*pi*sin(deg(x)))))^(2)*((sin(deg(4*pi*sin(deg(x)))))/(4*pi*sin(deg(x))))^(2)};
+        \end{axis}
+    \end{scope}
+    \draw[thin,densely dashed,blue] (2,0) -- (6.9,0);
+    \draw[thin,densely dashed,blue] (2,0) -- +(15:2.5);
+    \draw[thin,densely dashed,blue] (2,0) -- +(-15:2.5);
+    \draw[thin,densely dashed,blue] (2,0) -- +(32:2.5);
+    \draw[thin,densely dashed,blue] (2,0) -- +(-32:2.5);
+\end{tikzpicture}
+
+% Author: Izaak Neutelings (June 2020)
+% Inspiration:
+%  https://courses.physics.ucsd.edu/2011/Summer/session1/physics2c/diffraction.pdf
+%  https://tex.stackexchange.com/questions/201830/periodic-shading-in-tikz
+% \begin{tikzpicture}[
+%     nodal/.style={mylightgreen,dashed,very thin},
+%     declare function={
+%       %xnode(\n,\dn,\lam,\f) = sqrt( (\n^2+(\n+\dn)^2)*\lambd^2/2 - (\n^2-(\n+\dn)^2)^2*\lambd^4/(4*\a^2) - \a^2/4 );
+%       xnode(\n,\dn,\lam,\f) = \lam/\f*sqrt( \n^2*(\f^2-\dn^2)+\n*\dn*(\f^2-\dn^2)+\dn^2*\f^2/2-(\f^4+\dn^4)/4 );
+%       ynode(\n,\dn,\lam,\a) = (2*\n*\dn+\dn^2)*\lam/(2*\f);
+%       intensity(\y,\lam,\a,\L) = cos(180*\a*\y/(2*\lam*sqrt(\L*\L+\y*\y)))^2;
+%     }
+%   ]
+  
+%   \def\L{3.8}       % distance between walls
+%   \def\H{5.4}       % total wall height
+%   \def\h{2.8}       % plane wave height
+%   \def\t{0.15}      % wall thickness
+%   \def\a{1.15}      % slit distance
+%   \def\d{0.20}      % slit size
+%   \def\N{21}        % number of waves
+%   \def\lambd{0.20}  % wavelength
+%   \def\R{\N*\lambd} % wave radius
+%   \def\Nlines{3}    % number of nodal lines
+%   \def\A{1.6}       % amplitude
+%   %\def\r{0.06}      % point source radius
+%   %\def\nmax{10}
+%   \def\nsamples{100}
+%   \def\ang{62}
+  
+%   \begin{scope}
+%     \clip (-\t/2,-\H/2) rectangle (\L,\H/2);
+%     %\clip (-\t/2,0.7*\a) -- (0.6*\L,\H/2) -- (\L,\H/2) --
+%     %      (\L,-\H/2) -- (0.6*\L,-\H/2) -- (-\t/2,-0.7*\a) -- cycle;
+    
+%     % NODAL LINES
+%     \draw[nodal]
+%       (0.08*\N*\lambd,0) -- (1.06*\R,0);
+%     \coordinate (NP0) at (\L,0);  % to avoid "Dimension too large error"
+%     \foreach \dn [evaluate={
+%                    \f=\a/\lambd;
+%                    \nmin=2.5+0.2*\dn; %0.501*(-\dn+\f)
+%                    \nmax=10; %(NP0)
+%                    \c=int(\dn<\f);
+%                    \y=\L/sqrt((\a/(\lambd*\dn))^2-1);
+%                  }] in {1,...,\Nlines}{
+%       \coordinate (NP+\dn) at (\L,\y);  % to avoid "Dimension too large error"
+%       \coordinate (NP-\dn) at (\L,-\y); % to avoid "Dimension too large error"
+%       \ifnum\c=1
+%         \draw[nodal,variable=\n,samples=\nsamples,smooth]
+%           plot[domain=\nmin:\nmax] ({xnode(\n,\dn,\lambd,\f)},{ynode(\n,\dn,\lambd,\f)})
+%           -- (NP+\dn);
+%         \draw[nodal,variable=\n,samples=\nsamples,smooth]
+%           plot[domain=\nmin:\nmax] ({xnode(\n,\dn,\lambd,\f)},{-ynode(\n,\dn,\lambd,\f)})
+%           -- (NP-\dn);
+%       \fi
+%     }
+    
+%     % WAVES
+%     \foreach \i [evaluate={\R=\i*\lambd;}] in {1,...,\N}{
+%       \ifodd\i
+%         \draw[myblue,line width=0.8] (0,\a/2)++(\ang:\R) arc (\ang:-\ang:\R);
+%         \draw[myred,line width=0.8] (0,-\a/2)++(\ang:\R) arc (\ang:-\ang:\R);
+%       \else
+%         \draw[myblue!80,line width=0.1] (0,\a/2)++(\ang:\R) arc (\ang:-\ang:\R);
+%         \draw[myred!80,line width=0.1] (0,-\a/2)++(\ang:\R) arc (\ang:-\ang:\R);
+%       \fi
+%     }
+%   \end{scope}
+  
+%   % PLANE WAVES
+%   \foreach \i [evaluate={\x=-\i*\lambd;}] in {0,...,5}{
+%     \ifodd\i
+%       \draw[myblue,line width=0.8] (\x,-\h/2) -- (\x,\h/2);
+%     \else
+%       \draw[myblue,line width=0.1] (\x,-\h/2) -- (\x,\h/2);
+%     \fi
+%   }
+  
+%   % WALL
+%   \fill[wall]
+%     (\t/2,\a/2-\d/2) rectangle (-\t/2,-\a/2+\d/2)
+%     (\t/2,\a/2+\d/2) rectangle (-\t/2,\H/2)
+%     (\t/2,-\a/2-\d/2) rectangle (-\t/2,-\H/2)
+%     (\L,-\H/2) rectangle (\L+\t,\H/2);
+  
+%   % SHADES
+%   \begin{scope}[shift={(1.08*\L,0)}]
+%     \def\yz{\L/sqrt((\a/\lambd)^2-1)} % m = +- 1/2
+%     \def\yZ{\L/sqrt((\a/\lambd/2)^2-1)} % m = +- 1
+%     \clip (0,-\H/2) rectangle (1.1*\A,\H/2);
+%     \fill[white] (0,-\H/2) rectangle++ (\A,\H); % to fill seams
+%     \foreach \i [evaluate={\n=0.5*\i;\yn=\L/sqrt((\a/(2*\lambd*\n))^2-1);
+%                  }] in {1,...,\Nlines}{
+%       \ifodd\i % if even
+%         \fill[myshadow] (0,{-\yn-0.1}) rectangle++ (\A,0.2); % to fill seams
+%         \fill[myshadow] (0,{ \yn-0.1}) rectangle++ (\A,0.2); % to fill seams
+%       \fi
+%     }
+%     \path[left color=myshadow,right color=myshadow,middle color=white,shading angle={180}]
+%       (0,{-\yz}) rectangle (\A,{\yz});
+%     \foreach \i [evaluate={
+%                   \n=0.5*\i;
+%                   \m=0.5*(\i+1);
+%                   \yn=\L/sqrt((\a/(2*\lambd*\n))^2-1);
+%                   \ym=\L/sqrt((\a/(2*\lambd*\m))^2-1);
+%                   \dang=mod(\i,2)*180;
+%                  }] in {1,...,\Nlines}{
+%       \path[left color=myshadow,right color=white,shading angle={\dang}]
+%         (0,\yn) rectangle (\A,\ym);
+%       \path[left color=myshadow,right color=white,shading angle={180+\dang}]
+%         (0,-\yn) rectangle (\A,-\ym);
+%     }
+%   \end{scope}
+  
+%   % INTENSITY
+%   \begin{scope}[shift={(1.1*\L+1.1*\A,0)}]
+%     \draw[->,thick] (-0.08*\A,0) -- (1.3*\A,0) node[right=-2] {$\expval{I}$}; % I axis
+%     \draw[->,thick] (0,-0.52*\H) -- (0,0.54*\H) node[right] {$y$}; % y axis
+%     \draw[nodal] (NP0) --++ (0.15*\L+2.1*\A,0); % green nodal lines
+%     \foreach \i [evaluate={\y=\L/sqrt((\a/(\lambd*\i))^2-1)}] in {1,...,\Nlines}{ % green nodal lines
+%       \draw[nodal] (NP+\i) --++ ({0.15*\L+1.1*\A+\A*intensity(\y,\lambd,\a,\L)},0);
+%       \draw[nodal] (NP-\i) --++ ({0.15*\L+1.1*\A+\A*intensity(\y,\lambd,\a,\L)},0);
+%     }
+%     \draw[myred,thick,variable=\y,samples=\nsamples,smooth,domain=-\H/2:\H/2]
+%       plot({\A*intensity(\y,\lambd,\a,\L)},\y);
+%     \foreach \i [evaluate={ % ticks
+%         \modd=\i; %int(\i);
+%         \meven=int(\i-1);
+%         \y=\L/sqrt((\a/(\lambd*\i))^2-1);
+%                 }] in {1,...,\Nlines}{
+%       \ifodd\i
+%         \tick{0,-\y}{180} node[right=0,scale=0.85] {$m=-\frac{\modd}{2}$};
+%         \tick{0,\y}{180} node[right=0,scale=0.85] {$m=+\frac{\modd}{2}$};
+%       \else
+%         \tick{0,-\y}{180} node[right=0,scale=0.85] {$m=-\meven$};
+%         \tick{0,\y}{180} node[right=0,scale=0.85] {$m=+\meven$};
+%       \fi
+%     }
+%   \end{scope}
+  
+% \end{tikzpicture}
+
+% Schrödinger
+% -----------
+% The wave equation, Schrödinger equation has a general form 
+% \begin{align}
+%     i \hbar \frac{\partial}{\partial t} | \Psi \rangle &= \hat{H} | \Psi \rangle \ ,
+%     \label{eq:schrodinger_general}
+% \end{align}
+% where $i$ is the imaginary number, and $\hbar$ is Plancks constant. Here $| \Psi \rangle$ 
+% is the quantum state and $\hat{H}$ is a Hamiltonian operator. % which represent the energy of the system 
+% For two-dimensional position space, the quantum state can be expressed using the 
+% time-dependent complex-valued wave function $\Psi (x, y, t)$. The square modulus of the wave function $|\psi|^{2}$, predicts probability of finding the particle at position $(x, y)$ at time t. 
+% % Segue to Born
+% The modulus of the wave function, is related to the probability density function 
+% \begin{align}
+%     p(x, y \ | \ t) &= |\Psi(x, y, t)|^{2} = \Psi^{*}(x, y, t) \Psi(x, y, t)
+%     \label{eq:born_rule}
+% \end{align}
+% using the Born rule, where $^{*}$ denotes the complex conjugated wave function. 
+% When the potential is time-independent, the Schrödinger equation can be expressed as 
+% \begin{align*}
+%     i \hbar \frac{\partial}{\partial t} \Psi (x, y, t) &= - \frac{\hbar^{2}}{2m} \bigg( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \bigg) \Psi (x, y, t) \\
+%     & \quad + V(x, y, t) \Psi (x, y, t) \numberthis \ .
+%     \label{eq:schrodinger_special}
+% \end{align*}
+% When we scale Schrödinger equation by the dimensionful variables, we are left with 
+% a wave function $u$, potential $v$ and the dimensionless equation
+% \begin{align}
+%     i \frac{\partial u}{\partial t} &= - \frac{\partial^{2} u}{\partial x^{2}} - \frac{\partial^{2} u}{\partial y^{2}} + v(x, y) u \ .
+%     \label{eq:schrodinger_dimensionless}
+% \end{align}
+
+% Crank-Nicolson
+% --------------
+% To evaluate the position of a single particle, we have to consider partial differential 
+% equations (PDE). To solve these numerically, we have to discretize Equation \eqref{eq:schrodinger_dimensionless}. 
+% We use the Crank-Nicolson method (CN), which combines the forward and backward finite 
+% difference method. 
+% Include the general $\theta$-method 
+% Using the $\theta$-rule \footnote{For $\theta \in [0, 1]$ we can derive Forward Euler using $\theta = 1$, Backward Euler using $\theta = 0$, and $\theta = 1/2$ gives Crank-Nicolson} 
+% with $\theta = 1/2$, CN can be written as 
+% \begin{align}
+%     \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \ .
+%     \label{eq:crank_nicolson}
+% \end{align}
+% To simplify notation and avoid confusion of indices with the imaginary number $i$, 
+% we have used the notation $\ivec, \jvec$ in subscript to indicate the commonly named indices $i, j$ 
+% in x- and y-direction. In addition, the superscript $n, n+1$ indicate position in time. 
+% We use CN to derive the discretized Schrödinger equation 
+% \begin{align*}
+%     & u_{\ivec, \jvec}^{n+1} - \frac{i \Delta t}{2 \Delta x^{2}} \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
+%     & - \frac{i \Delta t}{2 \Delta y^{2}} \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+%     &= u_{\ivec, \jvec}^{n} + \frac{i \Delta t}{2 \Delta x^{2}} \big[ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n} \big] \\
+%     & \quad + \frac{i \Delta t}{2 \Delta y^{2}} \big[ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n} \big] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} \numberthis \ .
+%     \label{eq:schrodinger_discretized} 
+% \end{align*}
+% The full derivation of both Equation \eqref{eq:crank_nicolson} and Equation \eqref{eq:schrodinger_discretized} 
+% can be found in Appendix \ref{ap:crank_nicolson}. 
+
+% Implementation
+The implementation of CN is simplified usein Dirichlet boundary conditions, which 
+can be found in Table \ref{tab:dirichlet}. In addition, we use Gaussian wave packets 
+for the initial wave function 
+\begin{align}
+    u(x, y, t=0) = e^{- \frac{(x-x_{c})^{2}}{2 \sigma_{x}^{2}} - \frac{(y-y_{c})^{2}}{2 \sigma_{y}^{2}} + ip_{x}x + ip_{y}y}
+\end{align}
+\begin{algorithm}[H]
+    \caption{Crank-Nicolson scheme}
+    \label{algo:cn_scheme}
+    \begin{algorithmic}
+            \Procedure{Crank-Nicolson}{$args$}
+            \State Insert pseudo code $var \leftarrow \text{ some text}$
+            \EndProcedure
+    \end{algorithmic}
+\end{algorithm}
+
+
+\section{Notes}\label{sec:notes}
+\subsection*{Introduction - draft 2}
+In classical mechanics we study the kinematics and dynamics of physical objects, 
+ignoring their intrinsic properties for simplicity. It allows us to describe the 
+forces acting on an object as well as the motion of the object. We can describe 
+a planets orbital movement \cite{britannica:2023:kepler}, calculate the ... necessary 
+to launch satellites into orbit, or simply figure out where a ball is going to land 
+when you throw it... However, when want to study an object at a microscopic level, 
+e.g. a single atom, classical mechanics falls short. 
+
+% Thomas Young first performed the double-slit experiment in 1801, to demonstrate 
+% the principle of interference of light \cite{britannica:2023:young}, postulating 
+% light as waves. In the study of blackbodies, scientists were not able to describe 
+% the radiated intensity of increased frequencies using classical mechanichs, as they 
+% contradicted the principle of conservation of energy \cite{britannica:1998:planck}. 
+% Max Planck assumed that the radiated energy consist of discrete values, or quanta, 
+% to describe the peak in radiated energy. 
+
+Light as particles -> waves -> particles/packets
+Erwin Schrödinger wanted to find a mathematical description of the wave characteristics 
+of matter, supporting the wave-particle idea. He postulated a function which varies 
+with position, where the function swuared can be interpreted as the probability 
+of finding an electron at a given position. 
+
+\subsection*{Introduction - draft 1}
+In classical mechanics, we study the motion of particles at a macroscopic ,using 
+physical laws such as Newtons second law. These laws and concepts were used to solve 
+the motion of the planets. 
+
+However, when we zoom in on a single atom, at a microscopic level, classical mechanics 
+are no longer sufficient in describing the motion of particles. That is, elementary 
+particles do not behave in a way which can be described using classical mechanics. 
+Introducing quantum mechanics, where a quantum is the minimum amount of something, 
+such as an elementary particle of a given field. 
+
+Beginning with the study of blackbodies and radiation, radiated intensity as a 
+function of frequency, contradicting the principle of conservation of energy. Max 
+Planck derived an equation Plancks Law, electromagnetic energy takes the form of 
+quanta (discrete packets now known as photons). Energy of the photon is equal to 
+plancks constant times the freq of light.
+
+Thomas Young first performed the double-slit experiment in 1801, to demonstrate 
+the principle of interference of light \cite{britannica:2023:young}, postulating 
+light as waves.  
+
+In the 1860s James Clerk Maxwell predicted that light waves consisted of coupled 
+electric and magnetic fields.
+
+Erwin Schrödinger wanted to find a wave equation for matter when classical mechanics 
+fell short, the wave function. The wave function squared can be interpreted as the 
+probability of finding an electron at a given position. Schrödingers equation can 
+be solved for the hydrogen atom, wave equation of the energy levels for a hydrogen 
+atom. The time-dependent equation also shows how a quantum state evolves with time 
+\cite[p. 81]{wu:2023:quantum}.
+
+\subsection*{Research material}
+Simulating the double-slit experiment requires solving partial differential equations, 
+when the number of spatial variables increases so does the complexity of the problem. 
+
+Crank-Nicolson set out to find a numerical method to solve the diffusion equation, 
+with a higher stability than the explicit and implicit schemes. Combining both schemes, 
+they developed a finite difference method where a solution is stable for all dx and dt.
+
+To solve equations, such as Schrödingers, numerically, a discretized equation is 
+derived using difference methods of explicit and implicit scheme.
+
+Distribution of particle position.
+The wave function can be interpreted in several ways, not unity in agreement, does it apply to one electron, or a system of electrons. Use to find probabilities of measurements. Probability (squared) of finding the electron at a given position, need to do measurements many times for statistical reasons...
+
+energy defined in quanta, photoelectric effect
+
+Applying the wave-particle duality to matter, allows us to analyze the behavior of elementary particles more accurately.
+
+Heisenberg uncertainty principle
+
+diffraction and interference, describing waves through double slit experiment (Thomas Young). Diffraction pattern where light interfere constructively result in bright spots, whereas destructive interference result in dark spots. 
+
+Path difference: distance between the center of each slit times sine of the angle between the point on the screen and the slits.
+
+Progress in quantum mechanics without knowing how it works!
+- Quantum is the minimum amount of something, elementary particle of a given field, smallest unit of excitation in a fundamental field (inject energy).
+- Planck: father of quantum mechanics
+Quantum mechanics (MIT open course ware):
+- superposition principle: an object can simultaneously appear in two different places or have two different velocities \cite[p. 3]{wu:2023:quantum}
+- quantum randomness: the outcome of a measurement is random when a particle is in a superposition of two positions.
+- electron color: black or white, hardness: hard or soft
+- Identical particles: an absolute in quantum mechanics, two particles can't be distinguished from each other. 
+- Wave-Particle duality: Electrons are waves, and can form standing waves around protons.
+- Schrödinger equation: wave equation of the energy levels for a hydrogen atom. It also shows how a quantum state evolves with time \cite[p. 81]{wu:2023:quantum}. 
+\begin{align}
+    i\hbar \frac{\partial}{\partial t} \Psi (x, y, t) &= - \frac{\hbar^{2}}{2m} \bigg( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \bigg) \Psi (x, y, t) + V(x, y, t) \Psi (x, y, t)
+\end{align}
+- $\psi$ is a wave function and describes the quantum state of a particle at time t. The square modulus of the wave function $|\psi|^{2}$, predicts probability of finding the particle at position $(x, y)$ at time t. 
+- lhs: describe the rate of change of a wave function over time t. 
+- rhs: describe the variation of the wave function in space.
+- eq: describe the relation between the lhs and rhs, formulation of wave-particle duality.
+
+Partial differential equations
+- problems with many variables, constrained by boundary conditions and initial values.
+
+Diffusion equation: describe the density of a quantity as a function of time, where the flux density obeys the Gauss-Green theorem (div = 0).
+Crank-Nicolson: combines the implicit and explicit methods in solving a pde, rewriting the pde as a set of mat-vec multiplications.
+Originally the method was derived for the diffusion 
+equation, however, it can be used for the wave equation as well.
+
+Wave equation in two dimensions (lec. p. 322), discretize position and time, to obtain 
+\begin{align*}
+    u_{\ivec, \jvec}^{n+1} &- r [ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1}] - r [ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1}] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+    &= u_{\ivec, \jvec}^{n} + r [ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n}] + r [ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n}] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n}
+\end{align*}
+Goal: Simulate the two-dimensional time-dependent schrödinger equation, to study a double-slit-in-a-box setup.
+% Introduction
+The double-slit experiment, interference using slits observe result on canvas. An experiment by Thomas Young to demonstrate the wave-particle duality of light. Electrons are waves as well as particles. 
+
+When we consider an elementary particle, classical physics are not sufficient in describing the particles position in time. We introduce a probability density function for detecting the particle at a given position. 
+
+We can study the quantum state of a particle in two dimensions, using Schrödinger's equation and the wave function. The kinetic energy are found from the partial derivative of the wave function in respect to the position varables. To solve this numerically, we use the crank-nicolson method to discretize position and time, which combines the explicit and implicit method in solving pdes.
+\end{document}

latex/references.bib

@@ -0,0 +1,231 @@
+# Double-slit and physics
+@article{murray:2020:double_slits,
+  author          = {Andrew Murray},
+  journal         = {Physics World},
+  number          = {2},
+  title           = {Double slits with single atoms},
+  volume          = {33},
+  year            = {2020}
+}
+
+@misc{mit:2004:physics,
+  author          = {Sen-ben Liao and Peter Dourmashkin and John W. Belcher},
+  title           = {Course notes in Physics 8.02 at MIT},
+  year            = {2004},
+  url             = {https://web.mit.edu/8.02t/www/802TEAL3D/visualizations/coursenotes/modules/guide14.pdf},
+  urldate         = {2023-12-22}
+}
+
+@misc{britannica:1999:light,
+  author          = {The Editors of Encyclopaedia Britannica},
+  title           = {Light},
+  publisher       = {Britannica},
+  url             = {https://www.britannica.com/science/light},
+  urldate         = {2023-12-15},
+}
+
+@article{sinha:2010:multi_order_inference,
+  author          = {Urbasi Sinha and Christophe Couteau and Thomas Jennewein and Raymond Laflamme and Gregor Weihs},
+  journal         = {Science},
+  title           = {Ruling Out Multi-Order Interference in Quantum Mechanics},
+  volume          = {329},
+  number          = {5990},
+  pages           = {418--421},
+  year            = {2010},
+  doi             = {10.1126/science.1190545},
+  url             = {https://www.science.org/doi/abs/10.1126/science.1190545},
+  urldate         = {2023-12-16}
+}
+
+@article{young:1804:double_slit,
+  author          = {Thomas Young},
+  title           = {I. The Bakerian Lecture. Experiments and calculations relative to physical optics},
+  journal         = {Philosophical Transactions of the Royal Society of London},
+  volume          = {94},
+  pages           = {1--16},
+  year            = {1804},
+  doi             = {10.1098/rstl.1804.0001},
+  url             = {https://royalsocietypublishing.org/doi/abs/10.1098/rstl.1804.0001},
+  urldate         = {2023-12-16},
+}
+
+@book{benacquista:2018:classical_mechanics,
+  author         = {Matthew J. Benacquista and Joseph D. Romano},
+  title          = {Classical Mechanics},
+  publisher      = {Springer International Publishing : Imprint: Springer},
+  year           = {2018},
+  edition        = {1}
+}
+
+@book{griffiths:2018:quantum_mechanics,
+  author         = {David J. Griffiths and Darrell F. Schroeter},
+  title          = {Introduction to quantum mechanics},
+  publisher      = {Cambridge University Press},
+  year           = {2018},
+  edition        = {3}
+}
+
+@book{springer:2018:compendium_quantum_physics,
+  author    = {Daniel Greenberger and Klaus Hentschel and Friedel Weinert},
+  title     = {Compendium of Quantum Physics},
+  publisher = {Springer Berlin, Heidelberg},
+  year      = {2009},
+  edition   = {1}
+}
+
+@article{key,
+  author          = {Asher Peres and Daniel R. Terno},
+  title           = {Quantum Information and Relativity Theory},
+  journal         = {Reviews of Modern Physics},
+  number          = {1},
+  volume          = {76},
+  doi             = {10.1103/revmodphys.76.93},
+  year            = {2004},
+  url             = {http://dx.doi.org/10.1103/RevModPhys.76.93}
+}
+
+# Schrodinger and Crank-Nicolson
+@book{wu:2023:quantum,
+    author = {Biao Wu and Ying Hu},
+    title = {Quantum Mechanics: A Concise Introduction},
+    publisher = {Springer Singapore},
+    year = {2023},
+    edition = {1},
+}
+
+@misc{britannica:2023:young,
+    author = {Glenn Stark},
+    title = {Young’s double-slit experiment},
+    publisher = {Britannica},
+    year = {2023},
+    url = {https://www.britannica.com/science/light/Youngs-double-slit-experiment},
+    note = {Last accessed 2023-12-12}
+}
+
+@misc{britannica:2023:kepler,
+    author = {The Editors of Encyclopaedia Britannica}, 
+    title = {Kepler’s laws of planetary motion},
+    publisher = {Britannica},
+    url = {https://www.britannica.com/science/Keplers-laws-of-planetary-motion},
+    note = {Last accessed 2023-12-13}
+}
+
+@misc{britannica:1998:planck,
+    author = {The Editors of Encyclopaedia Britannica}, 
+    title = {Planck’s radiation law},
+    publisher = {Britannica},
+    url = {https://www.britannica.com/science/Plancks-radiation-law},
+    urldate = {2023-12-13}
+}
+
+@misc{britannica:1998:newton,
+    author = {The Editors of Encyclopaedia Britannica}, 
+    title = {Newton’s laws of motion},
+    publisher = {Britannica},
+    url = {https://www.britannica.com/science/Newtons-laws-of-motion},
+    urldate = {2023-12-13}
+}
+
+@article{crank:1947:numerical,
+    author = {John Crank and Phyllis Nicolson},
+    title = {A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type},
+    journal = {Mathematical proceedings of the Cambridge Philosophical Society},
+    doi = {10.1017/S0305004100023197},
+    year = {1947},
+    volume = {43},
+    number = {1},
+    pages = {50--67}
+}
+
+# Math and statistics
+@book{lindstrom:2016:kalkulus,
+  author    = {Tom Lindstrøm},
+  title     = {Kalkulus},
+  publisher = {Universitetsforlaget},
+  year      = {2016},
+  edition   = {4}
+}
+@book{gjevik:2014:appendix,
+  author    = {Bjørn Gjevik and Morten Wang Fagerland},
+  title     = {Feltteori og vektoranalyse},
+  publisher = {Farleia Forlag},
+  year      = {2018}
+}
+@book{springer:2012:modernstat,
+  author    = {Jay L. Devore and Kenneth N. Berk},
+  title     = {Modern Mathematical Statistics with Application},
+  publisher = {Cham: Springer International Publishing AG},
+  year      = {2021},
+  edition   = {3}
+}
+
+# Python libraries
+@article{hunter:2007:matplotlib,
+  author    = {Hunter, J. D.},
+  title     = {Matplotlib: A 2D graphics environment},
+  journal   = {Computing in Science \& Engineering},
+  volume    = {9},
+  number    = {3},
+  pages     = {90--95},
+  abstract  = {Matplotlib is a 2D graphics package used for Python for
+               application development, interactive scripting, and publication-quality
+               image generation across user interfaces and operating systems.},
+  publisher = {IEEE COMPUTER SOC},
+  doi       = {10.1109/MCSE.2007.55},
+  year      = 2007
+}
+@article{waskom:2021:seaborn,
+  author    = {Michael L. Waskom},
+  title     = {seaborn: statistical data visualization},
+  journal   = {Journal of Open Source Software},
+  volume    = {6},
+  number    = {60},
+  publisher = {The Open Journal},
+  year      = {2021},
+  doi       = {10.21105/joss.03021},
+  url       = {https://doi.org/10.21105/joss.03021},
+  pages     = {3021}
+}
+
+# C++ libraries
+@misc{openmp:2018,
+  author  = {OpenMP},
+  title   = {OpenMP Application Programming Interface},
+  year    = {2018},
+  url     = {https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5.0.pdf},
+  urldate = {2023-10-22}
+}
+@inproceedings{gabriel:2004:open_mpi,
+  author    = {Edgar Gabriel and Graham E. Fagg and George Bosilca
+               and Thara Angskun and Jack J. Dongarra and Jeffrey
+               M. Squyres and Vishal Sahay and Prabhanjan Kambadur
+               and Brian Barrett and Andrew Lumsdaine and Ralph
+               H. Castain and David J. Daniel and Richard L. Graham
+               and Timothy S. Woodall},
+  title     = {Open {MPI}: Goals, Concept, and Design of a Next Generation {MPI} Implementation},
+  booktitle = {Proceedings, 11th European PVM/MPI Users' Group Meeting},
+  address   = {Budapest, Hungary},
+  pages     = {97--104},
+  year      = {2004},
+  month     = {September}
+}
+
+# Miscellaneous things
+
+
+
+# Profiler
+@online{scalasca,
+  author  = {M. Geimer and F. Wolf and B.J.N. Wylie and E. Abraham and D. Becker and B. Mohr},
+  title   = {Scalasca},
+  url     = {https://www.scalasca.org/scalasca/about/about.html},
+  urldate = {2023-10-24},
+  note    = {Tool to support performance optimization of parallel programs, measuring and analyzing runtime behavior.}
+}
+@online{scorep,
+  title   = {Score-P: the Scalable Performance Measurement Infrastructure for Parallel Codes},
+  url     = {https://perftools.pages.jsc.fz-juelich.de/cicd/scorep/tags/latest/html/},
+  urldate = {2023-10-24},
+  note    = {Tool suite for profiling and event tracing.}
+}
+

latex/schrodinger_simulation.tex

@@ -0,0 +1,133 @@
+\documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-1}  % defines the basic parameters of the document
+% 
+% For preview: skriv i terminal: latexmk -pdf -pvc filnavn
+% If you want a single-column, remove "reprint"
+
+% Silence warning of revtex4-1
+\usepackage{silence}
+\WarningFilter{revtex4-1}{Repair the float}
+
+% Allows special characters (including æøå)
+\usepackage[utf8]{inputenc}
+\usepackage{fontawesome}
+\usepackage{multirow}
+\usepackage[table]{xcolor}
+% \rowcolors{2}{gray!15}{white}
+% \usepackage{tabularx}
+% \usepackage[english]{babel}
+
+%% Note that you may need to download some of these packages manually, it depends on your setup.
+%% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
+
+\usepackage{physics,amssymb}  % mathematical symbols (physics imports amsmath)
+\usepackage[no-test-for-array]{nicematrix}
+
+\usepackage{graphicx}         % include graphics such as plots
+\graphicspath{.images/}
+\usepackage{xcolor}           % set colors
+\usepackage{hyperref}         % automagic cross-referencing
+\usepackage{listings}         % display code
+% \usepackage{subfigure}        % imports a lot of cool and useful figure commands
+\usepackage{subcaption}
+% \usepackage{float}
+%\usepackage[section]{placeins}
+\usepackage{algorithm}
+\usepackage[noend]{algpseudocode}
+\usepackage{tikz}
+% \usepackage{pgfplots}
+% \pgfplotsset{compat=1.18}
+% \usetikzlibrary{3d}
+\usepackage{xcolor}
+% \usepackage{etoolbox} %ifthen
+% \usetikzlibrary{calc}
+% \usetikzlibrary{arrows,arrows.meta}
+% \usetikzlibrary{decorations.markings}
+% \usetikzlibrary{angles,quotes}
+% \usetikzlibrary{fadings}
+% \tikzset{>=latex}
+% \colorlet{wall}{blue!30!black}
+% \colorlet{myblue}{blue!70!black}
+% \colorlet{myred}{red!70!black}
+% \colorlet{mydarkred}{red!50!black}
+% \colorlet{mylightgreen}{green!60!black!70}
+% \colorlet{mygreen}{green!60!black}
+% \colorlet{myredgrey}{red!50!black!80}
+% \colorlet{myshadow}{blue!30!black!90}
+% \tikzstyle{wave}=[myblue,thick]
+% \tikzstyle{mydashed}=[black!70,dashed,thin]
+% \tikzstyle{mymeas}=[{Latex[length=3,width=2]}-{Latex[length=3,width=2]},thin]
+% \tikzstyle{mysmallarr}=[-{Latex[length=3,width=2]}]
+% \usetikzlibrary{quantikz}
+% defines the color of hyperref objects
+% Blending two colors:  blue!80!black  =  80% blue and 20% black
+\hypersetup{ % this is just my personal choice, feel free to change things
+    colorlinks,
+    linkcolor={red!50!black},
+    citecolor={blue!50!black},
+    urlcolor={blue!80!black}}
+
+% Biblio stuff
+% \def\biblio{\bibliographystyle{plain}\bibliography{../references/references}}
+
+\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}}
+\newcommand{\rvline}{\hspace*{-\arraycolsep}\vline\hspace*{-\arraycolsep}}
+
+% Defines indices i and j to avoid confusion with imaginary i
+\newcommand{\ivec}{\hat{\imath}}
+\newcommand{\jvec}{\hat{\jmath}}
+
+\usepackage{xr}
+\usepackage{subfiles}
+% \externaldocument[M-]{\subfix{main}}
+
+\begin{document}
+
+\title{Simulating the Schrödinger wave equation using the Crank-Nicolson method in 2+1 dimensions}  % self-explanatory
+\author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen \\ \faGithub \, \url{https://github.uio.no/FYS3150-G2-2023/Project-4}} % self-explanatory
+\date{\today}                             % self-explanatory
+\noaffiliation                            % ignore this, but keep it.
+
+% Abstract
+\subfile{sections/abstract}
+
+\maketitle
+
+% Introduction
+\subfile{sections/introduction}
+
+% Methods
+\subfile{sections/methods}
+
+% Results
+\subfile{sections/results}
+
+% Conclusion
+% \subfile{sections/conclusion}
+
+% Notes
+% \subfile{draft}
+
+\clearpage
+\newpage
+
+% Appendix
+\subfile{sections/appendices}
+
+\clearpage
+\onecolumngrid
+
+\bibliographystyle{unsrtnat}
+\bibliography{references}
+
+
+\end{document}
+
+% Methods
+% P1: Theory, imag i = i, index i, j = \hat{i}, \hat{j}
+% P2:
+% 
+
+% Results
+% P7:
+% P8:
+% P9:

latex/sections/abstract.tex

@@ -0,0 +1,18 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\begin{abstract}
+We have simulated the two-dimensional time-dependent Schrödinger equation, to study 
+variations of the double-slit experiment. To solve the partial differential equations 
+we have applied the Crank-Nicolson scheme in 2+1 dimensions, to derive a discretized 
+equation. In addition, we have used Dirichlet boundary conditions to express the 
+equation in matrix form and solve it using the sparse matrix solver \verb|superlu|. 
+Our implementation, and choice of solver method, resulted in conserved total 
+probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}=1$ for both the single and 
+double slit setup. To illustrate the time evolution of the probability function, 
+we created colormap plots at time steps $t = [0, 0.001, 0.002]$. We also included 
+separate plots for each time step of Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$.
+In addition, we determined the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, 
+for single-, double- and triple-slit.
+\end{abstract}
+\end{document}

latex/sections/appendices.tex

@@ -0,0 +1,97 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\appendix
+\section{The Crank-Nicholson method}\label{ap:crank_nicolson}
+The Crank-Nicolson (CN) approach consider both the forward difference, an explicit scheme,
+\begin{equation*}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} = F_{\ivec, \jvec}^{n} \ ,
+\end{equation*}
+and the backward difference, an implicit scheme, 
+\begin{equation*}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} = F_{\ivec, \jvec}^{n+1} \ .
+\end{equation*}
+The result is a linear combination of the explicit and implicit scheme, given by 
+\begin{align*}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \theta F_{\ivec, \jvec}^{n+1} + (1 - \theta) F_{\ivec, \jvec}^{n} \ .
+\end{align*}
+The parameter $\theta$ is introduced for a general approach, however, for CN $\theta = 1/2$. 
+\begin{align*}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \\
+\end{align*}
+
+We need the first derivative in respect to both time and position, as well as the second derivative in respect to position. Taylor expanding will result in a discretized version, assume this is known...
+
+Schrödinger contain $i$ at the lhs, factor it as
+\begin{align*}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2i} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \\
+    &= -\frac{i}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] & \text{, where $\frac{1}{i} = -i$}
+\end{align*}
+
+Using Equation \eqref{eq:schrodinger_dimensionless}, we get
+\begin{align*}
+    u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n} & -\frac{i \Delta t}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \\
+    &= -\frac{i \Delta t}{2} \bigg[ - \frac{u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1}}{2 \Delta x^{2}} \\
+    & \quad - \frac{u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1}}{2 \Delta y^{2}} + \frac{1}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+    & \quad - \frac{u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n}}{2 \Delta x^{2}} \\
+    & \quad - \frac{u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n}}{2 \Delta y^{2}} + \frac{1}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} \bigg]
+\end{align*}
+We rewrite the expression and gather all terms containing the $n+1$ time step on 
+the left hand side, and the terms containing $n$ time step on the right hand side. 
+\begin{align*}
+    & u_{\ivec, \jvec}^{n+1} - \frac{i \Delta t}{2 \Delta x^{2}} \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
+    & - \frac{i \Delta t}{2 \Delta y^{2}} \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+    &= u_{\ivec, \jvec}^{n} + \frac{i \Delta t}{2 \Delta x^{2}} \big[ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n} \big] \\
+    & \quad + \frac{i \Delta t}{2 \Delta y^{2}} \big[ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n} \big] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} 
+\end{align*}
+In addition, since we will use an equal step size $h$ in both $x$ and $y$ direction, 
+we can use 
+\begin{align*}
+    \frac{i \Delta t}{2 \Delta h^{2}} = \frac{i \Delta t}{2 \Delta x^{2}} = \frac{i \Delta t}{2 \Delta y^{2}} \ ,
+\end{align*}
+and define 
+\begin{align*}
+    r \equiv \frac{i \Delta t}{2 \Delta h^{2}}
+\end{align*}
+Now, the discretized Schrödinger equation can be written as 
+\begin{align*}
+    & u_{\ivec, \jvec}^{n+1} - r \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
+    & - r \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+    &= u_{\ivec, \jvec}^{n} + r \big[ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n} \big] \\
+    & \quad + r \big[ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n} \big] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} \ .
+\end{align*}
+
+
+\section{Matrix structure}\label{ap:matrix_structure}
+For $u$ vector of length $(M-2) = 3$, the matrices $A$ and $B$ have size 
+$(M-2)^{2} \times (M-2)^{2} = 9 \times 9$ given by 
+\begin{align*}
+    A = 
+    \begin{bmatrix}
+    a_{0} &  -r  &  0   &  -r   &   0   &  0    &  0   &  0    &  0    \\
+    -r  &  a_{1} &  -r  &  0    &   -r  &  0    &  0   &  0    &  0    \\
+    0   &  -r  &  a_{2} &  0    &   0   &  -r   &  0   &  0    &  0    \\
+    -r  &  0   &  0   &  a_{3}  &   -r  &  0    &  -r  &  0    &  0    \\
+    0   &  -r  &  0   &  -r   &   a_{4} &  -r   &  0   &  -r   &  0    \\
+    0   &  0   &  -r  &  0    &   -r  &  a_{5}  &  0   &  0    &  -r   \\
+    0   &  0   &  0   &  -r   &   0   &  0    &  a_{6} &  -r   &  0    \\
+    0   &  0   &  0   &  0    &   -r  &  0    &  -r  &  a_{7}  &  -r   \\
+    0   &  0   &  0   &  0    &   0   &  -r   &  0   &  -r   &  a_{8}  \\
+    \end{bmatrix}
+\end{align*}
+
+\begin{align*}
+    B = 
+    \begin{bmatrix}
+    b_{0} &  r   &  0   &  r    &   0   &  0    &  0   &  0    &  0    \\
+    r   &  b_{1} &  r   &  0    &   r   &  0    &  0   &  0    &  0    \\
+    0   &  r   &  b_{2} &  0    &   0   &  r    &  0   &  0    &  0    \\
+    r   &  0   &  0   &  b_{3}  &   r   &  0    &  r   &  0    &  0    \\
+    0   &  r   &  0   &  r    &   b_{4} &  r    &  0   &  r    &  0    \\
+    0   &  0   &  r   &  0    &   r   &  b_{5}  &  0   &  0    &  r    \\
+    0   &  0   &  0   &  r    &   0   &  0    &  b_{6} &  r    &  0    \\
+    0   &  0   &  0   &  0    &   r   &  0    &  r   &  b_{7}  &  r    \\
+    0   &  0   &  0   &  0    &   0   &  r    &  0   &  r    &  b_{8}  \\
+    \end{bmatrix}
+\end{align*}
+\end{document}

latex/sections/conclusion.tex

@@ -0,0 +1,25 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\section{Conclusion}\label{sec:conclusion}
+% Rewrite this section to differ from the abstract
+We have simulated the two-dimensional time-dependent Schrödinger equation, to study 
+variations of the double-slit experiment. To solve the partial differential equations 
+we have applied the Crank-Nicolson scheme in 2+1 dimensions, and derived a discretized 
+equation. In addition, we have used Dirichlet boundary conditions to express the 
+equation in matrix form and solve it using the sparse matrix solver \verb|superlu|. 
+Our implementation, and choice of solver method, resulted in conserved total 
+probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}=1$ for both the single and 
+double slit setup. 
+
+To illustrate the time evolution of the probability function $p_{\ivec, \jvec}^{n} = u_{\ivec, \jvec}^{n*} u_{\ivec, \jvec}^{n}$, 
+we created colormap plots at time steps $t = [0, 0.001, 0.002]$. Since we are working 
+with complex numbers, we included separate plots for each time step of Re$(u_{\ivec, \jvec})$ 
+and Im$(u_{\ivec, \jvec})$.
+% We observed something...
+
+In addition, we determined the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, 
+for single-, double- and triple-slit.
+% We observed something here as well...
+
+\end{document}

latex/sections/introduction.tex

@@ -0,0 +1,65 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\section{Introduction}\label{sec:introduction}
+The nature of light has long been a subject of interest and discussion. %
+% Important part of human behavior is observing and understanding our surroundings. 
+% Many big discoveries have been made through observations, verified by mathematical 
+% explanations. Classical physics is based on calculation predicting something we 
+% verify by observation etc. But what happens when we move down to the microscopic 
+% scale, can we still predict the position of a microscopic ball, also called an atom?
+In classical mechanics, we study the kinematics and dynamics of physical objects, 
+while ignoring their intrinsic properties for simplicity. Newton's second law can be 
+applied to an object to describe its trajectory. % It allows us to describe the 
+% forces acting on an object as well as the motion of the object. We can describe 
+% a planets orbital movement \cite{britannica:2023:kepler}, calculate the ... necessary 
+% to launch satellites into orbit, or simply figure out where a ball is going to land 
+% when you throw it... However, when want to study an object at a microscopic level, 
+% e.g. a single atom, classical mechanics falls short. 
+
+Elementary particles such as electrons, does not abide by the laws of classical mechanics. 
+For several years, scientists did not agree on whether light was a particle or a 
+wave. Through the study of interference of light, and radiation of ideal blackbodies, 
+it has been shown that light has both wavelike and particle-like characteristics. 
+This is known as the wave-particle duality, and was showed by Albert Einstein in 
+1905. %
+% Thomas Young studied the interference of light, and found that light to showed 
+% wavelike characteristics \cite{young:1804:double_slit}. This did not agree with  
+% Newtons particle-theory
+
+Erwin Schrödinger wanted to find a mathematical description of the wave characteristics 
+of matter, supporting the wave-particle idea. He postulated a wave function which varies 
+with position, where the function squared can be interpreted as the probability 
+of finding an electron at a given position. This resulted in the Schrödinger equation, 
+a wave eqution of the energy levels for a hydrogen atom. It also shows how a quantum 
+state evolves with time \cite[p. 81]{wu:2023:quantum}. 
+
+We will simulate the time-dependent Schrödinger equation in two dimensions, to 
+study the light wave interference in the double-slit experiment. In addition, we 
+will include variations of walls such as single- and triple-slit. To solve the equation, 
+we will apply the Crank-Nicolson method in 2+1 dimensions.
+
+% However, according to the Heisenberg uncertainty principle, we can't find dx and/or 
+% dp = 0. dx = sqrt{Var(x)} "spread in position", dp = hat{\Psi}(p) = sqrt{Var(p)} 
+% dx \cdot dp \geq \frac{\hbar}{2}
+% Fourier transform \Psi(x) \doublearrow hat{\Psi}(p)
+% \Psi(x) = \int_{infty}^{infty} (alt sum) hat{\Psi}(p) \cos(px) dp sum of different wave forms
+% Light - particle or wave 
+% - double-slit, blackbodies radiation
+% The nature of light.
+% Position space
+% - classical vs quantum mechanics
+% Intuitions of the behavior of physical objects, predictable. Not like the quantum 
+% when we scale down to the atom, our intuition are not as reliable and prediction 
+% are not as perfect. 
+% Instead of finding the path of a ball, we find all the possible paths a ball can take.
+% The world is not one-dimensional, and modelling it require partial diff eqs
+
+In Section \ref{sec:methods}, we will present the theoretical background for 
+this experiment, as well as the algorithms and tools used in the implementation. 
+Continuing with Section \ref{sec:results}, we will present our results and 
+discuss our findings. Lastly, we will conclude our findings in Section \ref{sec:conclusion}.
+\end{document}
+
+% crank-nicolson method!
+% wave equation

latex/sections/methods.tex

@@ -0,0 +1,280 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\section{Methods}\label{sec:methods} % 
+\subsection{The Schrödinger equation}\label{ssec:schrodinger} %
+% Add something that takes Planck to Schrödinger
+% In classical mechanics, we have Newton laws and conservation of energy. In quantum 
+% mechanics, we have Schrödinger equation.
+The Schrödinger equation has a general form 
+\begin{align}
+    i \hbar \frac{\partial}{\partial t} | \Psi \rangle &= \hat{H} | \Psi \rangle \ ,
+    \label{eq:schrodinger_general}
+\end{align}
+where $i$ is the imaginary number, and $\hbar$ is Plancks constant. $\hat{H}$ is 
+a Hamiltonian operator, which represent the energy for the system, and $| \Psi \rangle$ 
+is the quantum state. In two-dimensional position space, the quantum state can 
+be expressed using the time-dependent complex-valued wave function $\Psi (x, y, t)$. 
+Using Born rule, the square modulus of the wave function is proportional to the 
+probability density of finding a particle at position $(x, y)$ at time t. The relation 
+is given by 
+\begin{align}
+    p(x, y \ | \ t) &= |\Psi(x, y, t)|^{2} = \Psi^{*}(x, y, t) \Psi(x, y, t) \ ,
+    \label{eq:born_rule}
+\end{align}
+where $\Psi^{*}$ denotes the complex conjugated wave function. 
+% Add something about kinetic and potential energy, to introduce the potential V
+When the potential is time-independent, the Schrödinger equation can be expressed as 
+\begin{align*}
+    i \hbar \frac{\partial}{\partial t} \Psi (x, y, t) &= - \frac{\hbar^{2}}{2m} \bigg( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \bigg) \Psi (x, y, t) \\
+    & \quad + V(x, y, t) \Psi (x, y, t) \numberthis \ .
+    \label{eq:schrodinger_special}
+\end{align*}
+The partial derivatives (...) gives the kinetic energy, and the potental $V$ is 
+the external environment. In this experiment we will only consider the case where 
+the potential is time-independent, resulting in $V = V(x, y)$
+When we scale Schrödinger equation by the dimensionful variables, we are left with 
+a wave function $u$, potential $v$ and the dimensionless equation
+\begin{align}
+    i \frac{\partial u}{\partial t} &= - \frac{\partial^{2} u}{\partial x^{2}} - \frac{\partial^{2} u}{\partial y^{2}} + v(x, y) u \ .
+    \label{eq:schrodinger_dimensionless}
+\end{align} % 
+This gives us the Born rule 
+\begin{align}
+    p(x, y \ | \ t) &= |u(x, y, t)|^{2} = u^{*}(x, y, t) u(x, y, t) \ .
+    \label{eq:born_rule_scaled}
+\end{align}
+
+
+
+\subsection{The Crank-Nicolson scheme}\label{ssec:crank_nicolson} %
+When we evaluate a particle in position space, we have to consider partial differential 
+equations (PDE). To solve these numerically, we have to discretize Equation \eqref{eq:schrodinger_dimensionless}. 
+We use the $\theta$-rule \footnote{Using the $\theta$-rule, we can derive Forward Euler using $\theta = 1$, and Backward Euler using $\theta = 0$}, 
+to combine the forward (explicit) and backward (implicit) finite difference method. 
+The result is a linear combination of the explicit and implicit scheme, given by 
+\begin{align}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \theta F_{\ivec, \jvec}^{n+1} + (1 - \theta) F_{\ivec, \jvec}^{n} \ ,
+    \label{eq:theta_rule}
+\end{align} %
+where $\theta \in [0, 1]$.
+To simplify notation and avoid confusion of indices with the imaginary number $i$, 
+we have used the notation $\ivec, \jvec$ in subscript to indicate the commonly named indices $i, j$ 
+in x- and y-direction. In addition, the superscript $n, n+1$ indicate position in time. 
+We get the Crank-Nicolson method (CN) when $\theta = 1/2$ gives Crank-Nicolson
+\begin{align}
+    \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \ .
+    \label{eq:crank_nicolson_method}
+\end{align} %
+Using CN, we derive the discretized Schrödinger equation given by 
+\begin{align*}
+    & u_{\ivec, \jvec}^{n+1} - r \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
+    & - r \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
+    &= u_{\ivec, \jvec}^{n} + r \big[ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n} \big] \\
+    & \quad + r \big[ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n} \big] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} \numberthis \ ,
+    \label{eq:schrodinger_discretized}
+\end{align*} %
+where $r$ is defined as
+\begin{align*}
+    r \equiv \frac{i \Delta t}{2 \Delta h^{2}}
+\end{align*} %
+The full derivation of both Equation \eqref{eq:crank_nicolson_method} and Equation \eqref{eq:schrodinger_discretized} 
+can be found in Appendix \ref{ap:crank_nicolson}.
+
+
+\subsection{The double-slit experiment}\label{ssec:double_slit} % 
+Thomas Young first performed the double-slit experiment in 1801 to demonstrate the  
+principle of interference of light \cite{britannica:2023:young}, while postulating 
+light as waves rather than particles. The double-slit experiment result in a diffraction 
+pattern on a detector screen, where constructive interference of light result in 
+bright spots, and destructive interference result in dark spots as showed in Figure 
+\ref{fig:youngs_double_slit}.
+\begin{figure}
+    \centering 
+    \includegraphics[width=\linewidth]{images/youngs_double_slit.pdf}
+    \caption{The setup of Thomas Young's double slit experiment, where $S_{0}$ denotes 
+    the light source, $S_{1}$ and $S_{2}$ denotes the slits in the wall.}
+    \label{fig:youngs_double_slit}
+\end{figure}
+
+After the wave passes through the two slits, the pattern observed is determined by 
+the path difference determined by 
+\begin{align*}
+    \delta = d \sin (\theta) = m \lambda \ ,
+\end{align*}
+where $\lambda$ is the wavelength and $m$ is called the order number. $d$ is the 
+distance between the center of the two slits, assuming that the distance between 
+the wall and the detector screen $L >> \delta$ \cite[p. 6]{mit:2004:physics}. In 
+this case, we observe constructive interference when 
+\begin{align*}
+    \delta = m \lambda && m = 0, \pm 1, \pm 2 \dots \ ,
+\end{align*} 
+and destructive interference when
+\begin{align*}
+    \delta = (m + \frac{1}{2}) \lambda && m = 0, \pm 1, \pm 2 \dots \ ,
+\end{align*} 
+
+% Something about Heisenberg uncertainty principle 
+
+\subsection{Implementation}\label{ssec:implementation} % 
+% Add tables of parameters used, initial conditions, notation etc.
+A, B are sparse csc matrix
+- theory of csc matrix?
+% \begin{equation*}
+% \begin{pNiceArray}{ccc|ccc}
+%     \bullet & \bullet &  & \bullet &  &  &  &  &  \\
+%     \bullet & \bullet & \bullet &  & \bullet &  &  &  & \\
+%       & \bullet & \bullet &  &  & \bullet &  &  & 
+% \end{pNiceArray}
+% \end{equation*}
+We use Dirichlet boundary conditions, as given in Table \ref{tab:boundary_conditions}, 
+which allows us to express Equation \eqref{eq:schrodinger_discretized} as a matrix 
+equation 
+\begin{align}
+    A u^{n+1} = B u^{n} \ .
+\end{align}
+Here, both $u^{n+1}$ and $u^{n}$ are column vectors containing the internal points 
+of the $xy$ grid at time step $n+1$ and $n$, respectively. Since we have $M$ points 
+in $x$- and $y$-direction, we have $M-2$ internal points. Both $u$ vectors have 
+length $(M-2)^{2}$, and the matrices $A$ and $B$ have size $(M-2)^{2} \times (M-2)^{2}$. 
+The matrices can be decomposed as submatrices of size $(M-2) \times (M-2)$, with 
+the following pattern 
+\begin{align*}
+    A, B = 
+    \begin{bmatrix}
+        \begin{matrix}
+            \bullet & \bullet & \phantom{\bullet} \\
+            \bullet & \bullet & \bullet \\
+            \phantom{\bullet} & \bullet & \bullet 
+        \end{matrix} 
+        & \rvline &
+        \begin{matrix}
+            \bullet & \phantom{\bullet} & \phantom{\bullet} \\
+            \phantom{\bullet} & \bullet & \phantom{\bullet} \\
+            \phantom{\bullet} & \phantom{\bullet} & \bullet 
+        \end{matrix} 
+        & \rvline &
+        \begin{matrix}
+            \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
+            \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
+            \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet}
+        \end{matrix} \\
+        \hline
+        \begin{matrix}
+            \bullet & \phantom{\bullet} & \phantom{\bullet} \\
+            \phantom{\bullet} & \bullet & \phantom{\bullet} \\
+            \phantom{\bullet} & \phantom{\bullet} & \bullet 
+        \end{matrix} 
+        & \rvline &
+        \begin{matrix}
+            \bullet & \bullet & \phantom{\bullet} \\
+            \bullet & \bullet & \bullet \\
+            \phantom{\bullet} & \bullet & \bullet 
+        \end{matrix} 
+        & \rvline &
+        \begin{matrix}
+            \bullet & \phantom{\bullet} & \phantom{\bullet} \\
+            \phantom{\bullet} & \bullet & \phantom{\bullet} \\
+            \phantom{\bullet} & \phantom{\bullet} & \bullet 
+        \end{matrix} \\
+        \hline
+        \begin{matrix}
+           \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
+           \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
+           \phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet}
+       \end{matrix} 
+       & \rvline &
+       \begin{matrix}
+           \bullet & \phantom{\bullet} & \phantom{\bullet} \\
+           \phantom{\bullet} & \bullet & \phantom{\bullet} \\
+           \phantom{\bullet} & \phantom{\bullet} & \bullet 
+       \end{matrix} 
+       & \rvline &
+       \begin{matrix}
+           \bullet & \bullet & \phantom{\bullet} \\
+           \bullet & \bullet & \bullet \\
+           \phantom{\bullet} & \bullet & \bullet 
+       \end{matrix}
+    \end{bmatrix}
+\end{align*}
+To fill the matrices $A$ and $B$, we used 
+\begin{align*}
+    a_{k} &= 1 + 4r + \frac{i \Delta t}{2} v_{\ivec, \jvec} \\
+    b_{k} &= 1 - 4r - \frac{i \Delta t}{2} v_{\ivec, \jvec} \ .
+\end{align*} 
+An example of filled matrices can be found in Appendix \ref{ap:matrix_structure}.
+
+Notations:
+In addition, we use an equal step size in x- and y-direction, $h$ such that 
+\begin{align*}
+    x \in [0, 1] && x \rightarrow x_{\ivec} = \ivec h && \ivec = 0, 1, \dots, M-1 \\
+    y \in [0, 1] && y \rightarrow y_{\jvec} = \jvec h && \jvec = 0, 1, \dots, M-1 \\
+    t \in [0, T] && t \rightarrow t_{n} = n \Delta t && n = 0, 1, \dots, N_{t}-1 
+\end{align*}
+And simplify indices such that 
+\begin{align*}
+    u(x, y, t) \rightarrow u(\ivec h, \jvec h, n \Delta t) \equiv u_{\ivec, \jvec}^{n} \\
+    v(x, y) \rightarrow u(\ivec h, \jvec h) \equiv v_{\ivec, \jvec}
+\end{align*}
+which gives a matrix $U^{n}$ that contains elements $u_{\ivec, \jvec}^{n}$, and 
+a matrix $V$ that contains elements $v_{\ivec, \jvec}$.
+\begin{table}[H]
+    \centering
+    \begin{tabular}{l r} % @{\extracolsep{\fill}}
+        \hline
+        Position & Value  \\
+        \hline 
+        $u(x=0, y, t)$ & $0$  \\
+        $u(x=1, y, t)$ & $0$ \\
+        $u(x, y=0, t)$ & $0$  \\
+        $u(x, y=1, t)$ & $0$  \\
+        \hline
+    \end{tabular}
+    \caption{Boundary conditions in the xy-plane, also known as Dirichlet boundary conditions.}
+    \label{tab:boundary_conditions}
+\end{table}
+
+For the general setup of the wall, we used 
+\begin{table}[H]
+    \centering
+    \begin{tabular}{l r} % @{\extracolsep{\fill}}
+        \hline
+        Parameter & Value  \\
+        \hline 
+        Wall thickness & $0.02$  \\
+        Wall position & $0.5$ \\
+        Separator length & $0.05$  \\
+        Slit aperture & $0.05$  \\
+        \hline
+    \end{tabular}
+    \caption{Wall setup.}
+    \label{tab:wall_setup}
+\end{table}
+
+\begin{table}[H]
+    \centering
+    \begin{tabular}{l r r} % @{\extracolsep{\fill}}
+        \hline
+        Simulation & $1$ & $2$ \\
+        \hline 
+        $h$ & $0.005$ & $0.005$ \\
+        $\Delta t$ & $2.5 \times 10^{-5}$ & $2.5 \times 10^{-5}$ \\
+        $T$ & $0.008$ & $0.002$ \\
+        $x_{c}$ & $0.25$ & $0.25$ \\
+        $\sigma_{x}$ & $0.05$ & $0.05$ \\
+        $p_{x}$ & $200$ & $200$ \\
+        $y_{c}$ & $0.5$ & $0.5$ \\
+        $\sigma_{y}$ & $0.05$ & $0.20$ \\
+        $p_{y}$ & $0$ & $0$ \\
+        $v_{0}$ & $0$ & $1 \times10^{10}$ \\
+        \hline
+    \end{tabular}
+    \caption{Wall setup.}
+    \label{tab:sim_setup}
+\end{table}
+
+\subsection{Tools}\label{ssec:tools} %
+The double-slit experiment is implemented in C++. We use the Python library 
+\verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and 
+\verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures. 
+\end{document}

latex/sections/results.tex

@@ -0,0 +1,35 @@
+\documentclass[../schrodinger_simulation.tex]{subfiles}
+
+\begin{document}
+\section{Results}\label{sec:results}
+\subsection{Deviation}\label{ssec:deviation}
+% Problem 3: Discuss approaches to solve Au^{n+1} = b, dealing with sparse matrix...
+We used the superlu solver, which is a dedicated solver for sparse matrices. It is 
+generally used to solve nonsymmetric, sparse matrices. However, as the lapack solver 
+converts the sparse matrix to a dense matrix, it will increase memory usage compared 
+to superlu.
+% Problem 7: Consequenses of solver choice, in regards to accuracy of probability conserved
+%            Add plot of deviation for both single- and double-slit
+Since we use a solver for sparse matrices, we decrease number of computations performed 
+compared to solver using dense matrix. We check if the total probability is conserved 
+over time, by plotting the deviation $s$ as
+\begin{align*}
+    s^{n} = 1 - \sum_{\ivec , \jvec} p_{\ivec , \jvec}^{n} = 1 - \sum_{\ivec , \jvec} u_{\ivec , \jvec}^{n*} u_{\ivec , \jvec}^{n} \ .
+\end{align*}
+The deviation as a function of time is plotted in Figure \ref{fig:deviation}.
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/probability_deviation.pdf}
+    \caption{Deviation for $t \in [0, T]$ where $T=0.008$.}
+    \label{fig:deviation}
+\end{figure}
+
+\subsection{Time evolution}\label{ssec:time_evolution}
+% Problem 8: Colormap, include plot of both Re and Im for different time steps
+%            Account for color scale
+
+\subsection{Particle detection}\label{ssec:particle_detection}
+% Problem 9: Plot detection probability for single-, double- and triple-slit
+
+
+\end{document}

lib/WaveSimulation.cpp

@@ -17,56 +17,113 @@
 #include <complex>
 #include <cstdint>
 #include <cstdlib>
+#include <vector>
 
-WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
-                               double y_c, double sigma_x, double sigma_y,
-                               double p_x, double p_y, double thickness,
-                               double pos_x, double ap_sep, double ap,
-                               uint32_t slits)
+// Initializers
+void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
+                                  double sigma_y, double p_x, double p_y)
 {
-    this->dt = dt;
-    this->h = h;
-    this->T = T;
-    this->M = 1. / h;
-    this->N = M - 2;
-    this->build_V(thickness, pos_x, ap_sep, ap, slits);
-    this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
-    this->build_A();
-    this->build_B();
+    this->U.set_size(this->N, this->N);
+    double x, y, diff_x, diff_y;
+    std::complex<double> sum = 0.;
+    for (size_t j = 0; j < this->U.n_cols; j++) {
+        x = j * h;
+        diff_x = x - x_c;
+        for (size_t i = 0; i < this->U.n_rows; i++) {
+            y = i * h;
+            diff_y = y - y_c;
+            this->U(i, j) =
+                std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x)
+                         - (diff_y * diff_y) / (2. * sigma_y * sigma_y)
+                         + p_x * x * 1._i + p_y * y * 1._i);
+            sum += this->U(i, j) * std::conj(this->U(i, j));
+        }
+    }
+
+    // The imaginary part of sum should be 0.
+    if (std::abs(sum.imag()) > 1e-7) {
+        abort();
+    }
+
+    double norm = 1. / std::sqrt(sum.real());
+
+    // Normalize each element
+    this->U.for_each([norm](std::complex<double> &el) { el *= norm; });
 }
 
-WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c, double y_c,
-               double sigma_x, double sigma_y, double p_x, double p_y)
+void WaveSimulation::initialize_V()
 {
-    this->dt = dt;
-    this->h = h;
-    this->T = T;
-    this->M = 1. / h;
-    this->N = M - 2;
-    this->build_V(0.,0.,0.,0., 0);
-    this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
-    this->build_A();
-    this->build_B();
-
+    this->V.set_size(this->N, this->N);
+    this->V.fill(0.);
 }
 
-void WaveSimulation::solve(std::ofstream &ofile)
+void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
+                                  double ap, uint32_t slits)
 {
-    ofile << this->N << '\n';
-    uint32_t iterations = this->T / this->dt;
-    for (size_t i = 0; i < iterations; i++) {
-        this->write_U(ofile);
-        this->step();
+    this->V.set_size(this->N, this->N);
+    this->V.fill(0.);
+    if (slits == 0) {
+        return;
+    }
+
+    arma::cx_vec res;
+
+    // Differentiate between odd and even number of slits.
+    uint32_t ap_points, ap_sep_points;
+
+    if (slits % 2 == 0) {
+        ap_points = ap / this->h;
+        ap_sep_points = ap_sep / this->h
+                        + ((uint32_t)(ap_sep / this->h) % 2 != this->N % 2);
+
+        res = arma::cx_vec(ap_sep_points, arma::fill::value(1e10));
+
+        for (size_t i = 0; i < slits; i += 2) {
+            res = arma::join_cols(res,
+                                  arma::cx_vec(ap_points, arma::fill::zeros));
+            res = arma::join_cols(arma::cx_vec(ap_points, arma::fill::zeros),
+                                  res);
+            res = arma::join_cols(
+                res, arma::cx_vec(ap_sep_points, arma::fill::value(1e10)));
+            res = arma::join_cols(
+                arma::cx_vec(ap_sep_points, arma::fill::value(1e10)), res);
+        }
+    }
+    else {
+        ap_points =
+            ap / this->h + ((uint32_t)(ap / this->h) % 2 != this->N % 2);
+        ap_sep_points = ap_sep / this->h;
+
+        res = arma::cx_vec(ap_points, arma::fill::value(0));
+
+        for (size_t i = 0; i < slits - 1; i += 2) {
+            res = arma::join_cols(
+                res, arma::cx_vec(ap_sep_points, arma::fill::value(1e10)));
+            res = arma::join_cols(
+                arma::cx_vec(ap_sep_points, arma::fill::value(1e10)), res);
+            res = arma::join_cols(res,
+                                  arma::cx_vec(ap_points, arma::fill::zeros));
+            res = arma::join_cols(arma::cx_vec(ap_points, arma::fill::zeros),
+                                  res);
+        }
+    }
+
+    if (res.size() > this->N) {
+        abort();
+    }
+
+    uint32_t fill = (this->N - res.size()) / 2;
+    res = arma::join_cols(arma::cx_vec(fill, arma::fill::value(1e10)), res);
+    res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2),
+                                            arma::fill::value(1e10)));
+
+    uint32_t start = pos_x / this->h - thickness / this->h / 2;
+    for (size_t i = 0; i < thickness / this->h; i++) {
+        this->V.col(start + i) = res;
     }
 }
 
-void WaveSimulation::step()
-{
-    arma::cx_vec tmp = this->B * this->U.as_col();
-    arma::spsolve(this->U, this->A, tmp);
-}
-
-void WaveSimulation::build_A()
+void WaveSimulation::initialize_A()
 {
     // Create the diagonal
     arma::cx_vec diagonal(this->N * this->N);
@@ -96,7 +153,7 @@ void WaveSimulation::build_A()
     this->A.diag(this->N).fill(-r);
 }
 
-void WaveSimulation::build_B()
+void WaveSimulation::initialize_B()
 {
     std::complex<double> r = (1._i * this->dt) / (2 * h * h);
 
@@ -126,92 +183,149 @@ void WaveSimulation::build_B()
     this->B.diag(this->N).fill(r);
 }
 
-void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
-                                  double sigma_y, double p_x, double p_y)
+// Constructors
+WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
+                               double y_c, double sigma_x, double sigma_y,
+                               double p_x, double p_y, double thickness,
+                               double pos_x, double ap_sep, double ap,
+                               uint32_t slits)
 {
-    this->U.set_size(this->N, this->N);
-    double x, y, diff_x, diff_y;
-    std::complex<double> sum = 0.;
-    for (size_t j = 0; j < this->U.n_cols; j++) {
-        x = j * h;
-        diff_x = x - x_c;
-        for (size_t i = 0; i < this->U.n_rows; i++) {
-            y = i * h;
-            diff_y = y - y_c;
-            this->U(i, j) =
-                std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x)
-                         - (diff_y * diff_y) / (2. * sigma_y * sigma_y)
-                         + p_x * x * 1._i + p_y * y * 1._i);
-            sum += this->U(i, j) * std::conj(this->U(i, j));
-        }
-    }
-    if (std::abs(sum.imag()) > 1e-7) {
+    this->dt = dt;
+    this->h = h;
+    this->T = T;
+    this->M = 1. / h;
+    this->N = M - 2;
+    this->initialize_V(thickness, pos_x, ap_sep, ap, slits);
+    this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
+    this->initialize_A();
+    this->initialize_B();
+}
+
+WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
+                               double y_c, double sigma_x, double sigma_y,
+                               double p_x, double p_y)
+{
+    this->dt = dt;
+    this->h = h;
+    this->T = T;
+    this->M = 1. / h;
+    this->N = M - 2;
+    this->initialize_V();
+    this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
+    this->initialize_A();
+    this->initialize_B();
+}
+
+// Public methods
+void WaveSimulation::step()
+{
+    // Create the right hand side of Ax = b
+    arma::cx_vec tmp = this->B * this->U.as_col();
+
+    // Solve Ax = b
+    arma::spsolve(this->U, this->A, tmp);
+}
+
+void WaveSimulation::solve(std::string outfile, bool write_each_step)
+{
+    // Create path and proceed if successful.
+    if (!utils::mkpath(utils::dirname(outfile))) {
         abort();
     }
-    double norm = 1. / std::sqrt(sum.real());
 
-    this->U.for_each([norm](std::complex<double> &el) { el *= norm; });
+    // Open file
+    std::ofstream ofile;
+    ofile.open(outfile);
+
+    // Write the size of the matrix on the first line
+    ofile << this->N << '\n';
+
+    uint32_t iterations = this->T / this->dt;
+
+    if (write_each_step) {
+        for (size_t i = 0; i < iterations; i++) {
+            this->write_U(ofile);
+            this->step();
+        }
+    }
+    else {
+        for (size_t i = 0; i < iterations; i++) {
+            this->step();
+        }
+    }
+    this->write_U(ofile);
+    ofile.close();
+}
+void WaveSimulation::solve(std::string outfile, std::vector<double> &steps)
+{
+    // Create path and proceed if successful.
+    if (!utils::mkpath(utils::dirname(outfile))) {
+        abort();
+    }
+
+    // Open file
+    std::ofstream ofile;
+    ofile.open(outfile);
+
+    // Write the size of the matrix on the first line
+    ofile << this->N << '\n';
+    
+    uint32_t iterations;
+
+    for (size_t i=0; i < steps.size(); i++) {
+        if (i == 0) {
+            iterations = steps[i] / this->dt;
+        }
+        else {
+            iterations = (steps[i] - steps[i-1]) / this->dt;
+        }
+
+        for (size_t j=0; j < iterations; j++) {
+            this->step();
+        }
+        this->write_U(ofile);
+    }
+
+    ofile.close();
+}
+
+void WaveSimulation::probability_deviation(std::string outfile,
+                                           bool write_each_step)
+{
+    // Create path and proceed if successful.
+    if (!utils::mkpath(utils::dirname(outfile))) {
+        abort();
+    }
+
+    // Open file
+    std::ofstream ofile;
+    ofile.open(outfile);
+
+    uint32_t iterations = this->T / this->dt;
+
+    double sum;
+
+    if (write_each_step) {
+        for (size_t i = 0; i < iterations; i++) {
+            sum = arma::accu(this->U % arma::conj(this->U)).real();
+            ofile << i*this->dt << '\t' << 1. - sum << '\n';
+            this->step();
+        }
+    }
+    else {
+        for (size_t i = 0; i < iterations; i++) {
+            this->step();
+        }
+    }
+    sum = arma::accu(this->U % arma::conj(this->U)).real();
+    ofile << this->T << '\t' << 1. - sum << '\n';
+    ofile.close();
 }
 
 void WaveSimulation::write_U(std::ofstream &ofile)
 {
+    // Write each element to file in column-major order.
     this->U.for_each(
         [&ofile](std::complex<double> el) { ofile << el << '\t'; });
     ofile << '\n';
 }
-
-void WaveSimulation::build_V(double thickness, double pos_x,
-                             double aperture_separation, double aperture,
-                             uint32_t slits)
-{
-    this->V.set_size(this->N, this->N);
-    this->V.fill(0.);
-    if (slits == 0) {
-        return;
-    }
-    arma::cx_vec res;
-
-    if (slits % 2 == 0) {
-        res = arma::cx_vec(aperture_separation / this->h,
-                           arma::fill::value(1e10));
-        for (size_t i = 0; i < slits; i += 2) {
-            res = arma::join_cols(
-                res, arma::cx_vec(aperture / this->h, arma::fill::zeros));
-            res = arma::join_cols(
-                arma::cx_vec(aperture / this->h, arma::fill::zeros), res);
-            res =
-                arma::join_cols(res, arma::cx_vec(aperture_separation / this->h,
-                                                  arma::fill::value(1e10)));
-            res = arma::join_cols(arma::cx_vec(aperture_separation / this->h,
-                                               arma::fill::value(1e10)),
-                                  res);
-        }
-    }
-    else {
-        res = arma::cx_vec(aperture / this->h, arma::fill::value(0));
-        for (size_t i = 0; i < slits - 1; i += 2) {
-            res =
-                arma::join_cols(res, arma::cx_vec(aperture_separation / this->h,
-                                                  arma::fill::value(1e10)));
-            res = arma::join_cols(arma::cx_vec(aperture_separation / this->h,
-                                               arma::fill::value(1e10)),
-                                  res);
-            res = arma::join_cols(
-                res, arma::cx_vec(aperture / this->h, arma::fill::zeros));
-            res = arma::join_cols(
-                arma::cx_vec(aperture / this->h, arma::fill::zeros), res);
-        }
-    }
-    if (res.size() > this->N) {
-        abort();
-    }
-    uint32_t fill = (this->N - res.size()) / 2;
-    res = arma::join_cols(arma::cx_vec(fill, arma::fill::value(1e10)), res);
-    res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2),
-                                            arma::fill::value(1e10)));
-
-    uint32_t start = pos_x / this->h - thickness / this->h / 2;
-    for (size_t i = 0; i < thickness / this->h; i++) {
-        this->V.col(start + i) = res;
-    }
-}

lib/utils.cpp

@@ -45,14 +45,14 @@ bool mkpath(std::string path, int mode)
             cur_dir = path.substr(0, pos);
             if (mkdir(cur_dir.c_str(), mode) != 0
                 && stat(cur_dir.c_str(), &buf) != 0) {
-                return -1;
+                return false;
             }
         }
         else {
             break;
         }
     }
-    return 0;
+    return true;
 }
 
 std::string dirname(const std::string &path)

python_scripts/colormap.py

@@ -1,26 +1,76 @@
-import numpy as np
-import matplotlib.pyplot as plt
 import ast
+
+import matplotlib.pyplot as plt
+import numpy as np
 import seaborn as sns
 
 sns.set_theme()
+params = {
+    "font.family": "Serif",
+    "font.serif": "Roman",
+    "text.usetex": True,
+    "axes.titlesize": "large",
+    "axes.labelsize": "large",
+    "xtick.labelsize": "large",
+    "ytick.labelsize": "large",
+    "legend.fontsize": "medium",
+}
+plt.rcParams.update(params)
+
 
 def plot():
     with open("data/color_map.txt") as f:
-        lines = f.readlines();
+        lines = f.readlines()
         size = int(lines[0])
         for i, line in enumerate(lines[1:]):
-            fig, ax = plt.subplots()
+            # Create figures for each plot
+            fig1, ax1 = plt.subplots()
+            fig2, ax2 = plt.subplots()
+            fig3, ax3 = plt.subplots()
+
             arr = line.strip().split("\t")
-            arr = np.asarray(list(map(lambda x: ((a := complex(*ast.literal_eval(x)))*a.conjugate()).real, arr)))
+            arr = np.asarray(list(map(lambda x: complex(*ast.literal_eval(x)), arr)))
+            # Reshape and transpose array
+            arr = arr.reshape(size, size).T
 
-            arr = arr.reshape(size,size)
+            # Plot color maps
+            color_map1 = ax1.imshow(
+                np.multiply(arr, arr.conj()).real,
+                interpolation="nearest",
+                cmap=sns.color_palette("mako", as_cmap=True),
+            )
+            color_map2 = ax2.imshow(
+                arr.real,
+                interpolation="nearest",
+                cmap=sns.color_palette("mako", as_cmap=True),
+            )
+            color_map3 = ax3.imshow(
+                arr.imag,
+                interpolation="nearest",
+                cmap=sns.color_palette("mako", as_cmap=True),
+            )
 
-            ax.imshow(arr.T, interpolation="nearest")
-            plt.savefig(f"latex/images/color_map_{i}.pdf")
-            plt.close(fig)
+            # Create color bar
+            fig1.colorbar(color_map1, ax=ax1)
+            fig2.colorbar(color_map2, ax=ax2)
+            fig3.colorbar(color_map3, ax=ax3)
+
+            # Remove grids
+            ax1.grid(False)
+            ax2.grid(False)
+            ax3.grid(False)
+
+            # Save the figures
+            fig1.savefig(f"latex/images/color_map_{i}_prob.pdf")
+            fig2.savefig(f"latex/images/color_map_{i}_real.pdf")
+            fig3.savefig(f"latex/images/color_map_{i}_imag.pdf")
+
+            # Close figures
+            plt.close(fig1)
+            plt.close(fig2)
+            plt.close(fig3)
 
 
 if __name__ == "__main__":
     plot()
-    
+

python_scripts/detector.py

@@ -4,6 +4,17 @@ import ast
 import seaborn as sns
 
 sns.set_theme()
+params = {
+    "font.family": "Serif",
+    "font.serif": "Roman",
+    "text.usetex": True,
+    "axes.titlesize": "large",
+    "axes.labelsize": "large",
+    "xtick.labelsize": "large",
+    "ytick.labelsize": "large",
+    "legend.fontsize": "medium",
+}
+plt.rcParams.update(params)
 
 def plot():
     files = [

python_scripts/heat_map.py

@@ -4,6 +4,19 @@ import matplotlib
 from matplotlib.animation import FuncAnimation
 import ast
 
+import seaborn as sns 
+params = {
+    "font.family": "Serif",
+    "font.serif": "Roman",
+    "text.usetex": True,
+    "axes.titlesize": "large",
+    "axes.labelsize": "large",
+    "xtick.labelsize": "large",
+    "ytick.labelsize": "large",
+    "legend.fontsize": "medium",
+}
+plt.rcParams.update(params)
+
 wave_arr = []
 fig = plt.figure()
 ax = plt.gca()

python_scripts/normalization.py

@@ -1,23 +1,38 @@
 import numpy as np
 import matplotlib.pyplot as plt
-import matplotlib
-from matplotlib.animation import FuncAnimation
-import ast
+
+import seaborn as sns 
+
+sns.set_theme()
+params = {
+    "font.family": "Serif",
+    "font.serif": "Roman",
+    "text.usetex": True,
+    "axes.titlesize": "large",
+    "axes.labelsize": "large",
+    "xtick.labelsize": "large",
+    "ytick.labelsize": "large",
+    "legend.fontsize": "medium",
+}
+plt.rcParams.update(params)
 
 def plot():
-    with open("data/probability_deviation.txt") as f:
-        lines = f.readlines();
-        x = []
-        arr_narrow = []
-        arr_wide = []
-        for line in lines:
-            tmp = line.strip().split(",")
-            x.append(float(tmp[0]))
-            arr_narrow.append(float(tmp[1]))
-            arr_wide.append(float(tmp[2]))
+    files = [
+        "data/probability_deviation_no_slits.txt",
+        "data/probability_deviation_slits.txt",
+    ]
+    for file in files:
+        with open(file) as f:
+            lines = f.readlines();
+            x = []
+            arr = []
+            for line in lines:
+                tmp = line.strip().split("\t")
+                x.append(float(tmp[0]))
+                arr.append(float(tmp[1]))
+
+            plt.plot(x,arr)
 
-        plt.plot(x,arr_narrow)
-        plt.plot(x,arr_wide)
         plt.savefig("latex/images/probability_deviation.pdf")
 
 if __name__ == "__main__":

python_scripts/plot_v.py

@@ -4,6 +4,19 @@ import matplotlib
 from matplotlib.animation import FuncAnimation
 import ast
 
+import seaborn as sns 
+params = {
+    "font.family": "Serif",
+    "font.serif": "Roman",
+    "text.usetex": True,
+    "axes.titlesize": "large",
+    "axes.labelsize": "large",
+    "xtick.labelsize": "large",
+    "ytick.labelsize": "large",
+    "legend.fontsize": "medium",
+}
+plt.rcParams.update(params)
+
 
 def plot():
     with open("v.txt") as f:

src/main.cpp

@@ -11,36 +11,20 @@
  * */
 #include "WaveSimulation.hpp"
 #include "utils.hpp"
+
 #include <fstream>
 
 void probability_deviation()
 {
-    WaveSimulation sim_narrow(.005, 2.5e-5, .008, .25, .5, .05, .05, 200., 0.);
+    WaveSimulation sim_no_slits(.005, 2.5e-5, .008, .25, .5, .05, .05, 200.,
+                                0.);
 
-    WaveSimulation sim_wide(.005, 2.5e-5, .008, .25, .5, .05, .10, 200., 0.,
-                            0.02, .5, .05, .05, 2);
+    WaveSimulation sim_slits(.005, 2.5e-5, .008, .25, .5, .05, .10, 200., 0.,
+                             0.02, .5, .05, .05, 2);
 
-    std::ofstream ofile;
-    utils::mkpath("data");
-    ofile.open("data/probability_deviation.txt");
-
-    double sum_narrow, sum_wide;
-
-    for (size_t i = 0; i < 320; i++) {
-        sum_narrow = 0;
-        sum_wide = 0;
-        for (size_t j = 0; j < sim_narrow.U.n_elem; j++) {
-            sum_narrow += (sim_narrow.U(j) * std::conj(sim_narrow.U(j))).real();
-            sum_wide += (sim_wide.U(j) * std::conj(sim_wide.U(j))).real();
-        }
-
-        sim_narrow.step();
-        sim_wide.step();
-        ofile << i << ',' << utils::scientific_format(sum_narrow) << ','
-              << utils::scientific_format(sum_wide) << '\n';
-    }
-
-    ofile.close();
+    sim_no_slits.probability_deviation(
+        "data/probability_deviation_no_slits.txt", true);
+    sim_slits.probability_deviation("data/probability_deviation_slits.txt", true);
 }
 
 void color_map()
@@ -48,19 +32,8 @@ void color_map()
     WaveSimulation sim(.005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0., 0.02,
                        .5, .05, .05, 2);
 
-    std::ofstream ofile;
-    ofile.open("data/color_map.txt");
-    ofile << sim.N << '\n';
-    sim.write_U(ofile);
-    for (size_t i = 0; i < 40; i++) {
-        sim.step();
-    }
-    sim.write_U(ofile);
-    for (size_t i = 0; i < 40; i++) {
-        sim.step();
-    }
-    sim.write_U(ofile);
-    ofile.close();
+    std::vector<double> times{0., .001, .002};
+    sim.solve("data/color_map.txt", times);
 }
 
 void detector_screen()
@@ -68,53 +41,35 @@ void detector_screen()
     WaveSimulation *sim = new WaveSimulation(
         .005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0., 0.02, .5, .05, .05, 1);
 
-    std::ofstream ofile;
-    utils::mkpath("data/screen");
-    ofile.open("data/screen/single_slit.txt");
-    ofile << sim->N << '\n';
-    for (size_t i = 0; i < 80; i++) {
-        sim->step();
-    }
-    sim->write_U(ofile);
-    ofile.close();
+    sim->solve("data/screen/single_slit.txt");
 
     delete sim;
     sim = new WaveSimulation(.005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0.,
                              0.02, .5, .05, .05, 2);
 
-    ofile.open("data/screen/double_slit.txt");
-    ofile << sim->N << '\n';
-    for (size_t i = 0; i < 80; i++) {
-        sim->step();
-    }
-    sim->write_U(ofile);
-    ofile.close();
+    sim->solve("data/screen/double_slit.txt");
 
     delete sim;
     sim = new WaveSimulation(.005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0.,
                              0.02, .5, .05, .05, 3);
 
-    ofile.open("data/screen/triple_slit.txt");
-    ofile << sim->N << '\n';
-    for (size_t i = 0; i < 80; i++) {
-        sim->step();
-    }
-    sim->write_U(ofile);
-    ofile.close();
+    sim->solve("data/screen/triple_slit.txt");
 }
 
 int main()
 {
-    //probability_deviation();
-    //color_map();
+    DEBUG("Before probability deviation");
+    probability_deviation();
+    DEBUG("Before color map");
+    color_map();
+    DEBUG("Before detector screen");
     detector_screen();
-    
-    //std::ofstream ofile;
-    //ofile.open("test.txt");
-    //WaveSimulation sim(.005, 2.5e-5, .008, .25, .5, .05, .10, 200., 0.,
-                            //0.02, .5, .05, .05, 2);
-    //sim.solve(ofile);
 
+    // std::ofstream ofile;
+    // ofile.open("test.txt");
+    // WaveSimulation sim(.005, 2.5e-5, .008, .25, .5, .05, .10, 200., 0.,
+    // 0.02, .5, .05, .05, 2);
+    // sim.solve(ofile);
 
     return 0;
 }