Merge pull request 'janitaws/latex' (#14) from janitaws/latex into develop

Reviewed-on: #14
This commit is contained in:
Cory Balaton 2024-01-01 16:38:44 +00:00
commit bdecd6ad0c
47 changed files with 3005 additions and 28 deletions

View File

@ -1,2 +1,2 @@
CompileFlags: CompileFlags:
Add: [-I../include] Add: [-I../include, -std=c++11]

View File

@ -1,5 +1,5 @@
# The compiler # The compiler
CC=mpicxx CC=g++
# Flags # Flags
CFLAGS=-Wall -larmadillo -std=c++11 -O3 -fopenmp CFLAGS=-Wall -larmadillo -std=c++11 -O3 -fopenmp
@ -26,19 +26,26 @@ LATEXDIR=./latex
# Source directories # Source directories
SRC=./src SRC=./src
LIB=./lib
INCLUDE=./include INCLUDE=./include
# Source files and object file locations # Source files and object file locations
SRCFILES=utils.cpp testlib.cpp data_type.cpp IsingModel.cpp monte_carlo.cpp #SRCFILES=utils.cpp testlib.cpp
SRCS=$(addprefix $(SRC)/, $(SRCS)) #SRCS=$(addprefix $(SRC)/, $(SRCS))
BINOBJS=$(addprefix $(BINOBJDIR)/, $(SRCFILES:.cpp=.o)) #BINOBJS=$(addprefix $(BINOBJDIR)/, $(SRCFILES:.cpp=.o))
PROFOBJS=$(addprefix $(PROFOBJDIR)/, $(SRCFILES:.cpp=.o)) #PROFOBJS=$(addprefix $(PROFOBJDIR)/, $(SRCFILES:.cpp=.o))
DEBUGOBJS=$(addprefix $(DEBUGOBJDIR)/, $(SRCFILES:.cpp=.o)) #DEBUGOBJS=$(addprefix $(DEBUGOBJDIR)/, $(SRCFILES:.cpp=.o))
# Lib files
LIBSRCS=$(notdir $(shell find $(LIB) -type f))
LIBBINOBJS=$(addprefix $(BINOBJDIR)/, $(LIBSRCS:.cpp=.o))
LIBPROFOBJS=$(addprefix $(PROFOBJDIR)/, $(LIBSRCS:.cpp=.o))
LIBDEBUGOBJS=$(addprefix $(DEBUGOBJDIR)/, $(LIBSRCS:.cpp=.o))
# Location for Binaries # Location for Binaries
EXEC=main test_suite phase_transition phase_transition_mpi time pd_estimate mcmc_progression EXEC=$(basename $(notdir $(shell find $(SRC) -type f)))
BINS=$(addprefix $(BINDIR)/, $(EXEC)) BINS=$(addprefix $(BINDIR)/, $(EXEC))
PROFBINS=$(PROFDIR)/phase_transition_mpi PROFBINS=$(addprefix $(PROFDIR)/, $(EXEC))
DEBUGBINS=$(addprefix $(DEBUGDIR)/, $(EXEC)) DEBUGBINS=$(addprefix $(DEBUGDIR)/, $(EXEC))
# List phony targets # List phony targets
@ -56,23 +63,24 @@ latex:
$(MAKE) -C $(LATEXDIR) $(MAKE) -C $(LATEXDIR)
# Rule for binaries # Rule for binaries
$(BINDIR)/%: $(BINOBJDIR)/%.o $(BINOBJS) $(BINDIR)/%: $(BINOBJDIR)/%.o $(LIBBINOBJS)
$(MKDIR) $(BINDIR) $(MKDIR) $(BINDIR)
$(CC) $^ -o $@ $(CFLAGS) -I$(INCLUDE) $(CC) $^ -o $@ $(CFLAGS) -I$(INCLUDE)
# Rule for profiling binaries # Rule for profiling binaries
$(PROFDIR)/%: $(PROFOBJDIR)/%.o $(PROFOBJS) $(PROFDIR)/%: $(PROFOBJDIR)/%.o $(LIBPROFOBJS)
$(MKDIR) $(PROFDIR) $(MKDIR) $(PROFDIR)
$(INSTRUMENT) $(CC) $^ -o $@ $(CFLAGS) $(PROFFLAGS) -I$(INCLUDE) $(INSTRUMENT) $(CC) $^ -o $@ $(CFLAGS) $(PROFFLAGS) -I$(INCLUDE)
# Rule for debug binaries # Rule for debug binaries
$(DEBUGDIR)/%: $(DEBUGOBJDIR)/%.o $(DEBUGOBJS) $(DEBUGDIR)/%: $(DEBUGOBJDIR)/%.o $(LIBDEBUGOBJS)
$(MKDIR) $(DEBUGDIR) $(MKDIR) $(DEBUGDIR)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAGS) -I$(INCLUDE) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAGS) -I$(INCLUDE)
# Rule for object files # Rule for object files
$(BINOBJDIR)/%.o: $(SRC)/%.cpp $(BINOBJDIR)/%.o: $(SRC)/%.cpp
$(MKDIR) $(BINOBJDIR) $(MKDIR) $(BINOBJDIR)
echo $(LIBBINOBJS)
$(CC) -c $^ -o $@ $(CFLAGS) -I$(INCLUDE) $(CC) -c $^ -o $@ $(CFLAGS) -I$(INCLUDE)
# Rule for instrumented object files # Rule for instrumented object files
@ -85,6 +93,21 @@ $(DEBUGOBJDIR)/%.o: $(SRC)/%.cpp
$(MKDIR) $(DEBUGOBJDIR) $(MKDIR) $(DEBUGOBJDIR)
$(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAGS) -I$(INCLUDE) $(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAGS) -I$(INCLUDE)
# Rule for object files
$(BINOBJDIR)/%.o: $(LIB)/%.cpp
$(MKDIR) $(BINOBJDIR)
$(CC) -c $^ -o $@ $(CFLAGS) -I$(INCLUDE)
# Rule for instrumented object files
$(PROFOBJDIR)/%.o: $(LIB)/%.cpp
$(MKDIR) $(PROFOBJDIR)
$(INSTRUMENT) $(CC) -c $^ -o $@ $(CFLAGS) $(PROFFLAGS) -I$(INCLUDE)
# Rule for debug object files
$(DEBUGOBJDIR)/%.o: $(LIB)/%.cpp
$(MKDIR) $(DEBUGOBJDIR)
$(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAGS) -I$(INCLUDE)
# Cleaning # Cleaning
clean: objclean binclean latexclean clean: objclean binclean latexclean

4
data/color_map.txt Normal file

File diff suppressed because one or more lines are too long

View File

@ -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

View File

@ -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

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

116
include/WaveSimulation.hpp Normal file
View File

@ -0,0 +1,116 @@
/** @file WaveSimulation.hpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 0.1
*
* @brief The definition of the WaveSimulation class
*
* @bug No known bugs
* */
#ifndef __WAVE_SIMULATION__
#define __WAVE_SIMULATION__
#include "constants.hpp"
#include "literals.hpp"
#include <armadillo>
#include <cstdint>
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;
/* @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);
/* @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:
/* @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);
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

17
include/constants.hpp Normal file
View File

@ -0,0 +1,17 @@
/** @file constants.hpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief Library of constants
*
* @bug No known bugs
* */
#ifndef __CONST__
#define __CONST__
#define I std::complex<double>{0., 1.}
#endif

19
include/literals.hpp Normal file
View File

@ -0,0 +1,19 @@
/** @file literals.hpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief Useful literals
*
* @bug No known bugs
* */
#ifndef __LITERALS__
#define __LITERALS__
#include <complex>
std::complex<double> operator ""_i(long double magnitude);
#endif

View File

@ -117,10 +117,10 @@ std::string dirname(const std::string &path);
/** @brief Take 2 strings and concatenate them and make sure there is a /** @brief Take 2 strings and concatenate them and make sure there is a
* directory separator (/) between them. * directory separator (/) between them.
* *
* @details This function doesn't care whether or not the values given as * @details This function doesn't care whether or not the values given as
* parameters are valid path strings. It is the responsibility of the user to make * parameters are valid path strings. It is the responsibility of the user to
* sure that the values given are valid path strings. * make sure that the values given are valid path strings. The function only
* The function only guarantees that the output string is a valid path string. * guarantees that the output string is a valid path string.
* *
* @param left The left hand side of the result string * @param left The left hand side of the result string
* @param right The right hand side of the result string * @param right The right hand side of the result string
@ -129,6 +129,8 @@ std::string dirname(const std::string &path);
* */ * */
std::string concatpath(const std::string &left, const std::string &right); std::string concatpath(const std::string &left, const std::string &right);
// A function that prints the structure of a sparse matrix to screen.
void print_sp_matrix_structure(const arma::sp_cx_mat &A);
} // namespace utils } // namespace utils
#endif #endif

440
latex/draft.tex Normal file
View File

@ -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}

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -1,6 +1,141 @@
# 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 # 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 = {Youngs 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 = {Keplers 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 = {Plancks 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 = {Newtons 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 # Math and statistics
@book{lindstrom:2016:kalkulus, @book{lindstrom:2016:kalkulus,
@ -52,6 +187,29 @@
pages = {3021} pages = {3021}
} }
@article{harris:2020:numpy,
title = {Array programming with {NumPy}},
author = {Charles R. Harris and K. Jarrod Millman and St{\'{e}}fan J.
van der Walt and Ralf Gommers and Pauli Virtanen and David
Cournapeau and Eric Wieser and Julian Taylor and Sebastian
Berg and Nathaniel J. Smith and Robert Kern and Matti Picus
and Stephan Hoyer and Marten H. van Kerkwijk and Matthew
Brett and Allan Haldane and Jaime Fern{\'{a}}ndez del
R{\'{i}}o and Mark Wiebe and Pearu Peterson and Pierre
G{\'{e}}rard-Marchant and Kevin Sheppard and Tyler Reddy and
Warren Weckesser and Hameer Abbasi and Christoph Gohlke and
Travis E. Oliphant},
year = {2020},
month = sep,
journal = {Nature},
volume = {585},
number = {7825},
pages = {357--362},
doi = {10.1038/s41586-020-2649-2},
publisher = {Springer Science and Business Media {LLC}},
url = {https://doi.org/10.1038/s41586-020-2649-2}
}
# C++ libraries # C++ libraries
@misc{openmp:2018, @misc{openmp:2018,
author = {OpenMP}, author = {OpenMP},

View File

@ -20,6 +20,7 @@
%% I recommend downloading TeXMaker, because it includes a large library of the most common packages. %% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
\usepackage{physics,amssymb} % mathematical symbols (physics imports amsmath) \usepackage{physics,amssymb} % mathematical symbols (physics imports amsmath)
\usepackage[no-test-for-array]{nicematrix}
\usepackage{graphicx} % include graphics such as plots \usepackage{graphicx} % include graphics such as plots
\graphicspath{.images/} \graphicspath{.images/}
@ -33,6 +34,29 @@
\usepackage{algorithm} \usepackage{algorithm}
\usepackage[noend]{algpseudocode} \usepackage[noend]{algpseudocode}
\usepackage{tikz} \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} % \usetikzlibrary{quantikz}
% defines the color of hyperref objects % defines the color of hyperref objects
% Blending two colors: blue!80!black = 80% blue and 20% black % Blending two colors: blue!80!black = 80% blue and 20% black
@ -44,7 +68,13 @@
% Biblio stuff % Biblio stuff
% \def\biblio{\bibliographystyle{plain}\bibliography{../references/references}} % \def\biblio{\bibliographystyle{plain}\bibliography{../references/references}}
\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}} \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{xr}
\usepackage{subfiles} \usepackage{subfiles}
@ -52,33 +82,36 @@
\begin{document} \begin{document}
\title{Simulating the Schrödinger wave equation using the Crank-Nicolson method in 2+1 dimensions} % self-explanatory \title{Simulating the Schrödinger wave equation using the Crank-Nicolson scheme 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 \author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen \\ \faGithub \, \url{https://github.uio.no/FYS3150-G2-2023/Project-5}} % self-explanatory
\date{\today} % self-explanatory \date{\today} % self-explanatory
\noaffiliation % ignore this, but keep it. \noaffiliation % ignore this, but keep it.
% Abstract % Abstract
% \subfile{sections/abstract} \subfile{sections/abstract}
\maketitle \maketitle
% Introduction % Introduction
% \subfile{sections/introduction} \subfile{sections/introduction}
% Methods % Methods
% \subfile{sections/methods} \subfile{sections/methods}
% Results % Results
% \subfile{sections/results} \subfile{sections/results}
% Conclusion % Conclusion
% \subfile{sections/conclusion} \subfile{sections/conclusion}
% Notes
% \subfile{draft}
\clearpage \clearpage
\newpage \newpage
% Appendix % Appendix
% \subfile{sections/appendices} \subfile{sections/appendices}
\clearpage \clearpage
\onecolumngrid \onecolumngrid
@ -89,4 +122,14 @@
\end{document} \end{document}
% Abstract: OK!
% Introduction: Hyugen, Heisenberg
% Methods: OK!
% Results: OK
% Conclusion: OK!
% Appendices:
% Results
% P7: OK!
% P8: OK!
% P9: OK

View File

@ -2,6 +2,19 @@
\begin{document} \begin{document}
\begin{abstract} \begin{abstract}
We have simulated the two-dimensional time-dependent Schrödinger equation, to study
variations of the double-slit experiment. To derive a discretized equation
we applied the Crank-Nicolson scheme in 2+1 dimensions. 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 a deviation from conserved total probability on the
scale $10^{-14}$, for both the single and double slit setup. To illustrate the time
evolution of the probability function, we created colormap plots for 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{abstract}
\end{document} \end{document}
% $| \sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n} - 1 | \approx $

View File

@ -2,6 +2,159 @@
\begin{document} \begin{document}
\appendix \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.
The Schrödinger equation contains $i$ on the left hand side, we rewrite 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] \ ,
\end{align*}
where $1/i = -i$. Using Equation \eqref{eq:schrodinger_dimensionless}, which is
found in Section \ref{ssec:schrodinger}, 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*}
Since we 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 rewritten 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 a $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*}
\section{Figures}\label{ap:figures}
We created colormap plots of the real and imaginary part of $u_{\ivec, \jvec}$,
in Figure \ref{fig:colormap_real_imag}.
\begin{figure*}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_0_real.pdf}
\caption{Re$(u_{\ivec, \jvec})$ at time $t=0$.}
\label{fig:colormap_0_real}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_1_real.pdf}
\caption{Re$(u_{\ivec, \jvec})$ at time $t=0.001$.}
\label{fig:colormap_1_real}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_2_real.pdf}
\caption{Re$(u_{\ivec, \jvec})$ at time $t=0.002$.}
\label{fig:colormap_2_real}
\end{subfigure}
\newline
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_0_imag.pdf}
\caption{Im$(u_{\ivec, \jvec})$ at time $t=0$.}
\label{fig:colormap_0_imag}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_1_imag.pdf}
\caption{Im$(u_{\ivec, \jvec})$ at time $t=0.001$.}
\label{fig:colormap_1_imag}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{images/color_map_2_imag.pdf}
\caption{Im$(u_{\ivec, \jvec})$ at time $t=0.002$.}
\label{fig:colormap_2_imag}
\end{subfigure}
\caption{The time evolution of Re($u_{\ivec, \jvec}$) and Im($u_{\ivec, \jvec}$),
for time steps $t=[0, 0.001, 0.002]$.}
\label{fig:colormap_real_imag}
\end{figure*}
% \begin{figure*}
% \centering
% \includegraphics[width=0.9\textwidth]{images/color_map_all.pdf}
% \caption{The time evolution of the probability function $p_{\ivec, \jvec}^{n}$ (top row),
% Re($u_{\ivec, \jvec}^{n}$) (middle row), and Im($u_{\ivec, \jvec}^{n}$) (bottom row).
% Time step $t=0$ in the left column, $t=0.001$ in the middle column, and $t=0.002$
% in the right column.}
% \label{fig:colormap_all}
% \end{figure*}
\end{document} \end{document}

View File

@ -2,6 +2,23 @@
\begin{document} \begin{document}
\section{Conclusion}\label{sec:conclusion} \section{Conclusion}\label{sec:conclusion}
We simulated the two-dimensional time-dependent Schrödinger equation, and studied
variations of the double-slit experiment. To solve the partial differential equations
we applied the Crank-Nicolson scheme in 2+1 dimensions, and derived a discretized
equation. We used Dirichlet boundary conditions to simplify the equation,
expressed the equation in matrix form and solved it using the sparse matrix solver
\verb|superlu|. The total probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}$
deviated from $1.0$ by a factor of $10^{-14}$ for both the single and double slit
setup. % Add something about computational accuracy?
We illustrated the time evolution of the probability function $p_{\ivec, \jvec}^{n} = u_{\ivec, \jvec}^{n*} u_{\ivec, \jvec}^{n}$,
using colormap plots for time steps $t = \{0, 0.001, 0.002\}$. In addition, we included
separate plots for each time step of Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$,
to show the components of the complex values. This resulted in visible diffraction
patterns for the double-slit experiment.
In addition, we studied the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$,
for single-, double- and triple-slit setups. We found that increasing the number of slits
in the barrier, resulted in an increased number of areas of both high and low probability
for particle detection. It also increased the variance of particle detection probability.
\end{document} \end{document}

View File

@ -2,6 +2,88 @@
\begin{document} \begin{document}
\section{Introduction}\label{sec:introduction} \section{Introduction}\label{sec:introduction}
% Light: wave particle
The nature of light has long been a subject of interest and discussion. In classical
mechanics, we study the kinematics and dynamics of physical objects, while ignoring
their intrinsic properties for simplicity. Elementary particles, such as photons
and electrons, do not abide by the laws of classical mechanics. A solution was
proposed by Max Planck with the radiation law, which he later derived using Boltzmann's
statistical interpretation of the second law of thermodynamics. Planck's findings
gave rise to the quantum hypothesis, and later Einstein's wave-particle duality \cite{britannica:1998:planck}.
The particle theory was the leading theory in the beginning of the 1800s, when
Thomas Young demonstrated the interference of light, through his double-slit experiment,
while postulating light as waves \cite{young:1804:double_slit}. The study of
interference of light, and Gustav R. Kirchhoff's study of ideal blackbodies,
showed that light exhibits both wavelike and particle-like characteristics. 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 particles, resulting in the wave equation.
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 barriers with single-slit and triple-slits. To solve the equation,
we will apply the Crank-Nicolson method in 2+1 dimensions.
In Section \ref{sec:methods}, we will present the theoretical background for
this experiment, as well as the methods 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} \end{document}
% 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
% crank-nicolson method!
% wave equation

View File

@ -1,7 +1,298 @@
\documentclass[../schrodinger_simulation.tex]{subfiles} \documentclass[../schrodinger_simulation.tex]{subfiles}
\begin{document} \begin{document}
\section{Methods}\label{sec:methods} \section{Methods}\label{sec:methods} %
\subsection{The Schrödinger equation}\label{ssec:schrodinger} %
Erwin Schrödinger wanted to find a mathematical description of the wave characteristics
of matter, which supported 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}. 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 unit, and $\hbar$ is the reduced Planck's constant. $\hat{H}$ is
a Hamiltonian operator, which represents 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 the Born rule, the square modulus of the wave function is proportional to the
probability density of detecting 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 conjugate of the wave function. When the potential
is time-independent, and the particle is non-relativistic, 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 are expressions of the kinetic energy, and the potental $V$
encodes 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 the Schrödinger equation by the dimensionful variables, we are left with
the wave function $u$ and the potential $v$. The dimensionless equation is given by
\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} %
As a result of working in position space, the Born rule is given by
\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}
where we assume a normalized wave function $u(x, y, t)$. We will initialize the wave
function, using a Gaussian wavepacket, given by
\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*}
$x_{c}$ and $y_{c}$ are the coordinates of the center of the wavepacket, $\sigma_{x}$
and $\sigma_{y}$ are the width of the wavepacket. The wave packet momenta are
given by $p_{x}$ and $p_{y}$.
\subsection{The Crank-Nicolson scheme}\label{ssec:crank_nicolson} %
When we evaluate a particle's position, we have to consider partial differential
equations (PDEs). To solve these numerically, we have to discretize Equation \eqref{eq:schrodinger_dimensionless}.
We use the $\theta$-rule \footnote{We can derive Forward Euler using $\theta = 1$, and Backward Euler using $\theta = 0$, in the $\theta$-rule.},
to combine the forward (explicit) and backward (implicit) finite difference methods.
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 the notation, and avoid any confusion of the indices with the imaginary unit $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 derive the Crank-Nicolson scheme (CN) by using $\theta = 1/2$, given by
\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_scheme}
\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_scheme} 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 the interference of light \cite{britannica:2023:young}, while postulating
light as waves rather than particles. The double-slit experiment results in a diffraction
pattern on a detector screen, where constructive interference of light result in
bright spots, and destructive interference result in dark spots. An illustration
of Thomas Young's setup can be found 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 barrier \cite[p. 4]{mit:2004:physics}.}
\label{fig:youngs_double_slit}
\end{figure}
After the wave passes through the barrier, the pattern observed is determined by
the path difference given by
\begin{align}
\delta = d \sin (\theta) = m \lambda \ ,
\label{eq:interference}
\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, while 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} %
In this experiment, we set up the grid with an equal step size in the x- and y-direction $h$,
and step size in the t-direction $\Delta t$, 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*}
In addition, we simplify the 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 results in a matrix $U^{n}$ that contains elements $u_{\ivec, \jvec}^{n}$, and
a matrix $V$ that contains elements $v_{\ivec, \jvec}$. We used Dirichlet boundary
conditions, given by
\begin{align*}
u(x=0, y, t) &= 0 & u(x=1, y, t) &= 0 \\
u(x, y=0, t) &= 0 & u(x, y=1, t) &= 0 \ ,
\end{align*}
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 are sparse and 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 use
\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 a pair of filled matrices can be found in Appendix \ref{ap:matrix_structure}.
For the general setup of the barrier, we used the values in Table \ref{tab:barrier_setup},
and for the simulations, we used the parameter settings in Table \ref{tab:sim_settings}.
% Insert Heisenberg uncertainty here? Or refer to it?
\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{Barrier parameters and values.}
\label{tab:barrier_setup}
\end{table}
\begin{table}[H]
\centering
\begin{tabular}{l r r} % @{\extracolsep{\fill}}
\hline
Parameter & Setting 1 & Setting 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{Simulation settings used in the double slit experiment. Setting 1 is
applied when the barrier is switched off and setting 2 is applied when the barrier
switched on.}
\label{tab:sim_settings}
\end{table}
To check if the total probability is conserved over time, and that the implementation
was correct, we computed the deviation from $1.0$ given by
\begin{align*}
s^{n} &= |1.0 - \sum_{\ivec , \jvec} p_{\ivec , \jvec}^{n}| \\
&= |1.0 - \sum_{\ivec , \jvec} u_{\ivec , \jvec}^{n*} u_{\ivec , \jvec}^{n}| \ .
\end{align*}
\subsection{Tools}\label{ssec:tools} %
The double-slit experiment is implemented in C++. We use the Python library
\verb|NumPy| \cite{harris:2020:numpy}, \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} \end{document}
% \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}

View File

@ -2,6 +2,127 @@
\begin{document} \begin{document}
\section{Results}\label{sec:results} \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 \verb|superlu| solver, which is a solver for sparse matrices. It is
generally used to solve nonsymmetric, sparse matrices. However, as the alternative
solver \verb|lapack| converts a sparse matrix to a dense matrix, it will increase
memory usage compared to \verb|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 used a solver for sparse matrices, we decrease the number of computations performed
compared to the number of computations using a solver for dense matrices.
We checked if the total probability was conserved over time, by plotting the deviation
from $1.0$.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/probability_deviation.pdf}
\caption{Deviation of total probability, for time $t \in [0, T]$ where $T=0.008$.}
\label{fig:deviation}
\end{figure}
We simulated the wave equation with the barrier switched off, using setting 1 in
Table \ref{tab:sim_settings} found in Section \ref{ssec:implementation}. When the
barrier was switched on, we used setting 2 in \ref{tab:sim_settings}. We observed
a larger deviation of total probability for a barrier with double slits compared
to no barrier. The result can found in Figure \ref{fig:deviation}. When the wave interacts
with the barrier, it results in a larger change in kinetic energy. The result is more prone
to computational errors, than if the wave propagates without interacting with a
barrier. No interaction results in a more stable deviation from the total probability.
In addition, we have to consider the limitations of the computer, therefore some computational
error is to be expected.
\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
We studied the time evolution of the probability function, using setting 2 in
Table \ref{tab:sim_settings}, found in Section \ref{ssec:implementation}. To visualize
the time evolution, we created colormap plots for different time steps. Figure \ref{fig:colormap_0_prob},
Figure \ref{fig:colormap_1_prob}, and Figure \ref{fig:colormap_2_prob} show the
results for time steps $t=\{0, 0.001, 0.002\}$, respectively. In addition, we created
separate plots for the real and imaginary part of $u_{\ivec, \jvec}$, for the same
time steps. The results can be found in Appendix \ref{ap:figures}, in Figure \ref{fig:colormap_real_imag}.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/color_map_0_prob.pdf}
\caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0$.}
\label{fig:colormap_0_prob}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/color_map_1_prob.pdf}
\caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0.001$.}
\label{fig:colormap_1_prob}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/color_map_2_prob.pdf}
\caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0.002$.}
\label{fig:colormap_2_prob}
\end{figure}
At time step $t=0.001$, Figure \ref{fig:colormap_1_prob}, when the wave interacts
with the double slit barrier, we observe a clear diffraction pattern in the
probability function. At time step $t=0$ (Figure \ref{fig:colormap_0_prob}) and
$t=0.002$ (Figure \ref{fig:colormap_2_prob}), the diffraction pattern is not as
clear. It is, however, more visible when we observe the real and imaginary part
separately in Figure \ref{fig:colormap_real_imag}. Since the probability function
is a product of $u_{\ivec, \jvec}$ and its conjugate $u_{\ivec, \jvec}^{*}$,
initialized by a Gaussian wavepacket, the result is a sum of the real and imaginary part.
% This can be found using Euler's formula, and the diffraction pattern is determined by interference given by \eqref{eq:interference}
In Figure \ref{fig:colormap_2_prob}, the probability function result in positive
areas at both sides of the barries. Some of the probability function is reflected
by the barrier, while the the rest spread out after passing the barrier. This is
a consequence of the wave-particle duality.
To compare the probability function $p_{\ivec, \jvec}$ for all the time steps in
a colormap plot, with the real and imaginary part of $u_{\ivec, \jvec}$, we created
Figure \ref{fig:colormap_all}. Where we excluded all x- and y-ticks, and labels, to
better visualize the diffraction pattern.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/color_map_all.pdf}
\caption{The time evolution of the probability function $p_{\ivec, \jvec}^{n}$ (top row),
Re($u_{\ivec, \jvec}$) (middle row), and Im($u_{\ivec, \jvec}$) (bottom row).
Time step $t=0$ in the left column, $t=0.001$ in the middle column, and $t=0.002$
in the right column.}
\label{fig:colormap_all}
\end{figure}
\subsection{Particle detection}\label{ssec:particle_detection}
% Problem 9: Plot detection probability for single-, double- and triple-slit
We simulated the wave equation using setting 2 in Table \ref{tab:sim_settings},
and assumed a detector screen located at $x=0.8$. To visualize the pattern of constructive
and destructive interference, we plotted the probability of particle detection,
along the screen, at time $t=0.002$. We adjusted the parameters to include single-, double-, and triple-slit
barriers. The results are found in Figure \ref{fig:particle_detection_single},
Figure \ref{fig:particle_detection_double}, and Figure \ref{fig:particle_detection_triple},
respectively.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/single_slit_detector.pdf}
\caption{Probability of particle detection along a detector screen at time $t=0.002$,
when using a single-slit barrier.}
\label{fig:particle_detection_single}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/double_slit_detector.pdf}
\caption{Probability of particle detection along a detector screen at time $t=0.002$,
when using a double-slit barrier.}
\label{fig:particle_detection_double}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/triple_slit_detector.pdf}
\caption{Probability of particle detection along a detector screen at time $t=0.002$,
when using a triple-slit barrier.}
\label{fig:particle_detection_triple}
\end{figure}
When the barrier has a single slit, there is no destructive interference and we
observe a single peak in the probability of particle detection. Adding another slit
result in more peaks, as there is both constructive and destructive interference.
When we used a triple-slit barrier, we observed an increase in interference which
resulted in narrow peaks. In addition, the probability of detecting a particle at
the ends of the screen increased with the number of slits.
\end{document} \end{document}

331
lib/WaveSimulation.cpp Normal file
View File

@ -0,0 +1,331 @@
/** @file WaveSimulation.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 0.1
*
* @brief Implementation of the WaveSimulation class.
*
* @bug No known bugs
* */
#include "WaveSimulation.hpp"
#include "utils.hpp"
#include <cmath>
#include <complex>
#include <cstdint>
#include <cstdlib>
#include <vector>
// Initializers
void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
double sigma_y, double p_x, double p_y)
{
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; });
}
void WaveSimulation::initialize_V()
{
this->V.set_size(this->N, this->N);
this->V.fill(0.);
}
void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
double ap, uint32_t slits)
{
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::initialize_A()
{
// Create the diagonal
arma::cx_vec diagonal(this->N * this->N);
// Set diagonal values
std::complex<double> r = (1._i * this->dt) / (2 * h * h);
for (size_t i = 0; i < diagonal.size(); i++) {
diagonal(i) = 1. + 4. * r + (1._i * this->dt / 2.) * this->V(i);
}
// Create the submatrix
arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros);
sub_matrix.diag(-1).fill(-r);
sub_matrix.diag(1).fill(-r);
// Set the size of A
this->A.set_size(this->N * this->N, this->N * this->N);
// Fill in the values in the submatrix diagonal
for (size_t i = 0; i < this->A.n_cols; i += this->N) {
this->A.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix;
}
// Fill the last sub/sup-diagonals
this->A.diag() = diagonal;
this->A.diag(-this->N).fill(-r);
this->A.diag(this->N).fill(-r);
}
void WaveSimulation::initialize_B()
{
std::complex<double> r = (1._i * this->dt) / (2 * h * h);
// Create the diagonal
arma::cx_vec diagonal(this->N * this->N);
for (size_t i = 0; i < diagonal.size(); i++) {
diagonal(i) = 1. - 4. * r - (1._i * this->dt / 2.) * this->V(i);
}
// Create the submatrix
arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros);
sub_matrix.diag(-1).fill(r);
sub_matrix.diag(1).fill(r);
// Set the size of B
this->B.set_size(this->N * this->N, this->N * this->N);
// Fill in the values in the submatrix diagonal
for (size_t i = 0; i < this->B.n_cols; i += this->N) {
this->B.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix;
}
// Fill the last sub/sup-diagonals
this->B.diag() = diagonal;
this->B.diag(-this->N).fill(r);
this->B.diag(this->N).fill(r);
}
// 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->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();
}
// 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';
}

17
lib/literals.cpp Normal file
View File

@ -0,0 +1,17 @@
/** @file literals.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief The implementation of the literals.
*
* @bug No known bugs
* */
#include "literals.hpp"
std::complex<double> operator""_i(long double magnitude)
{
return std::complex<double>(0.,magnitude);
}

View File

@ -45,14 +45,14 @@ bool mkpath(std::string path, int mode)
cur_dir = path.substr(0, pos); cur_dir = path.substr(0, pos);
if (mkdir(cur_dir.c_str(), mode) != 0 if (mkdir(cur_dir.c_str(), mode) != 0
&& stat(cur_dir.c_str(), &buf) != 0) { && stat(cur_dir.c_str(), &buf) != 0) {
return -1; return false;
} }
} }
else { else {
break; break;
} }
} }
return 0; return true;
} }
std::string dirname(const std::string &path) std::string dirname(const std::string &path)
@ -70,4 +70,47 @@ std::string concatpath(const std::string &left, const std::string &right)
} }
} }
void print_sp_matrix_structure(const arma::sp_cx_mat &A)
{
using namespace std;
using namespace arma;
// Declare a C-style 2D array of strings.
string S[A.n_rows][A.n_cols];
// Initialise all the strings to " ".
for (int i = 0; i < A.n_rows; i++) {
for (int j = 0; j < A.n_cols; j++) {
S[i][j] = " ";
}
}
// Next, we want to set the string to a dot at each non-zero element.
// To do this we use the special loop iterator from the sp_cx_mat class
// to help us loop over only the non-zero matrix elements.
sp_cx_mat::const_iterator it = A.begin();
sp_cx_mat::const_iterator it_end = A.end();
int nnz = 0;
for (; it != it_end; ++it) {
S[it.row()][it.col()] = "";
nnz++;
}
// Finally, print the matrix to screen.
cout << endl;
for (int i = 0; i < A.n_rows; i++) {
cout << "| ";
for (int j = 0; j < A.n_cols; j++) {
cout << S[i][j] << " ";
}
cout << "|\n";
}
cout << endl;
cout << "matrix size: " << A.n_rows << "x" << A.n_cols << endl;
cout << "non-zero elements: " << nnz << endl;
cout << endl;
}
} // namespace utils } // namespace utils

View File

@ -0,0 +1,96 @@
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():
ticks = [0, 0.25, 0.5, 0.75, 1.0]
with open("data/color_map.txt") as f:
lines = f.readlines()
size = int(lines[0])
for i, line in enumerate(lines[1:]):
# 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: complex(*ast.literal_eval(x)), arr)))
# Reshape and transpose array
arr = arr.reshape(size, size).T
# Plot color maps
color_map1 = ax1.imshow(
np.multiply(arr, arr.conj()).real,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True),
extent=[0, 1.0, 0, 1.0]
)
color_map2 = ax2.imshow(
arr.real,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True),
extent=[0, 1.0, 0, 1.0]
)
color_map3 = ax3.imshow(
arr.imag,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True),
extent=[0, 1.0, 0, 1.0]
)
# 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)
# Set custom ticks
ax1.set_xticks(ticks)
ax1.set_yticks(ticks)
ax2.set_xticks(ticks)
ax2.set_yticks(ticks)
ax3.set_xticks(ticks)
ax3.set_yticks(ticks)
# Set labels
ax1.set_xlabel("x-axis")
ax1.set_ylabel("y-axis")
ax2.set_xlabel("x-axis")
ax2.set_ylabel("y-axis")
ax3.set_xlabel("x-axis")
ax3.set_ylabel("y-axis")
# Save the figures
fig1.savefig(f"latex/images/color_map_{i}_prob.pdf", bbox_inches="tight")
fig2.savefig(f"latex/images/color_map_{i}_real.pdf", bbox_inches="tight")
fig3.savefig(f"latex/images/color_map_{i}_imag.pdf", bbox_inches="tight")
# Close figures
plt.close(fig1)
plt.close(fig2)
plt.close(fig3)
if __name__ == "__main__":
plot()

View File

@ -0,0 +1,62 @@
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,
"xtick.bottom": False,
"xtick.labelbottom": False,
"ytick.left": False,
"ytick.labelleft": False,
"axes.grid": False,
# "figure.autolayout": True
# "figure.constrained_layout.use": False,
# "figure.constrained_layout.h_pad": 0.04167,
# "figure.constrained_layout.w_pad": 0.04167
# "figure.subplot.wspace": 0,
# "figure.subplot.hspace": 0
}
plt.rcParams.update(params)
def plot():
fig, axes = plt.subplots(figsize=(6,6), nrows=3, ncols=3)
with open("data/color_map.txt") as f:
lines = f.readlines()
size = int(lines[0])
for i, line in enumerate(lines[1:]):
arr = line.strip().split("\t")
arr = np.asarray(list(map(lambda x: complex(*ast.literal_eval(x)), arr)))
# Reshape and transpose array
arr = arr.reshape(size, size).T
# Plot color maps
axes[0, i].imshow(
np.multiply(arr, arr.conj()).real,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True)
)
axes[1, i].imshow(
arr.real,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True)
)
axes[2, i].imshow(
arr.imag,
interpolation="nearest",
cmap=sns.color_palette("mako", as_cmap=True)
)
# Create tight subplots and save figure
fig.subplots_adjust(wspace=0, hspace=0)
fig.savefig(f"latex/images/color_map_all.pdf", bbox_inches="tight")
if __name__ == "__main__":
plot()

View File

@ -0,0 +1,55 @@
import numpy as np
import matplotlib.pyplot as plt
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 = [
"data/screen/single_slit.txt",
"data/screen/double_slit.txt",
"data/screen/triple_slit.txt",
]
outputs = [
"latex/images/single_slit_detector.pdf",
"latex/images/double_slit_detector.pdf",
"latex/images/triple_slit_detector.pdf",
]
colors = sns.color_palette("mako", n_colors=2)
for file, output in zip(files, outputs):
with open(file) as f:
lines = f.readlines();
size = int(lines[0])
column = int(.8 * size)
fig, ax = plt.subplots()
arr = lines[1].strip().split("\t")
arr = np.asarray(list(map(lambda x: complex(*ast.literal_eval(x)), arr)))
arr = arr.reshape(size,size)
x = np.linspace(0,1, size)
slice = arr[column, :]
norm = 1. / np.sqrt(sum([(i*i.conjugate()).real for i in slice]))
slice *= norm
slice = np.asarray([i * i.conjugate() for i in slice])
ax.plot(x, slice, color=colors[0])
ax.set_xlabel("Detector screen (y-axis)")
ax.set_ylabel("Detection probability")
fig.savefig(output, bbox_inches="tight")
plt.close(fig)
if __name__ == "__main__":
plot()

View File

@ -0,0 +1,56 @@
import numpy as np
import matplotlib.pyplot as plt
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()
img = ax.imshow([[]])
def plot():
with open("test.txt") as f:
lines = f.readlines();
size = int(lines[0])
for line in lines[1:]:
arr = line.strip().split("\t")
arr = np.asarray(list(map(lambda x: ((a := complex(*ast.literal_eval(x)))*a.conjugate()).real, arr)))
# print(sum(arr))
arr = arr.reshape(size,size)
wave_arr.append(arr.T)
# print(arr)
# plt.imshow(arr, cmap="hot", interpolation="nearest")
# plt.show()
def animation(i):
norm = matplotlib.cm.colors.Normalize(vmin=0, vmax=np.max(wave_arr[i]))
img.set_norm(norm)
img.set_data(wave_arr[i])
if __name__ == "__main__":
plot()
norm = matplotlib.cm.colors.Normalize(vmin=0, vmax=np.max(wave_arr[0]))
img = ax.imshow(wave_arr[0], extent=[0,1,0,1], cmap=plt.get_cmap("viridis"), norm=norm)
anim = FuncAnimation(fig, animation, interval=1, frames=np.arange(0,len(wave_arr)), repeat=True, blit=0)
# plt.show()
anim.save("./animation.mp4", writer="ffmpeg", bitrate=10000, fps=15)

View File

@ -0,0 +1,51 @@
import numpy as np
import matplotlib.pyplot as plt
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 = [
"data/probability_deviation_no_slits.txt",
"data/probability_deviation_slits.txt",
]
labels = [
"none",
"double-slit"
]
colors = sns.color_palette("mako", n_colors=2)
fig, ax = plt.subplots()
for file, label, color in zip(files, labels, colors):
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]))
ax.plot(x, arr, label=label, color=color)
ax.legend(title="Barrier")
ax.set_xlabel("Time (t)")
ax.set_ylabel("Deviation")
fig.savefig("latex/images/probability_deviation.pdf", bbox_inches="tight")
if __name__ == "__main__":
plot()

40
python_scripts/plot_v.py Normal file
View File

@ -0,0 +1,40 @@
import numpy as np
import matplotlib.pyplot as plt
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:
lines = f.readlines();
size = int(lines[0])
for line in lines[1:]:
arr = line.strip().split("\t")
arr = np.asarray(list(map(lambda x: ((a := complex(*ast.literal_eval(x)))*a.conjugate()).real, arr)))
# print(sum(arr))
arr = arr.reshape(size,size)
# print(arr)
plt.imshow(arr.T, cmap="hot", interpolation="nearest")
plt.show()
if __name__ == "__main__":
plot()

View File

@ -0,0 +1,75 @@
/** @file main.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief Implementation of the testing library
*
* @bug No known bugs
* */
#include "WaveSimulation.hpp"
#include "utils.hpp"
#include <fstream>
void probability_deviation()
{
WaveSimulation sim_no_slits(.005, 2.5e-5, .008, .25, .5, .05, .05, 200.,
0.);
WaveSimulation sim_slits(.005, 2.5e-5, .008, .25, .5, .05, .10, 200., 0.,
0.02, .5, .05, .05, 2);
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()
{
WaveSimulation sim(.005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0., 0.02,
.5, .05, .05, 2);
std::vector<double> times{0., .001, .002};
sim.solve("data/color_map.txt", times);
}
void detector_screen()
{
WaveSimulation *sim = new WaveSimulation(
.005, 2.5e-5, .002, .25, .5, .05, .20, 200., 0., 0.02, .5, .05, .05, 1);
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);
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);
sim->solve("data/screen/triple_slit.txt");
}
int main()
{
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);
return 0;
}

View File

@ -0,0 +1,4 @@
int main()
{
return 0;
}