"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%javascript\n",
"IPython.tab_as_tab_everywhere = function(use_tabs) {\n",
" if (use_tabs === undefined) {\n",
" use_tabs = true; \n",
" }\n",
" // apply setting to all current CodeMirror instances\n",
" IPython.notebook.get_cells().map(\n",
" function(c) { return c.code_mirror.options.indentWithTabs=use_tabs; }\n",
" );\n",
" // make sure new CodeMirror instances created in the future also use this setting\n",
" CodeMirror.defaults.indentWithTabs=use_tabs;\n",
" };\n",
"IPython.tab_as_tab_everywhere()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting Makefile\n"
]
}
],
"source": [
"%%file Makefile\n",
"hello: hello.c\n",
"\tcc -o hello hello.c"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"cc -o hello hello.c\n",
"Hello, world! atan(3)=1.249046\n"
]
}
],
"source": [
"%%bash\n",
"rm -f hello\n",
"make\n",
"./hello"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This was super simple, and for such a simple \"project\" of only one source-code file probably unnecessary, but we can very quickly make this Makefile more robust and demonstrate some of the more advanced features of make. In larger projects, makefiles save a lot of time.\n",
"\n",
"For details, consult this makefile tutorial."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting Makefile\n"
]
}
],
"source": [
"%%file Makefile\n",
"\n",
"OBJ = hello.o\n",
"\n",
"# these are pre-defined, but we can also change the default behaviour, if we have multiple compilers\n",
"CC = gcc\n",
"#CC = icc\n",
"CFLAGS = -O\n",
"LIBS = -lm\n",
"\n",
"# special abbreviations: $@ is the file to be made; $? is the changed dependent(s).\n",
"\n",
"hello: $(OBJ)\n",
"#\t@echo -n \" 2.link: \"\n",
"\t$(CC) $(CFLAGS) -o $@ $(OBJ) $(LIBS)\n",
"\n",
"# this is the default implicit rule to convert xxx.c into xxx.o \n",
".c.o:\n",
"#\t@echo -n \" 1.compile: \"\n",
"\t$(CC) $(CFLAGS) -c $<\n",
"\n",
"clean:\n",
"#\t@echo -n \" 0.cleanup: \"\n",
"\trm -f hello a.out *.o *~ core"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rm -f hello a.out *.o *~ core\n",
"gcc -O -c hello.c\n",
"gcc -O -o hello hello.o -lm\n",
"Hello, world! atan(3)=1.249046\n"
]
}
],
"source": [
"%%bash\n",
"make clean; make\n",
"./hello"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gaussian packet evolution: a C project"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"More homework! This is a week-long project to graphically illustrate what happens to a Gaussian wave packet (we will make it up by adding a few individual waves) impacting on a barrier. We want to observe how some of the packet is reflected, and some undergoes quantum-mechanical tunneling, so both a reflected and a transmitted packet emerges after the impact. We will run the program repeatedly, changing the number of wave making up the packet, the number of points to plot, the initial energy of the incoming packet, etc. The physics of this scattering process will get discussed next, but first let's set up the mechanics of the new project:\n",
"\n",
"\t- make a project subdirectory\n",
"\t
- create a basic skeleton of a program that can read and parse user input (though does not compute anything at first)\n",
"\t
- have a basic Makefile in place\n",
"
\n",
"\n",
"The second part of the homework is to try to figure out what is going on in an old Fortran program that does the same thing. We do not know either Fortran or C syntax, but hopefully the computational content, i.e. the algorithmic core, will be self-explanatory. The task is to \"port\", or convert this old code to C, and learn some basic C programming in the process."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"%%bash\n",
"rm -rf packet\n",
"mkdir packet"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/esternin/5P10/Lab2/packet\n"
]
}
],
"source": [
"%cd packet"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing packet.c\n"
]
}
],
"source": [
"%%file packet.c\n",
"/*\n",
" * packet.c\n",
" * packet - generate a Gaussian wavepacket impacting on an energy barrier.\n",
" *\n",
" * Completed: January.2018 (c) E.Sternin\n",
" * Revisions: \n",
" *\n",
" */\n",
"\n",
"#ifndef VERSION /* date-derived in Makefile */\n",
"#define VERSION \"2018.01\" /* default, that's when we first wrote the program */\n",
"#endif\n",
"\n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"\n",
"#define MAX_STR 256\n",
"\n",
"/* Global variables */\n",
"\n",
"static char whoami[MAX_STR] ; /* argv[0] will end up here */\n",
"static int verbosity; /* Show detailed information */\n",
"static char options[] = \"Vvhp:t:\"; /* command-line options, : means takes a value */\n",
"static char help_msg[] = \"\\\n",
"%s []\\n\\\n",
"\\t-V report version information\\n\\\n",
"\\t-v increase verbosity, may combine several\\n\\\n",
"\\t-p # number of points for the packet\\n\\\n",
"\\t-t # time since the beginning\\n\\\n",
"\\t-h help message\\n\\\n",
"\\n\\\n",
"e.g.\\tpacket -v -p 601 -t 20\\n\" \n",
";\n",
"\n",
"/*************************************************service routines************/\n",
"\n",
"void __attribute__((noreturn)) die(char *msg, ...) {\n",
" va_list args;\n",
" va_start(args, msg);\n",
" vfprintf(stderr, msg, args);\n",
" fputc('\\n', stderr);\n",
" exit(1);\n",
" }\n",
"\n",
"/***************************************************************** main *********/\n",
"int main(int argc, char **argv) {\n",
" int i,p;\n",
" double t;\n",
"\n",
"/*\n",
" * default values, may get changed by the command-line options\n",
" */\n",
" verbosity = 0; /* 0 = quiet by default, 1 = info, 2 = debug */\t\n",
" p = 25; /* default to 25 points in the packet */\n",
" t = 0;\n",
"\n",
" strncpy(whoami, argv[0], MAX_STR);\n",
"\n",
" while ((i = getopt(argc, argv, options)) != -1)\n",
" switch (i) {\n",
" case 'V':\n",
" printf(\" %s: version %s\\n\",whoami,VERSION);\n",
" break;\n",
" case 'v':\n",
" verbosity++;\n",
" if (verbosity > 1) printf(\" %s: verbosity level set to %d\\n\",whoami,verbosity);\n",
" break;\n",
" case 'h':\n",
" die(help_msg,whoami);\n",
" break;\n",
" case 'p':\n",
" if ( ((p = atoi(optarg)) > 10000) || p < 10 )\n",
" die(\" %s: -p %d is not a valid number of points (10..10000)\\n\",whoami,p);\n",
" if (verbosity > 1) printf(\" %s: Number of points = %d\\n\",whoami,p);\n",
" break;\n",
" case 't':\n",
" if ( ((t = atof(optarg)) < 0) )\n",
" die(\" %s: -t %d is not a valid time\\n\",whoami,t);\n",
" if (verbosity > 1) printf(\" %s: time = %f\\n\",whoami,t);\n",
" break;\n",
" default:\n",
" if (verbosity > 0) die(\" try %s -h\\n\",whoami);\t/* otherwise, die quietly */\n",
" return 0;\n",
" }\n",
"\n",
"/*\n",
" * when we get here, we parsed all user input, and are ready to calculate things\n",
" */\n",
"\n",
" for (i=0; i < p; i++) { \n",
" if (verbosity > 0) printf(\"%d\\t%f\\n\",i,i*t);\n",
" }\n",
"\n",
" return 0;\n",
" }\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" ./packet: verbosity level set to 2\n",
" ./packet: verbosity level set to 3\n",
" ./packet: version 2023-Sep-12-12:52\n",
" ./packet: Number of points = 10\n",
"0\t0.000000\n",
"1\t0.000000\n",
"2\t0.000000\n",
"3\t0.000000\n",
"4\t0.000000\n",
"5\t0.000000\n",
"6\t0.000000\n",
"7\t0.000000\n",
"8\t0.000000\n",
"9\t0.000000\n"
]
}
],
"source": [
"%%bash\n",
"cc -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" -o packet packet.c\n",
"./packet -vvvV -p10"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing Makefile\n"
]
}
],
"source": [
"%%file Makefile\n",
"\n",
"OBJ = packet.o\n",
"\n",
"# these are pre-defined, but we can also change the default behaviour, if we have multiple compilers\n",
"CC = cc\n",
"#CC = icc\n",
"CFLAGS = -O -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" \n",
"LIBS = -lm\n",
"\n",
"# special abbreviations: $@ is the file to be made; $? is the changed dependent(s).\n",
"\n",
"packet: $(OBJ)\n",
"\t$(CC) $(CFLAGS) -o $@ $(OBJ) $(LIBS)\n",
"\n",
"# this is the default implicit rule to convert xxx.c into xxx.o \n",
".c.o:\n",
"\t$(CC) $(CFLAGS) -c $<\n",
"\n",
"clean:\n",
"\trm -f packet *.o *~ core"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rm -f packet *.o *~ core\n",
"cc -O -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" -c packet.c\n",
"cc -O -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" -o packet packet.o -lm\n",
" ./packet: version 2023-Sep-12-12:52\n",
" ./packet: verbosity level set to 2\n",
" ./packet: verbosity level set to 3\n",
" ./packet: Number of points = 11\n",
" ./packet: time = 25.000000\n",
"0\t0.000000\n",
"1\t25.000000\n",
"2\t50.000000\n",
"3\t75.000000\n",
"4\t100.000000\n",
"5\t125.000000\n",
"6\t150.000000\n",
"7\t175.000000\n",
"8\t200.000000\n",
"9\t225.000000\n",
"10\t250.000000\n"
]
}
],
"source": [
"%%bash\n",
"make clean; make\n",
"./packet -Vvvv -p 11 -t 25"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have a skeleton of a C program working, let's examine the old Fortran code to see how it works and what its output looks like."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing packet.f\n"
]
}
],
"source": [
"%%file packet.f\n",
"C packet.f -- propagation of a Gaussian wave packet through a barrier in 1D\n",
"C\tx = coordinate (in units of a, the barrier is from -a to +a)\n",
"C\tP = wavefunction (not normalized)\n",
"C\tE = energy of the wavefunction (in units of V, the barrier height)\n",
"C\n",
"C\tWritten by:\tE.Sternin\n",
"C\tCompleted:\t30.03.93\n",
"C=======================================================================\t\n",
"\treal*4\t\tx(601),x0,p0,sp,p_now,p_step,Pi,A0,t\n",
"\tcomplex*8\tP(601),I,A\n",
"\tinteger\t\tj,Nargs,Nwaves\n",
"\tcharacter*80 \tbuffer\n",
"\n",
"\tNwaves = 20\t! add up (2*Nwaves+1) eigenwaves to form a packet\n",
"\tx0 = -25.\t! start the packet far away from the barrier\n",
"\tp0 = 0.5\t! keep incident energy low, E=0.5**2 < 1 (units of V)\n",
"\tsp = 0.1\t! p0 +/- 4 sigma_p > 0,\tall eigenwaves move to the right\n",
"\t \n",
"\tNargs = iargc() \t! arguments supplied on the command line ?\n",
"\tif (Nargs .lt. 1) then \n",
"\t write(*,'(A,$)') ' Error: need time t: '\n",
"\t read(*,*) t\n",
"\telse\n",
"\t call getarg(1,buffer)\n",
"\t read(buffer,*) t\n",
"\tend if\n",
"\n",
"\tI = cmplx(0.,1.)\t\t!complex i\n",
"\tPi = 2.*acos(0.)\n",
"\tA0 = 1./sqrt( sqrt(2.*Pi) * sp ) /float (2*Nwaves+1)\n",
"\n",
"\tdo j = 1,601\n",
"\t x(j) = -75 + (j-1)*0.25\n",
"\t P(j) = cmplx(0.,0.)\n",
"\tend do\n",
"\n",
"\tp_step = 4.*sp/float(Nwaves)\t! calculate over +/- 4*sigma_p\n",
"\tdo j = -Nwaves,Nwaves\t\t! in Nwaves steps\n",
"\t p_now = p0 + j * p_step \n",
"\t E = p_now**2\n",
"\t A = A0*exp( - ( (p_now-p0)/(2*sp) )**2 - I*(p_now*x0 + E*t) )\n",
"\t call add_wave(x,P,A,E)\n",
"\tend do\n",
"\n",
"C..output the saved values as is, without normalization\n",
" 10\tdo j=1,\t601\n",
"\t write(*,*) x(j),abs(P(j))**2,real(P(j)),aimag(P(j))\n",
"\tend do\n",
"\tend\n",
"\n",
"C=======================================================================\t\n",
"\tsubroutine add_wave(x,P,A,E)\n",
"C..adds the specified wave to the the supplied complex storage array P\n",
"C the values of Psi are only calculated at points given in real array x \n",
"\treal*4 \t\tx(601),E\n",
"\tcomplex*8\tP(601),I,A,k1,K2,C1,b,c,d,f\n",
"\n",
"\tI = cmplx(0.,1.)\t!complex i\n",
"\n",
"\tk1 = sqrt(cmplx(E,0.))\n",
"\tK2 = sqrt(cmplx(E-1.,0.))\n",
"\n",
"\tC1 = -2*k1*K2-k1**2-2*k1*K2*exp(4*I*K2)+exp(4*I*K2)*k1**2\n",
" *\t-K2**2+K2**2*exp(4*I*K2)\n",
"\n",
"\tb = (k1**2-K2**2)/ C1 *(exp(2*I*(2*K2-k1))-exp(-2*I*k1))\n",
"\tc = -2*k1*(K2+k1)/ C1 *exp(I*(K2-k1))\n",
"\td = (-2*K2+2*k1)*k1/ C1 *exp(I*(-k1+3*K2))\n",
"\tf = -4*k1*K2/ C1 *exp(2*I*(K2-k1))\n",
"\n",
"\tdo j = 1,296\t\t! left of the barrier, x < -1\n",
"\t P(j) = P(j) + A * ( exp(I*k1*x(j)) + b*exp(-I*k1*x(j)) )\n",
"\tend do\n",
"\tdo j=297,304\t\t! inside the barrier, -1 < x < 1\n",
"\t P(j) = P(j) + A * ( c*exp(I*K2*x(j)) + d*exp(-I*K2*x(j)) )\n",
"\tend do\n",
"\tdo j=305,601\t\t! to the right of the barrier, x > 1\n",
"\t P(j) = P(j) + A * f*exp(I*k1*x(j))\n",
"\tend do\n",
"\t \n",
" \treturn\n",
"\tend"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" -75.0000000 8.52197530E-08 1.16693496E-04 2.67586205E-04\n",
" -74.7500000 3.20557263E-07 3.57286714E-04 4.39207768E-04\n",
" -74.5000000 7.47516651E-07 6.27484871E-04 5.94793586E-04\n",
" -74.2500000 1.36120696E-06 9.13982512E-04 7.25150225E-04\n",
" -74.0000000 2.12284544E-06 1.20234303E-03 8.22931761E-04\n",
" -73.7500000 2.96266194E-06 1.47765456E-03 8.82722496E-04\n",
" -73.5000000 3.78836876E-06 1.72579114E-03 9.00007668E-04\n",
" -73.2500000 4.50573361E-06 1.93375722E-03 8.75395373E-04\n",
" -73.0000000 5.02329522E-06 2.09036935E-03 8.08486657E-04\n",
" -72.7500000 5.27878183E-06 2.18690187E-03 7.04444654E-04\n",
" -72.5000000 5.24243342E-06 2.21807300E-03 5.67965733E-04\n",
" -72.2500000 4.92594882E-06 2.18175026E-03 4.07325802E-04\n",
" -72.0000000 4.37627432E-06 2.07901793E-03 2.32291524E-04\n",
" -71.7500000 3.67039229E-06 1.91512029E-03 5.20246103E-05\n",
" -71.5000000 2.90051889E-06 1.69875333E-03 -1.21474615E-04\n",
" -71.2500000 2.15335467E-06 1.44088222E-03 -2.77872314E-04\n",
" -71.0000000 1.50178641E-06 1.15568750E-03 -4.07642918E-04\n",
" -70.7500000 9.90534545E-07 8.59129708E-04 -5.02424780E-04\n",
" -70.5000000 6.29514943E-07 5.67211304E-04 -5.54784900E-04\n",
" -70.2500000 4.03137477E-07 2.97086081E-04 -5.61139314E-04\n",
" -70.0000000 2.74815335E-07 6.43452368E-05 -5.20264381E-04\n",
" -69.7500000 2.00807236E-07 -1.17278614E-04 -4.32496192E-04\n",
" -69.5000000 1.46606908E-07 -2.35196581E-04 -3.02141474E-04\n",
" -69.2500000 9.81526753E-08 -2.82348512E-04 -1.35764451E-04\n",
" -69.0000000 6.73033540E-08 -2.52886268E-04 5.78955514E-05\n",
" -68.7500000 9.34146485E-08 -1.47164857E-04 2.67875264E-04\n",
" -68.5000000 2.34040129E-07 3.21054831E-05 4.82710428E-04\n",
" -68.2500000 5.51280777E-07 2.76413339E-04 6.89112814E-04\n",
" -68.0000000 1.09518010E-06 5.73823112E-04 8.75161262E-04\n",
" -67.7500000 1.88719878E-06 9.09983180E-04 1.02914008E-03\n",
" -67.5000000 2.90695652E-06 1.26707926E-03 1.14081835E-03\n",
" -67.2500000 4.09146287E-06 1.62594335E-03 1.20323373E-03\n",
" -67.0000000 5.33215916E-06 1.96641800E-03 1.21052016E-03\n",
" -66.7500000 6.49898357E-06 2.26926315E-03 1.16164889E-03\n",
" -66.5000000 7.44904264E-06 2.51625036E-03 1.05713122E-03\n",
" -66.2500000 8.05988930E-06 2.69161747E-03 9.02820029E-04\n",
" -66.0000000 8.24550898E-06 2.78319465E-03 7.06637395E-04\n",
" -65.7500000 7.98025758E-06 2.78387824E-03 4.79874376E-04\n",
" -65.5000000 7.28935311E-06 2.68957182E-03 2.35704472E-04\n",
" -65.2500000 6.26574547E-06 2.50312570E-03 -1.03411730E-05\n",
" -65.0000000 5.03992305E-06 2.23173853E-03 -2.43447372E-04\n",
" -64.7500000 3.76202502E-06 1.88760250E-03 -4.46073827E-04\n",
" -64.5000000 2.57725833E-06 1.48712832E-03 -6.04737783E-04\n",
" -64.2500000 1.60647710E-06 1.05198694E-03 -7.06965802E-04\n",
" -64.0000000 9.16446481E-07 6.04671310E-04 -7.42171891E-04\n",
" -63.7500000 5.25519908E-07 1.70712476E-04 -7.04540405E-04\n",
" -63.5000000 4.00786774E-07 -2.23876239E-04 -5.92170749E-04\n",
" -63.2500000 4.74573596E-07 -5.55795734E-04 -4.07019281E-04\n",
" -63.0000000 6.67948825E-07 -8.02147144E-04 -1.56552880E-04\n",
" -62.7500000 9.16697616E-07 -9.45987180E-04 1.47668005E-04\n",
" -62.5000000 1.18827541E-06 -9.73457354E-04 4.90567181E-04\n",
" -62.2500000 1.49675054E-06 -8.77579034E-04 8.52411729E-04\n",
" -62.0000000 1.90388926E-06 -6.57910947E-04 1.21286535E-03\n",
" -61.7500000 2.50078847E-06 -3.20466526E-04 1.54857663E-03\n",
" -61.5000000 3.39026064E-06 1.21302844E-04 1.83726603E-03\n",
" -61.2500000 4.65240237E-06 6.47406210E-04 2.05749064E-03\n",
" -61.0000000 6.31526882E-06 1.23210717E-03 2.19024671E-03\n",
" -60.7500000 8.33706690E-06 1.84477249E-03 2.22123414E-03\n",
" -60.5000000 1.05908302E-05 2.45194626E-03 2.13981071E-03\n",
" -60.2500000 1.28799593E-05 3.01780133E-03 1.94237835E-03\n",
" -60.0000000 1.49629741E-05 3.50721832E-03 1.63168437E-03\n",
" -59.7500000 1.65833972E-05 3.88620817E-03 1.21687388E-03\n",
" -59.5000000 1.75303558E-05 4.12542885E-03 7.14978145E-04\n",
" -59.2500000 1.76757167E-05 4.20163572E-03 1.48241408E-04\n",
" -59.0000000 1.70022122E-05 4.09815973E-03 -4.55302768E-04\n",
" -58.7500000 1.56223559E-05 3.80731258E-03 -1.06147316E-03\n",
" -58.5000000 1.37751495E-05 3.33144492E-03 -1.63603912E-03\n",
" -58.2500000 1.17840091E-05 2.68272613E-03 -2.14172597E-03\n",
" -58.0000000 1.00079678E-05 1.88314915E-03 -2.54199072E-03\n",
" -57.7500000 8.79676918E-06 9.65245708E-04 -2.80447328E-03\n",
" -57.5000000 8.41458677E-06 -3.08033777E-05 -2.90062721E-03\n",
" -57.2500000 9.01049771E-06 -1.05794892E-03 -2.80913548E-03\n",
" -57.0000000 1.05928020E-05 -2.06299732E-03 -2.51730881E-03\n",
" -56.7500000 1.30348853E-05 -2.99088843E-03 -2.02224450E-03\n",
" -56.5000000 1.61138578E-05 -3.78680765E-03 -1.33189501E-03\n",
" -56.2500000 1.95610373E-05 -4.39817924E-03 -4.65892372E-04\n",
" -56.0000000 2.31367121E-05 -4.77899890E-03 5.45784831E-04\n",
" -55.7500000 2.66794759E-05 -4.89073480E-03 1.66138215E-03\n",
" -55.5000000 3.01715318E-05 -4.70670173E-03 2.83169374E-03\n",
" -55.2500000 3.37367819E-05 -4.21285350E-03 3.99858039E-03\n",
" -55.0000000 3.76395037E-05 -3.41027766E-03 5.09995222E-03\n",
" -54.7500000 4.22198063E-05 -2.31518969E-03 6.07121922E-03\n",
" -54.5000000 4.78313013E-05 -9.60440026E-04 6.84900396E-03\n",
" -54.2500000 5.47439013E-05 6.05673355E-04 7.37408036E-03\n",
" -54.0000000 6.30851209E-05 2.32096715E-03 7.59593491E-03\n",
" -53.7500000 7.27457809E-05 4.11053933E-03 7.47323502E-03\n",
" -53.5000000 8.34101229E-05 5.89031586E-03 6.97956327E-03\n",
" -53.2500000 9.45698048E-05 7.57003017E-03 6.10446092E-03\n",
" -53.0000000 1.05608735E-04 9.05702077E-03 4.85583209E-03\n",
" -52.7500000 1.15908770E-04 1.02605270E-02 3.26042413E-03\n",
" -52.5000000 1.24996470E-04 1.10966908E-02 1.36378640E-03\n",
" -52.2500000 1.32657777E-04 1.14920577E-02 -7.68357306E-04\n",
" -52.0000000 1.39047828E-04 1.13893524E-02 -3.05458345E-03\n",
" -51.7500000 1.44682970E-04 1.07490271E-02 -5.39827580E-03\n",
" -51.5000000 1.50450665E-04 9.55420081E-03 -7.69206742E-03\n",
" -51.2500000 1.57490605E-04 7.81247579E-03 -9.82119329E-03\n",
" -51.0000000 1.67064078E-04 5.55775221E-03 -1.16694244E-02\n",
" -50.7500000 1.80347270E-04 2.84943869E-03 -1.31235654E-02\n",
" -50.5000000 1.98269292E-04 -2.27141310E-04 -1.40789812E-02\n",
" -50.2500000 2.21380440E-04 -3.56273819E-03 -1.44460145E-02\n",
" -50.0000000 2.49718258E-04 -7.02941464E-03 -1.41529366E-02\n",
" -49.7500000 2.82896624E-04 -1.04816696E-02 -1.31541332E-02\n",
" -49.5000000 3.20124294E-04 -1.37642613E-02 -1.14310738E-02\n",
" -49.2500000 3.60413804E-04 -1.67173911E-02 -8.99681263E-03\n",
" -49.0000000 4.02785983E-04 -1.91835742E-02 -5.89715503E-03\n",
" -48.7500000 4.46478196E-04 -2.10139658E-02 -2.21166899E-03\n",
" -48.5000000 4.91187500E-04 -2.20769364E-02 1.94842578E-03\n",
" -48.2500000 5.37160086E-04 -2.22638827E-02 6.44046534E-03\n",
" -48.0000000 5.85187867E-04 -2.14950033E-02 1.10974200E-02\n",
" -47.7500000 6.36558922E-04 -1.97261460E-02 1.57301631E-02\n",
" -47.5000000 6.92870934E-04 -1.69525705E-02 2.01365650E-02\n",
" -47.2500000 7.55716639E-04 -1.32109318E-02 2.41078399E-02\n",
" -47.0000000 8.26478878E-04 -8.58242996E-03 2.74375789E-02\n",
" -46.7500000 9.05986526E-04 -3.19017912E-03 2.99300738E-02\n",
" -46.5000000 9.94412228E-04 2.80074077E-03 3.14096808E-02\n",
" -46.2500000 1.09125383E-03 9.18914936E-03 3.17303203E-02\n",
" -46.0000000 1.19533867E-03 1.57422144E-02 3.07818335E-02\n",
" -45.7500000 1.30525301E-03 2.22040173E-02 2.84997337E-02\n",
" -45.5000000 1.41946936E-03 2.83027254E-02 2.48681568E-02\n",
" -45.2500000 1.53699284E-03 3.37631181E-02 1.99259799E-02\n",
" -45.0000000 1.65757514E-03 3.83147635E-02 1.37678627E-02\n",
" -44.7500000 1.78213743E-03 4.17049080E-02 6.54508593E-03\n",
" -44.5000000 1.91278907E-03 4.37084027E-02 -1.53771020E-03\n",
" -44.2500000 2.05288385E-03 4.41394821E-02 -1.02269184E-02\n",
" -44.0000000 2.20666802E-03 4.28607017E-02 -1.92257129E-02\n",
" -43.7500000 2.37896340E-03 3.97920758E-02 -2.82055680E-02\n",
" -43.5000000 2.57440354E-03 3.49168628E-02 -3.68132629E-02\n",
" -43.2500000 2.79705203E-03 2.82865763E-02 -4.46869284E-02\n",
" -43.0000000 3.04991286E-03 2.00244263E-02 -5.14678098E-02\n",
" -42.7500000 3.33440048E-03 1.03225214E-02 -5.68141378E-02\n",
" -42.5000000 3.65031790E-03 -5.59103151E-04 -6.04152754E-02\n",
" -42.2500000 3.99616593E-03 -1.22994892E-02 -6.20071627E-02\n",
" -42.0000000 4.36929567E-03 -2.45239753E-02 -6.13829829E-02\n",
" -41.7500000 4.76662349E-03 -3.68160494E-02 -5.84054962E-02\n",
" -41.5000000 5.18542528E-03 -4.87304702E-02 -5.30166626E-02\n",
" -41.2500000 5.62385516E-03 -5.98065779E-02 -4.52441014E-02\n",
" -41.0000000 6.08159089E-03 -6.95862025E-02 -3.52044180E-02\n",
" -40.7500000 6.56005880E-03 -7.76285455E-02 -2.31055636E-02\n",
" -40.5000000 7.06258509E-03 -8.35293457E-02 -9.24300961E-03\n",
" -40.2500000 7.59381475E-03 -8.69352892E-02 6.00590743E-03\n",
" -40.0000000 8.15933105E-03 -8.75608996E-02 2.21905336E-02\n",
" -39.7500000 8.76477081E-03 -8.52020681E-02 3.87991965E-02\n",
" -39.5000000 9.41498019E-03 -7.97474012E-02 5.52750565E-02\n",
" -39.2500000 1.01132030E-02 -7.11869895E-02 7.10324943E-02\n",
" -39.0000000 1.08607309E-02 -5.96186668E-02 8.54771659"
]
},
{
"data": {
"text/html": [
"limit_output extension: Maximum message size of 10000 exceeded with 41469 characters"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%bash\n",
"gfortran -o packet packet.f\n",
"./packet 40"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Perhaps it's better to just see things graphically. gnuplot is a very good, universally available plotting program, and with just a few keystrokes we can generate good-looking graphs."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing Vo.dat\n"
]
}
],
"source": [
"%%file Vo.dat\n",
"-75\t0\n",
"-1\t0\n",
"-1\t1\n",
"1\t1\n",
"1\t0\n",
"75\t0"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"%%bash\n",
"gnuplot --persist\n",
"plot 'Vo.dat' with lines"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"%%bash\n",
"gnuplot --persist\n",
"set style data line\n",
"plot [-80:80] [-0.1:1.1] 'Vo.dat' title \"energy barrier V\",'<./packet 0' title \"packet at time = 0\",'<./packet 20' title \"packet at time = 20\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can even create an \"animation\" of sorts, by using a macro file."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing packet.gnu\n"
]
}
],
"source": [
"%%file packet.gnu\n",
"set style data line\n",
"set yrange [0:2]\n",
"set xrange [-75:75]\n",
"set xlabel \"x/a\"\n",
"set ylabel \"Psi^2\"\n",
"set title \"Scattering of a Gaussian wave packet (41 waves), E=V/4\"\n",
"plot '<./packet 0' t \"time = 0\",'Vo.dat' t \"V(x)\"\n",
"pause 2 \"initially, the packet is at x = -25 a\"\n",
"plot '<./packet 5' t \"time = 5\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"as time goes on...\"\n",
"plot '<./packet 10' t \"time = 10\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"...the packet approaches the barrier\"\n",
"plot '<./packet 15' t \"time = 15\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"the interference is observed\"\n",
"plot '<./packet 20' t \"time = 20\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"the interference is observed\"\n",
"plot '<./packet 25' t \"time = 25\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"the interference is observed\"\n",
"plot '<./packet 30' t \"time = 30\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"the interference is observed\"\n",
"plot '<./packet 35' t \"time = 35\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"transmitted and reflected packets form\"\n",
"plot '<./packet 40' t \"time = 40\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"transmitted and reflected packets form\"\n",
"plot '<./packet 45' t \"time = 45\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"and move away from the barrier\"\n",
"plot '<./packet 50' t \"time = 50\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"and move away from the barrier\"\n",
"plot '<./packet 55' t \"time = 55\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"note how they spread out with time\"\n",
"plot '<./packet 60' t \"time = 60\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"note how they spread out with time\"\n",
"plot '<./packet 65' t \"time = 65\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"note how they spread out with time\"\n",
"plot '<./packet 70' t \"time = 70\",'Vo.dat' t \"V(x)\" \n",
"pause 2 \"note how they spread out with time\"\n",
"plot '<./packet 80' t \"time = 80\",'Vo.dat' t \"V(x)\" \n",
"pause 2\n",
"plot '<./packet 90' t \"time = 90\",'Vo.dat' t \"V(x)\"\n",
"pause 2"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"initially, the packet is at x = -25 a\n",
"as time goes on...\n",
"...the packet approaches the barrier\n",
"the interference is observed\n",
"the interference is observed\n",
"the interference is observed\n",
"the interference is observed\n",
"transmitted and reflected packets form\n",
"transmitted and reflected packets form\n",
"and move away from the barrier\n",
"and move away from the barrier\n",
"note how they spread out with time\n",
"note how they spread out with time\n",
"note how they spread out with time\n",
"note how they spread out with time\n"
]
}
],
"source": [
"%%bash\n",
"gnuplot\n",
"load 'packet.gnu'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another way of using gnuplot is to load it right into the notebook itself."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"%load_ext gnuplot_kernel"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"%gnuplot inline pngcairo font \"Arial,16\" size 640,480"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAIAAAC6s0uzAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1wUd/4/8PcAS9llVzqCgkoxoGBA7CCCBVs05qL56sXYUX9RicZ4ycUYiTlN0YgmlhjbqWchsRfEgr2gYsMoKAoqKlV6h935/TF3e3sLLAPsMri8nn/4kM+09+zs7mtn5jMzDMuyBAAAAE3LQOgCAAAAWiIEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAWnoA37t3LzQ01N3dXSwWS6VST0/PWbNmPXz4UBfLevz48c6dOzW3NMzs2bMZhjl69GjjZ6UV169f79Onj0QisbS0/Prrr7U+/+Li4m3btg0bNszV1dXMzKxVq1be3t5z5859/Pix1pelLc1tG7Uc2nrl+X9aN2/ebGFhkZmZWX3Qw4cPJRJJQECAhsmnT59ubm7ewCqbvQ0bNtjb26enp9d3woMHDzIaDRkyhOeshg0bxjDMjz/+WNsId+7cYRjG3d2dZVllY0hIiL+/f/WR+WzT2rToAN66dauPj8+mTZsKCwu9vb09PDwyMzPXrVvn5eW1Y8cO7S7r7t27Xl5eZ8+e1dCiH+Ry+ciRI69evdqxY8e+fft26tRJu/O/evWqt7f3pEmTTp48aWlpOWjQoB49erx69Wr16tWenp4bNmzQ7uIAqD6f1uTk5LCwsPDwcDs7O7VBVVVV48ePLykp0TD5qlWrNm7c2Kham7fQ0NA2bdpMmTKlYZNbWVn1q0WXLl14zmT69OlEpOFLfuvWrVypDMNwLUVFRefPnx8+fLjamHy2qSZsS/Xy5UsTExNTU9Pff/9doVBwjZWVlevXrzcwMBCJRMnJyVpcHPfRnTp1qoaWBisoKEhLSysrK2v8rBrv+fPnRNSpUyddzDw+Pt7ExISIPvnkk6ysLGV7VVXVli1bzMzMGIa5cOGCLhbdSM1qG7Uos2bNIqIjR440Zib8P61Dhgzp0KFDRUVF9UGLFi3ivnX9/f2rDy0rK+NKJSKJRNKYapu506dPE1FkZGS9pjpw4AARDR8+vPEFVFZWtm7dmohu375dfWh5ebm1tbVIJMrIyFA27t+/n4ju3LmjNrLmbVqnlrsHfPLkyfLy8smTJ48ZM0b5M8fIyGjmzJnTp0+vrKzcvn27sBXyJ5VKW7duzSWT4MrLy4nIyclJ63NWKBR//etfy8vLv//++1WrVtnY2CgHGRoaTp48+ddff2VZ9ptvvtH6ohuvWW0j0JFz585FR0eHhYWJRCK1QbGxscuWLRs1alSNE0ZFRXl6eq5duzYwMFAqleq+UiENGDDAx8dn4cKFcrlckAKMjIwmT55MRDV+yR85cuT169ejRo1SPYZx7NixNm3avP3226pjat6mvDQgtPXDli1biGjs2LHVB927d++nn35S3ZGqqKj46aeffH19zc3NbW1tBwwYcOLECdVJnj17NmfOHA8PD7FYbGJi0r59+xkzZrx8+ZIbOnHiRNXXfPHixdVblAtatWqVr6+vWCw2Nzf39/ffs2eP6oKmTp1KRLdv3+7Xr5+JiUm7du1u3bql9ht/xowZRPTw4cMFCxZ06NDB2NjY0dFx5syZmZmZqrPKycn57LPPOnToYGJi4u7uvmLFioMHDxJRRESEhtetpKRk6dKlXl5epqamUqnU399/+/btyqF+fn6q69WuXbva5qP5FasR98PZw8NDLpfXOIJcLg8NDV2/fn29FvThhx8S0cWLF1Wnmj9/PhH98ccfypaNGzf26dPH2tra1NT0rbfeWrBgwevXr1Un0TxC9f2wOgvjuR1V9ejRg4geP36sWhUR+fn5qY4WFBRkYGDA/cDXXMZnn31GRD///LPagt555x0iunz5Mvdnne/b6ri1i4uLmzx5sq2trZmZmZeX1+rVq9U2Lp/3ieaPp9orX1paOmjQICKaOHGiclma66/t01rdsGHDjIyM1N4YLMsWFRW5ubl17NiR619SfW/p/ffft7S0XLlyZVVVlbW1dZ17wM15Q9f5SWFZdvXq1Wqfrzrx3wNWi0k1p06dYln2yZMnDMO0bt26qqpKbXLuODM3GkehUDg4OEyfPl11tDq3KR8tN4DT09PFYjERzZgx48aNG8qj0NWVlZX17t2biNzd3T/++OOZM2daWFgQ0datW7kRHj16ZG1tbWBgMGzYsPnz50+fPr1Dhw5E5Orqyh1yPHPmzLx587iNtH79+hs3blRvYVm2vLx8wIAB3IQff/zxJ5980r59eyKaN2+eshgugN3d3du1a/fee+95enqWl5fXGMBubm7u7u7Lly//5z//OWzYMCLq2bOncj7Z2dmenp5E1LVr1zlz5nD9Fzp27Kg5gHNycnx8fLgCZsyYMX78eG43dMKECdwI+/btW7JkCRF16tRp/fr1//rXv2qcT52vWI1mzpyp+euvYQviE8DLli3jsv+TTz5ZsGAB9/XXpUuXyspKniOobSM+hfHZjmqWLl1KRGvXrlW2jBkzhogMDAxycnK4lry8PCMjo4CAAD5lPHjwgIi6d++uupTMzEwjI6OOHTtyf/J531bHrZ2dnV3Xrl0PHTp04cIF7r09ZsyYem2+Oj+eqq98aWlpSEgIEc2cOVP5ka+z/ho/rdXl5uYaGhoGBgZWHxQaGmpoaBgbG5uamlrjl/XZs2fz8vK4//MJ4Ga7oev8IHC4/pLvv/++5tVUxT+A58yZM7x2ysPI3LocP35cddpXr14ZGhq6urqqJkJcXBwRHT58WHXMOrcpHy03gFmWjY6Otra25n4WWVtbjxgx4ocffrh586baaFw/3rFjxyo/848fP5bJZK1atSopKWFZdty4cUSkuuNVWVnJ7QuePn2aa+FzDnjx4sVE9O6775aWlnItJSUl/fr1U32XcF9SnTp1Ki4uVk5YYwB36dJFOR+WZXv16kVEV69e5f4MDQ0lojlz5ijfZ7///jv3UmgIYO64zUcffaR8KdLS0ri+D5s3b+ZakpKSiGjw4MG1zYTnK1Zd//79iejgwYMa5tyABfEJYHNz8zZt2ijXWqFQDB06lIiOHTvGcwS1bcSnMD7bUc2ff/5JRCNHjuT+lMvlVlZW3CHNffv2cY179uwhop9++olnGdwSExMTleOsWrWKiL777jvuTz7v2+q4tWvbtm1hYaGykesdozw7yKe8Oj+eyle+rKxs8ODB1X8Z8Kmfzzngffv2EdHf/vY3tfYjR44Q0cKFC1mW5fNlzSeAm+2GrvODoGRjY2NpaVnb0azquAC2sbEZUIvnz5/znJXq6zNu3DjVxh9++EF1lTnh4eEmJiaqX7n13aa1adEBzLJsTk5ORERE//79TU1N6T88PT1VvzicnZ0NDQ3Vjvtt3br1hx9+4I7tHD9+PCIiQq3bBfcNvnv3bu7POgNYoVDY2dkxDKN65p9l2cuXLxPRqFGjuD+5AFZ7f9QYwGo5+vnnnxMRd7i4tLRULBbb2NgUFRWpjjNy5EgNAVxUVGRsbCyTyVTfiCzLnjt3joi6devG/ckngPm8YtV17tyZiM6fP6/a+PDhw+ofReWPXD4LqjOAFQqFubm5mZmZapeNrKws5ZaqcwS22jbiU1id27FGrq6u5ubm3JyvXbtGRF9++SURzZo1ixth/PjxRPTkyROeZfz2229E9NVXXylH8PPzMzQ05I5e8nzfVsetndoxT+7NM2TIEP6vUp0fT+6V37t3L5cEX375peqYPOvnE8BffPFF9U2TmZlpb2/v6+vLrYW2AphtlhuazwdBKTg4mIgePnxY55pyuADWICEhgeesOOXl5TY2NmZmZgUFBcpGT09PkUiUnp6uOmb37t2V70m2Qdu0NkaaV0nvWVpazp07d+7cuRUVFTdu3Dh37ty+fftu3749fPjwyMjI0aNHZ2dnP3/+3NPT09bWVnXCSZMmKf8/ZMiQIUOGlJaW3r1798mTJ0+ePImPj+dOWPLvZfD06dPMzEyZTLZu3TrV9oqKCiLijoEoeXh41DlDV1dX1T+5Kwu5ucXHx5eUlPj7+0skEtVxgoKCDh8+XNsM7927V1FRERQUxB26VwoICBCJRNxHTtmdTbOGvWLc4Yrs7GzVxoKCgpiYGLUxc3NzG7MgNQzDhIWFLVu2zNfX18vLa8CAAYMGDQoODla+DnWO0JhXQMN2rNG77767cuXKy5cvBwUFnTx5kohmzJixbds27lVSKBTR0dHe3t4uLi48yxg7duzcuXN37ty5ZMkShmESEhJu3rw5bNgwR0dHquf7tjq1TgNubm5isfjWrVs8XyU+H0/OJ5988vLlS4ZhuF+ZSo2sX1VGRgYRqV19FBoampeXFxMTU71bViM1ww1drw8C90JlZGRwZ754Gj58eJ2XdIeFhSUnJ9c2dOnSpdxJYmNj4wkTJqxcuXLfvn3cGyY2NjYhIWH06NH29vbK8TMyMuLi4n7++Wdliza3aQNCW2+oHvtSxXUQcHZ2ZlmWO8Heq1cvDfN5/fp1aGiomZkZ95KamZkFBAR069aNiHbs2MGNU+cesOaPulgs5kbj9oBVOwiwtewBqx39+/bbb4lo48aNLMseO3aMVE7cKkVGRlLte8DR0dFE9MEHH1QfxPV55o748dkD5vOKVcddO6j5HDB3jPHs2bP8F8SzE9auXbuCg4OVHzmJRDJnzhzVM9aaR1DbRnwKq3M71uj8+fNE9MUXX7AsGxAQ4ObmxrLsRx99REQvX768dOkSEX399df8y2D/0wuJe4m4XfC9e/dyg3i+b6vj1q76xX7W1tbGxsY8y+Pz8VRe2/P+++8TkYeHh+ohfZ7189kD5nrDqr6RNm3aRETLly9XtmhxD7jZbug6Pykc7quM/xkl7XbC4iQkJBBR//79uT+5MyAnT55UnRvXXVf5Rm3YNq1NC70MiWVZS0tLW1vbGvckwsLCHBwcnj9/npeXx51Wyc/PVxunrKxM+ctxxIgRGzdu7Nev38GDB5OSkgoLCy9evFjfG6NwC/L29q5xOxUXFzdkPWshk8mIqKCgQK29eouqVq1aEdGLFy/U2lmWzcvLE4vFyo93nRr2inHdTP71r39VVVVpcUHcXrvafmdRUZHarMaNG3fmzJmcnJzo6OhPP/1ULBb/8ssvCxcu5D9C418BPvz9/a2trU+cOFFUVBQbGztw4EAi4v49c+YMt/egvHCCZxncT5+dO3eyLLtz504bG5sRI0Zwgxr5vlX7ZFVUVOTn5yt3Z+ssj8/Hk7No0aK9e/d+8MEHiYmJX331ldoctPK5497/eXl5ypbdu3cT0YIFC5R3a+J+ql6+fJlhGC8vL/4zr67ZbmieHwTuMJWGo0QNVv1qXVXcS8Tx8PAICAg4d+5campqaWlpZGSki4uL6ghEdPTo0U6dOnF91kjb27SFHoJmGMbb2/vixYtbt27lfomrysvLy8vLMzMzk8lkFhYWdnZ2SUlJr1+/VvbYIqIvv/zyl19+OXjwoIuLy5UrV1xcXKKiolQPwCYmJhIR+587mVU/NqvW4urqKpPJHjx4kJaW5uDgoGxPTU1duXKlr6/vhAkTtLHqREQ+Pj7GxsaxsbFyudzQ0FDZzp3OqY23t7dIJLp3715OTo6VlZWyPS4urrCwkOvBwUdCQgKfV6y6kJCQzp07379//7PPPuO6h6hhWbawsLC+C+KuzVX93iQirpMLJzU1dePGjc7OztOmTTM3Nx88ePDgwYPHjRvXvXt3bseozhG09QrwYWho+M4772zfvn3fvn1VVVVcV0/l9/L169fbtWvn6+tbrzICAwPd3NwOHDgwZsyYFy9efPLJJ8bGxtygRr5vL126xPWr51y4cKGqqqpPnz48y3NwcND88VTeuojrjvvLL7/ExMRERES899573G0FedbP59yKs7MzEWVlZam+blyXbKXS0tKoqChra+ugoCBu/AZrhhu6Xh8E7oVq165dY16ExgsNDb106dIff/zh5OSUn5//+eefq75ElZWVp06d4q6/4Gh5m2r4paDfoqOjGYYRiUTLli1TdtxnWfbRo0dc7wBlXwbuUMzEiROVXRhSUlIsLS2lUmlhYeGTJ0+IyN7eXrVrUmRkJLcVlYcKuWwbPXq0cpzqLdylDkOGDFH2jaqsrOR+gSqPu2rlEDT7n8O5X3zxhbIX9OnTp7kPm4Ze0NyZkg8//FD5UmRlZXFfbb/99hvXUuchaJ6vWI0ePHjAnbceNmzY9evXle2lpaX79+/v2bMnERkaGnJd2XkuKCIiQm1D7Nq1i/t0cIegc3NzTUxMbG1tVS8/5TqNc9dR1DkC+7/biGdhDTsEzf7nrj3t27c3MDBQXoLp6enJfXF88sknXEu9NgR33UvXrl2p2h4Gn/dtddzatWnThuslxLJsXl4eF8Zcv1ye5Wn+eLLVPh3cDQjd3d2Vs+VTf/VPa3XcGRzly1sjLR6CZpvfhubzQeDI5XILCwuZTMZdhiuXy7OyslTvbVedFu+EpaqkpMTCwsLf33/06NFGRkZpaWmqQ0+dOkVE586d0zAH9IJuoJ9//tnIyIj7yvbw8OjZsyd3WRsRBQcHK9+sxcXF3bt3J6LOnTuHhYVNmzZNJpMxDHPgwAFuBO7yGG9v7/Dw8CVLlnCnIbnT+MruytwNGg0MDPr167dmzZoaW4qLi7n9SGdn5+nTp8+bN++tt94iooCAAOU7XlsBnJWVxfXu6datW1hY2MiRIw0MDLhuEdUvxlfKzc3lLjry8PCYPXv25MmTudX88MMPlUHO5xwwn1esNtxdebnNZGtr27Vr186dOyuPfr/99tuxsbH1WlB2djZ3TL5r165Tp0719/dnGIa7/4DyHDB3cYKlpeWUKVO++OKL0aNHi0QimUwWHx/PcwS1bcSnsAYHcFFREderv2vXrsrGOXPmcC+R8gR5vTbEy5cvuYMlqvPk8HnfVsetnY2NTatWrSZPnjx37ty2bdsS0fz58+tVXp0fz+q3QOEupw4LC+Nff/VPa3VZWVlGRkY+Pj61rTKr7QBuhhu6zg8C5/bt26o/aLiXhYhUT8+r4QJYw72g+/Xrp7yWul5mzZrFMIypqelf/vIXtUFz585t1aqV2kXMahDADZeQkBAWFubt7W1lZSUSiezt7YcOHbpjxw61+3KUlJR8++23nTp1MjExkUqlgwYNOnPmjHJobm7up59+6uLiYmJi0qZNm8DAwI0bN3LdQ/r06aMcbcWKFQ4ODiYmJpMnT66tpays7Mcff/Tx8eEOgPv4+KxYsUL1d6u2Aphl2YyMjJkzZzo6OopEoo4dO65evZobR3lFb42Kioq++eYb7qWwtLQMCgratWuX6gh8ApjnK1ab8vLyvXv3jh49umPHjhKJRCKRuLu7T5o06dixY2objueC7t27N2rUKAsLC7FY3KdPn+PHj3NnelQ7Yf3zn//09/eXSqVGRkZt27adPHmy6n2I6hxBbRvxKazBAcz+5wZGqtekHjp0iIisra1Vb/1Trw3B5dYvv/xSfXF1vm+r49bu4MGD8+bNs7W1lUgkvXv3VrutEs/yNH88qwfw8+fPpVIpwzDKPRs+9Vf/tFb37rvvGhgYpKam1jaCdgOYbZYbus5PCvuf/WzlRuEfwJpp3oeuzZ07d7jJo6Oj1Qa5u7vX2O1UVWMCmGEbd8IJ3lBJSUlt27ZV6zY1e/bstWvXnjp1Sq0bAoDWzZw5c8OGDcePH+f/FLnm78aNGz169Pj666+b5w3JmwmFQuHq6mphYXHr1i3lCdcXL15wt+XieSmjfmihvaBh6NChrVq1evbsmbIlOTl5165dUqmUf3cqAFDVvXv3kSNHbtq0iXskCdTo6NGjT58+5a42VjZu3769b9++LSp9CQHcYs2aNauyspI76/nll19OmDDBx8enoKBg/fr1evwwcABdi4iIKCwsrLGXPhCRXC7/8ssvR44cqbzAiYgKCgrOnz+v3w9CrlELvQwJ5s2b5+rqum7duuPHj3OXcAwePHjevHncFSAA0DAuLi6//PLLnDlzJk2apHpDJeD8+uuv2dnZalclyWSyEydOCFWSgHAOGAAAQAA4BA0AACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAggJYYwMXFxeHh4Z07dzYzM5NIJN26dVu7dq1CoRC6LgAAaEEYlmWFrqFJFRYWBgYG3rlzx8PDIzg4uKKi4tixY+np6X/5y1/27t3LMIzQBQIAQItgJHQBTe3HH3+8c+fO2LFjd+zYYWRkREQFBQVBQUH79+/fs2fPuHHjhC4QAABahBZ3CPr3338nouXLl3PpS0QymWzhwoVEdOzYMSErAwCAlqTF7QEvXLgwMzOzbdu2qo1isZiISktLBSoKAABanBYXwBMmTKjeuHfvXiLy8/Nr8nIAAKCFanGdsKo7c+bMoEGDZDJZUlKSjY2N0OUAAECL0OL2gNVcuXLlvffeI6IdO3bwT18zM7OysjJd1gUAAMJosv3SFh3Ahw4d+utf/1pRUbFt27Z33nmH/4RlZWXN8MgBw+jJ8QysSLOiH2tBWJHmp3muSFNejNriekErff/99++9956hoeHhw4fHjx8vdDkAANCytMQ94KqqqmnTpm3btq19+/aHDx/29vYWuiIAAGhxWlwAsyw7fvz4yMjI3r17Hzp0yNbWVuiKAACgJWpxAbx8+fLIyMiuXbuePn2au/wXAACg6TXHc+C6k5eX5+TkVFRU5OnpaWdnpzbUz8/vp59+4jOfZtt3oBlW1QBYkWZFP9aCsCLNT/NckaasqmXtAV+/fr2oqIiIEhISEhIS1IYqb04JAACga83xB0jzhx9uOoUVaVb0Yy0IK9L8NM8VacqqWu5lSAAAAAJCAOuPxYsXC12CdmBFmhX9WAvCijQ/erMiDdYcjwA0f83zyAkAADQSDkEDAADoOQQwAACAABDAAAAAAkAAAwAACAABDABAjx8/3rlzp/LP2bNnMwxz9OhRAUvSij179vTq1Usmk1lZWY0YMSIuLk7oiuC/EMAA0NLdvXvXy8vr7NmzQheiZUuWLBk3blxGRsbUqVNHjBhx+vTpPn366N9qvrlwOU1D4DIkAH1y7ty54ODgqVOnbtq0iWspLCwsLi62tLQ0MTERtrYGu3//vre3d+fOna9evWpubk5EN27cCAwMdHBwSExMNDY2FrrAZgqXIQEACEkqlbZu3frNTV8iWrNmDcuyixYt4tKXiLp37z5x4sSUlJSoqChhawMOAhgAWrRJkyYFBwcT0ebNmxmGCQ8Pp2rngGfOnMkwTFJS0rx585ydnU1MTFxcXJYuXapQKK5cuTJw4ECZTGZnZzdq1Kjk5GTVmVdWVq5evbpr164SiUQqlQYEBERGRmqux8fHh6nd6dOnea7XxYsXiSgkJES1cfDgwUTEfyagU3j+DwC0aBMnTrSysoqIiPD39x8/fny3bt1qGzMkJMTCwuLvf/87y7LLly//6quv4uLioqOjp0yZMmHChNjY2PXr1z969Oju3bsikYiIKioqhg0bFhMT4+rqOmnSJJFIdOjQobFjx167dm3lypW1LSUwMLBt27a1DbW1teWzUizLJiYm2traWlhYqLa7uroSUWJiIp+ZgK4hgAGgRQsODmYYJiIiwsPDY+bMmRrGlEgkV69eNTU1JaJOnToFBwcfPHhw/fr13FQTJkzIyMjYv3//9evX/f39iWjZsmUxMTHvvvvunj17uKm+++67oUOHRkREhISEDBkypMal/Pzzz41fqaKiIrlcbmVlpdZuaWlJRHl5eXxmcj0888Y3mY0vpml0X2zXI1z9Ke/NHAIYAICXSZMmcTlKRH5+fkRkZGQ0depU5QhdunTZv39/SkqKv78/y7Lr169nGOa3335TTmVmZrZs2TJ/f/8NGzbUFsBaUVxcTETVe1pxZ7XLysr4zKRH+JsXaW8WBDAAAC8uLi7K/5uZmRGRo6Mjd7SZoxpvT58+zczMlMlk69atU51JRUUFEWm4HjcsLEztRLKqpUuXvv3223WWykU+tyxV5eXlRKTslgXCQgADAPAiFovVWlTTV01OTg4RFRQUfPPNN7UNrdGFCxfu3r1b29C5c+fWXSiRVCoViUTVDzXn5uYSkdqJYRAKAhgAQPukUikReXt7x8fH12vCO3fuNH7phoaG7u7uCQkJRUVFqvu7T548IaLOnTs3fhHQeLgMCQBaOoZhtD5PV1dXmUz24MGDtLQ01fbU1NR58+Zt375d60tUExQUxLJsTEyMauPJkyeJKDAwUNdLBz4QwADQ0nFHkvPz87U4T0NDw6lTp8rl8ilTpnBdooioqqpq1qxZq1at0nCWV1umTJnCMMyiRYuU63Xr1q1t27Z16NBh+PDhul468IFD0ADQ0jk5ORHR/v37g4KCxowZM2vWLK3M9h//+MfVq1ejo6M7deo0ZMgQiUQSFRX18OHDgICABQsWaGURGvj5+c2fP3/FihXe3t5jxozJz8/fvXu3QqHYsmWLkRG++ZsF7AEDQEvn5OS0YsUKe3v72NjYmzdvamu2YrH43LlzP/74o5WV1Y4dOzZv3mxmZrZixYoTJ05IJBJtLUWD5cuXb9q0ycbGZt26dUeOHAkODr506VJQUFATLBr4wEMFGgIPYwAA0Et4GAMAAICeQwADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAA5ckqBwAACAASURBVAAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAOAnnv8+PHOnTuVf86ePZthmKNHjwpYknY9fPhQIpEEBARUH7Rnz55evXrJZDIrK6sRI0bExcU1fXlQGwQwAOizu3fvenl5nT17VuhCdKWqqmr8+PElJSXVBy1ZsmTcuHEZGRlTp04dMWLE6dOn+/Tpo8cvxRsHAQwA+iw3N7e8vFy15bvvvktLSxs0aJBQJWnXkiVLatyvvX//fnh4uJeX17179yIiIrZt23bhwgXuIYkVFRVNXydUhwAGgJZFKpW2bt3axMRE6EK0IDY2dtmyZaNGjao+aM2aNSzLLlq0yNzcnGvp3r37xIkTU1JSoqKimrZMqBkCGAD01qRJk4KDg4lo8+bNDMOEh4dTtXPAM2fOZBgmKSlp3rx5zs7OJiYmLi4uS5cuVSgUV65cGThwoEwms7OzGzVqVHJysurMKysrV69e3bVrV4lEIpVKAwICIiMjNdfj4+PD1O706dP1Wrvi4uKPPvrI1dX1hx9+qD704sWLRBQSEqLaOHjwYCKq74JAR/BYZgDQWxMnTrSysoqIiPD39x8/fny3bt1qGzMkJMTCwuLvf/87y7LLly//6quv4uLioqOjp0yZMmHChNjY2PXr1z969Oju3bsikYiIKioqhg0bFhMT4+rqOmnSJJFIdOjQobFjx167dm3lypW1LSUwMLBt27a1DbW1ta3X2s2bNy8lJeXy5ctisVhtEMuyiYmJtra2FhYWqu2urq5ElJiYWK8FgY4ggAFAbwUHBzMMExER4eHhMXPmTA1jSiSSq1evmpqaElGnTp2Cg4MPHjy4fv16bqoJEyZkZGTs37//+vXr/v7+RLRs2bKYmJh33313z5493FTffffd0KFDIyIiQkJChgwZUuNSfv75Z22t2tGjRzdu3Lhw4cKePXu+ePFCbWhRUZFcLreyslJrt7S0JKK8vDw+iwi/ffKbOye1Um0TWOwTEu4bUvd4zQkCGACAJk2axOUoEfn5+RGRkZHR1KlTlSN06dJl//79KSkp/v7+LMuuX7+eYZjffvtNOZWZmdmyZcv8/f03bNhQWwBrS1ZW1rRp03x9fRcvXlzjCMXFxURkbGys1s6d+S4rK+OzlHDfNy/S3iwIYAAAcnFxUf7fzMyMiBwdHbmjzRzV6Hr69GlmZqZMJlu3bp3qTLjexRqutQ0LC1M7kaxq6dKlb7/9Np9qQ0ND8/LyYmJiVCtUxf0sqN7bmesQruyWBcJCAAMAUPXTqLVlGxHl5OQQUUFBwTfffFPb0BpduHDh7t27tQ2dO3du3YUSbd68+dChQ8uXL+/cuXNt40ilUpFIVP1Qc25uLhGpnRgGoaAXNABA/UilUiLy9vZma8Id/q3RnTt3apyEM3DgQD5L3717NxEtWLBA2X3aycmJiC5fvswwjJeXFxEZGhq6u7tnZmYWFRWpTvvkyRMi0pDc0JSwBwwA+oxhGK3P09XVVSaTPXjwIC0tzcHBQdmempq6cuVKX1/fCRMmaH2hSoGBgWq7sKWlpVFRUdbW1kFBQc7OzlxjUFDQgwcPuJ5iyjFPnjzJzUF35QF/CGAA0GfckeT8/HwtzpO7n1RERMSUKVP27t0rkUiIqKqqatasWUeOHKmtY5S2fP3112otL168cHJy8vDw2Lt3r7JxypQp69evX7RoUVBQUKtWrYjo1q1b27Zt69Chw/Dhw3VaIfCEAAYAfcYdnt2/f39QUNCYMWNmzZqlldn+4x//uHr1anR0dKdOnYYMGSKRSKKioh4+fBgQELBgwQKtLKKR/Pz85s+fv2LFCm9v7zFjxuTn5+/evVuhUGzZssXICN/8zQLOAQOAPnNyclqxYoW9vX1sbOzNmze1NVuxWHzu3Lkff/zRyspqx44dmzdvNjMzW7FixYkTJ7gd4uZg+fLlmzZtsrGxWbdu3ZEjR4KDgy9duhQUFCR0XfBvDMuyQtfw5mEYvG4AAHqoKb/esQcMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAgMCGDRvGMMyPP/5Y2wh37txhGMbd3V31Jk2bN2+2sLDIzMzks4gNGzbY29unp6droVzQEgQwAIDApk+fTkQ7duyobYStW7cSUWhoqPLpisnJyWFhYeHh4XZ2dnwWERoa2qZNmylTpmijXtAO3NO4IXAvaADQoqqqKicnp/T09Nu3b/v4+KgNraiocHR0LCgoePHihTJuhw4d+vDhw4cPH3LPW+QjJiZm4MCBkZGRH3zwgTar1y+4FzQAQAtiZGQ0efJkItq+fXv1oUeOHHn9+vWoUaOU6Xvu3Lno6OiwsDD+6UtEAwYM8PHxWbhwoVwu10rZ0EgIYAAA4U2bNo1hmN27d1dPR+74M3eYmrN8+XIjI6MJEyZwf16/fl0kEpmZmSUmJqpOxTCMj49PWVmZsnHy5MmPHz8+cOCADtcEeEMAAwAIz8XFpX///unp6adOnVJtT0tLi46OdnV1HTBgANeSl5d34sSJPn36WFlZcS09evRYvHhxWVnZlClTFAoFESUlJc2ZM0cikURGRpqamirnNnz4cCLas2dPE60VaGQkdAEAAEBEFBoaGhMTs3379iFDhigbd+zYIZfLuf1jruXMmTNyubxXr16q03755ZcnT568ePHi6tWrZ8+ePW7cuOLi4m3btr311luqo7m6utrY2Jw5c0ahUBgY1LED9sej4r1JxVpaOZ0b7S4Z01EidBX1xEL94XUDAK0rLy+3sbExMzMrKChQNnp6eopEovT0dGXLF198QUTbt29Xm/zp06etWrWSSCRcV+eJEyfWuJTg4GAievjwoQ7WQB805dc7DkEDADQLxsbGEyZMKC0t3bdvH9cSGxubkJDw7rvv2tvbK0fLyMggoupXH7Vr127dunXFxcVbtmzx8PBYu3ZtjUvhJuRmAsLSnwCePn26ubk5nzE3btzI1OTo0aO6LhIAQIPQ0FBSuSC4evcrIsrNzSUiiaSGw63Dhg2TyWRE5OfnV+MIRMR9T+bk5GizbmgQPTkHvGrVqo0bN9b2hlNz584dIurVq5fa+La2tjopDgCAHw8Pj4CAgHPnzqWmptrY2ERGRrq4uAwcOFB1HDMzMyLKy8urPvnMmTMLCgpsbGx27tw5atSo0aNHVx+Hy2+xWKybNYB6eOMDuLy8fP78+bUdbKnR7du3GYY5ceIE91MRAKD5CA0NvXTp0h9//OHk5JSfn//5558ru19xnJ2diSgrK0ttwu3bt0dGRgYGBm7YsMHX13fGjBm9e/du06aN2mjchO3atdPlSgAvb/Yh6KioKE9Pz7Vr1wYGBkqlUj6TKBSK+Ph4Nzc3pC8ANENjxoyxsLDYv3//77//rrxBh6quXbsS0d27d1UbU1JSuOuOuBPA3377bU5ODtcPS3U0hUJx7949mUzm6uqq6xWBOr3ZAbxly5a8vLyVK1eeOXPG2NiYzyRJSUnFxcXcOxgAoLkxMzP78MMPr1y5cvTo0ZEjR7Zu3VpthP79+xsZGZ0/f17ZIpfLP/roo4KCgu+//55L1k8//bRXr14xMTErV65UnTY+Pj4vLy8kJMTQ0LAJ1gU0e7MDePbs2SkpKfPmzeP/Zrp9+zYRWVhYTJ06tX379qampl5eXsuXL6+qqtJlpQAAfIWGhrIsW1ZWptb9imNjYzN8+PD4+PgXL15wLcuWLbt8+XK/fv1mzZrFtRgYGGzZssXExOTLL79U3VeOiooiookTJ+p+JYCHJrvgSdesra0lEkmdo33++efciru5uY0bN27AgAFcj4aQkJCKigqey9L8ki5evLhRawIAoNH169eJ6Ouvv67XVHK5vH379j4+PgqFQkeFvREWL17cTGKxxQXw1KlTTU1NlyxZonwLJiUlcQdtvv/+e57L0qcfLgDwJho5cqSjo2NZWRn/SQ4dOkREhw8f1l1VeqApv97157F6NjY2ZWVlRUVFdY7Jsqxar8KoqKjhw4d36tTp/v37fJaFxxECgLCSk5O5Rxspj+ppJpfL3377bVdXVy6GoTZ4HKFuqaUvEQUEBBDRkydPhCgHAKDeXFxcfvnll6VLl/K8p9Wvv/6anZ29adMmXRcG/OnPnhyfPWCFQvHo0aPCwsLu3burtmdnZ9va2lpYWHCXqNcJe8AAAHqpKb/e3/gbcdSLXC738fGpqKh49eqVauf+CxcuEJHa00UAAAB0R88PQbMsm52dnZ2dzf2iEYlEo0ePZll2/vz53FMziejZs2efffYZEX366adC1goAAC2J/hxKrfEQdF5enqWlJRFlZWXZ2NgQUXp6ekBAwJMnT7y9vQcOHJiVlXX48OGCgoLFixeHh4fzXBYOQQMA6KWm/HrXnyDhGcBc47Jlyw4cOPDs2TOJROLn5zdv3rzhw4fzXxYCGABALyGAmzsEMACAXsJlSAAAAHoOAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIoNYAzsrKWrFixYIFC/bu3atQKFQHnTlzZubMmbqvDQAAQG8xLMtWb01NTe3Ro0d6ejrDMCzL9u7d+8CBA/b29tzQVatWzZs3r8YJWwjuZRG6CgAA0LKm/HqveQ/473//u4mJya1bt0pKSnbt2vX48ePg4OCsrKymqQkAAEDv1RzA586dW7x4sa+vr6mp6bhx465cuZKfnz9y5MjS0tImrg8AAEAv1RzAeXl5Dg4Oyj/d3NyOHz9+//79SZMmNVFdAAAAeq3mAHZzc4uOjlZt6dKly65du/bu3btw4cImKQwAAECfGdXYOnPmzI8//rigoGDy5Ml9+/blGt95553Vq1fPmTOnTZs2TVghAACAHqo5gKdPn56Tk/P999/b2NgoA5iIZs+ebWpqGhYW1lTlAQAA6CdN/a3lcnl+fr6VlZVa+6tXr44cOTJjxgwd19Z84TIkAAC91JRf7wiShkAAAwDopab8eq/5ELSao0ePEtGgQYNMTEy4/2vwzjvvaKEuAAAAvcYr6hmGIaK0tLTWrVtz/9egJewaYg8YAEAvNbs94FmzZhGRRCJR/l9vTJ8+fdeuXUVFRUIXAgAALUuL3pPjbmotkUjqG8DYAwYA0EvC3wta75WXl8+ePXvevHlCFwIAAC1UQwK4srIyKiqK+395eflnn33m5+c3evToe/fuabU2XYmKivL09Fy7dm1gYKBUKhW6HAAAaInqva/9+vXroKCgP//8MzMz09bWdvbs2WvXruUGyWSyW7duubq66qBObRo9evSZM2cWLVoUFhZmb29fVlaGQ9AAAEDN/BD0999//+effw4dOpRhmJKSki1bttja2t64cSMyMrK0tHTZsmW6qFK7Zs+enZKSMm/ePENDQ6FrAQCAFqreUe/h4dGqVavY2FiGYY4cOTJy5Mj58+evWLGCiMaMGXPt2rXnz5/rplSdsLGxwR4wAABwmt1lSKqePXs2Z84c7mrg06dPE9GgQYO4QW5ubocOHdJufc2WhuuhFy9eHB4e3oS1AGjf/x3LrLE9crhdE1cCoF3h4eHffPON0FUQNSCATU1NKyoquP/HxMSIRCJ/f3/uz4yMDDMzM21W14xhDxj0XvWsrS2VAd4g4eHhGvaR6rzZlBbVO4A7dux48uTJ8vLyhISE+/fvBwcHm5ubE9GLFy/279/fpUsXHRQJAACgb+rdCWvatGkJCQmurq59+vQhounTpxPRxo0bu3btmp+fHxoaqv0aAQAA9E69Azg0NHTZsmX5+fmGhoZfffXV2LFjiej58+cFBQU//PDDhAkTdFAkAACAvml4dy+WZZXHylNTU6VSqYWFhfYKayLoBQ1Qo/87llnjOWB0wgL91kx7QbMse/fu3aysLCcnJw8PD9Uz1U5OTjqoDQAAQG/xPQQdHR3t4uLi6+sbEhLi6enp6+t769YtnVYGAACgx3gFcFxc3MiRI58+fWppaent7S0Wi+/cudO/f/+UlBRd16dr2dnZeBYhAAA0PV4BvHLlysrKyjlz5qSlpcXHx2dmZo4dOzY/Pz8iIkLX9QEAAOglXgF8+fJlOzu7n376ycTEhIgkEslvv/1mamp69uxZHZcHAACgn3gFcEZGhqurq0gkUrZIpVIPD49nz57prDAAAAB9xiuAy8vLjY2N1RplMhnOngIAADRMvW/EoYRrYQEAABqs4QEMAAAADYYABgAAEADfO2Hdv39/9OjRai1EpNZIRHv37tVKZQAAAHqMbwBnZ2fv27evenuNjQAAAKAZrwDesWOHrusAAABoUXgF8Pjx43VdBwAAQIuCTlgAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAgAAQwAACAABDAAAIAAEMAAAAACQAADAAAIAAEMAAAgAAQwAACAABDAAAAAAkAAAwAACAABDAAAIAAEMAAAgAAQwAAAAAJAAAMAAAjASOgCAJq7rNtlt3/Iykkot/U19fvS1qKjidAVAYA+wB4wgCZJu/OPDn3qECDuv7mNpafJ/r4pqSeLhC4KAPQB9oABavU8uujS/LR3z3Sw6mRCRHbdzBz7SqLee/7O0XZ23c2Erg4A3mzYAwaoWXme/Gzoy8F7nLj05bTuIw5c43DqoxfyMlbA2gBADyCAAWp2Y0lm+3ekjoEStXa3Ma2sOpncXZUtSFUAoDcQwAA1KEipeLgjr8cSuxqH9v6h9e2fssvz5E1cFQDoEwQwQA1uLs3y/tjazLbmThIW7sYdRsju/fK6iasCAH2CAAZQV5JelXygoMsn1hrG6fqFTfyaHJwJBoAGQwADqPvz1xy3/2tlamWoYRyLjiZ2fmZJe/KarCoA0DMIYID/wSoocWuu10yrOsf0+n9WDzblNkFJAKCXEMAA/yPtYrGxhaF1F9M6x3Qeap7/uKIgpaIJqgIA/YMABvgfj3bnd/zQgs+YBkaMy/uypN35ui4JAPQSAhjgvxSV7JN9+e4fyHiO3/GvFg//hdPAANAQCGCA/0o9VWT5lom0vTHP8R36iKuKFTl/lum0KgDQSwhggP9K3l/gOqZVPSZgyHVMqyf7CnRWEQDoLQQwwL+xcjblSGGHkdJ6TeUySpZ8EAEMAPWGAAb4t/SrpeZtjGQd+B5/5rTuIy5Jqyp8VqmjqgBAXyGAAf7t2fHCdsPrt/tLRIwBOYWYPz9RqIuSAECPIYAB/u15dJHz4HoHMBE5hZinnizSej0AoN8QwABERCUZVQUpFfa9zBowrXOI+YszxYoq3BcaAOoBAQxARJR6qqhtf4mBEdOAac3sjGQdjDOulWq9KgDQYwhgACKi1BMNPP7McR5s/vwEjkIDQD0ggAGIWEo9VeQUYt7gGTgPNk9FPywAqA8EMABlx5cZtzKUthM1eA72vcW5ieVlOXItVgUA+g0BDEAvThe1HShpzBwMjRmHAMnLM8XaKgkA9B4CGIBexBQ7DWj48WeO0yDz1FM4DQwAfCGAoaVj5Wx6bIljv0btARNR2wGSl2cRwADAFwIYWrqs22XmbUWm1oaNnI+1l2lFgaLoBe5JCQC8IIChpUu7WOzQV6yFGTHk0Ff88ixOAwMALwhgaOleXShxDGzs8WdOm36SVxcQwADACwIYWjaWXl0sduyrnQB2CJCkXSrRyqwAQO8hgKFFe32vzNTGSOJopJW5WXcxLUmvKs2s0srcAEC/vfEBvGfPnl69eslkMisrqxEjRsTFxdU5ycaNG5maHD16tAkKhmbl1YVix0BtnAAmIiLGgFr3Fqddxk4wANRNOz/8hbJkyZLFixe3b99+6tSpOTk5v//++4kTJ06cOBEcHKxhqjt37hBRr169JJL/OfBoa2ur23Kh+Xl5rtjlLzItztAxUPzqYrHLe9qcJwDopTc4gO/fvx8eHu7l5XX16lVzc3Mimj17dmBg4NSpUxMTE42NjWub8Pbt2wzDnDhxQibDt2TLxtKriyUBEQ5anKVjP8mF2WlanCEA6Ks3+BD0mjVrWJZdtGgRl75E1L1794kTJ6akpERFRdU2lUKhiI+Pd3NzQ/pCTkK5sdTA3Knht4Cuzq6bWf7j8vI83BQaAOrwBgfwxYsXiSgkJES1cfDgwUR0+vTp2qZKSkoqLi7u2rWrrsuD5u/VhWJtXYCkZCBi7HuI03EaGADq8qYGMMuyiYmJtra2FhYWqu2urq5ElJiYWNuEt2/fJiILC4upU6e2b9/e1NTUy8tr+fLlVVXoudripF0s0WIPLCWHvuJXuBgJAOrypgZwUVGRXC63srJSa7e0tCSivLy82ibkemBt2LDhwoULffr0CQgISE5O/tvf/jZ8+PDKStxEsGV5dbHYQUtXAKty7IvbcQBA3d7UAC4uLiai6j2tTExMiKisrKy2CbOzs01NTZcsWfLo0aNdu3adPn06Pj7e1dX15MmTK1eu5F9AjRcyccLDwxuyStC0Cp9WKKrYVq61dtZrMPteZq/jy6pKFVqfMwA0Xnh4uIYv8Kas5E0NYFNTUyKqqKhQay8vLyciZbes6jZt2lRSUrJo0SLlC+3m5vbzzz8T0fbt2/kXwNYOAfxGSL9a6tBH+8eficjIzMCqk0lmXKkuZg4AjRQeHq7hC7wpK3lTA1gqlYpEouqHmnNzc4lI7cSwmuq/cQICAojoyZMnWq0RmrX0qyX2vXQSwETUuo84/QpOAwOAJm9qABsaGrq7u2dmZhYV/c8TWLkQ7dy5c41TKRSKxMTEGzduqLVze9JmZma6KRaao4xrJa11F8C9xelXsQcMAJq8qQFMREFBQSzLxsTEqDaePHmSiAIDA2ucRC6X+/j49OzZMz09XbX9woULRNSrVy+dFQvNS1WJIud+ua2fqY7m/+894CY9mgUAb5g3OICnTJnCMMyiRYvy8/O5llu3bm3btq1Dhw7Dhw/nWliWzc7Ozs7O5o7si0Si0aNHsyw7f/58heLffWSePXv22WefEdGnn34qxHqAANJjS63fNjUy09X737ytyEhikPuwXEfzBwA98AbfitLPz2/+/PkrVqzw9vYeM2ZMfn7+7t27FQrFli1bjIz+vV75+fncHZ6zsrJsbGyIaMWKFbGxsbt27bp3797AgQOzsrIOHz5cUFCwePHiQYMGCbk+0ITSLmntEYS1cQwQp10qsfQw0elSAODN9QbvARPR8uXLN23aZGNjs27duiNHjgQHB1+6dCkoKEjDJK1bt46Li1uwYEFpaemaNWuOHj3avXv3o0ePoutyi5J2qcTBX1cngDkOfcVpl3A1MADUimniXtf6gWHwur3BFFXsZuuEj1LeMrUy1N1SchPKj414Nv5xx0bOp7Cy/PeUu+3MLQc6umulMCJKzM+ceHFPUGvXZX5DDZmaf4X/37HMyOF2fBoB9ElTfr2/2XvAAA2QfafM3NlYp+lLRJYeJhWFiqIXjbq9WnppYdfDEUdS74ddO/jJtYNaKYwl9q/nd45y9jqf/mT745tamScANAACGFqctIvFjn11e/yZiIghB39xWiOeyqBg2ffPbJvk1u3ggMnX3gmLepEY9SKh8XUdS00goi+6BEf0eHfJnVMKHMsBEAgCGFqctCs6PwHMcQgQp11s+Gng9YlXjA0Mv3x7ABFJRSYruo/46lY02+hrm3Y8ufmxRx+GmN527SRGxjeyUxs5QwBoGAQwtDjpV0tb6+YmlGoaczuOMnnV0viY1T1HMfTvG7eNdO6kYNnTr5IaU1JxVcXJl4/+0s6b+3OEc6cjqQ8aM0MAaDAEMLQsRamVikpW1kH7z2CozrarWd7D8sqihjyV4Z9JN7pZt+1i5aBsYYgJ7dhza5L6fdzq5eTLRz1tna1M/v37Y1hbz+gXtT67EwB0CgEMLUt6bEnr3k2x+0tEhiaMlbdp1s167wQrWHbVgwsLvIPV2se5+B5/kVhQWevDvuoUk5Y0yPG/HbN72Dgl5mcWVeKGIQACQABDy9IEVwCrcvAXv7pU735YUS8SZCLTvvYd1NqtTMS97dqdePmwwfXEvErq7+Cm/NPE0MjHyvE6TgMDCAEBDC1L2sVix0Dd3gNLlWPfhvTDWpd4Zbanf42DRjh1OvK8gWdt00oKssqK3rZyVG3sY9/+ckZKw2YIAI2BAIYWpKJAkfe4wsZXV89gqM6hryT9aomiqh5dl5MLX9/MfvFBB58ah77j1Cn65cOG9YW+kZ3aw9bZ4H8fx9nTxjnu9YsGzA0AGgkBDC1I2uVi+x5iQ2P1B0LrjqmVobmT6HV8Pc7a/vYwdoJbN1PDmu/T7iSxsDA2/TM3vcahmsW9ftHDxlmt0ce6zd2cVw2YGwA0EgIYWpBXF0qa4hYc/8sxUPLqIt/TwBUK+T8fx01/S9OTMYNau55Lf9KASuKyU7tat1Fr7GBulVNekleBpxcDNDUEMLQgTXwCmOPYV8L/NPDh5/c7Wdi7y2w0jBPk4HY+PbkBlcRlv+hm46TWaMAw3pYO2AkGaHoIYGgpqkoV2XfLn2zDmwAAIABJREFU7HuaNfFyHfqKX10s4XnSdtOja9M69tQ8Tl/7DhfTk+t7Gvh5cZ4hwziKZdUHvW3lGJ+bVq+5AUDjIYChpciMK7XqZGIkbur3vHlbkcjcIDex7mttHxVk3cl5pbxNVW2cJBamhkZJBdn1KiMuO7X67i/nrVa2D/Oz6jU3AGg8BDC0FGmXShwCmvr4M6dNP8nL83UfhV59/+LMt3rX1v1KVWBr14vp9bt26Pbrl77VTgBzOspsHyGAAZocAhhaivSrTXcPLDWO/SSv6grgrLKiPSl3Pvbsw2eG/nbtr2Q+rVcNt16/9LNuW+Mgd5lNffenAaDxEMDQMrCUfrWkaZ7BUJ1joPjVhToC+JeEy2Pav21nas5nhr3s2sVmPatXDbdzat0Dbm9ulV5aUCavqtcMAaCREMDQIuTcLzO1NJQ41n10VxdkHYwZA6YgpaK2EXIrStcnXvmii/rNn2vjZdE6tTgvv4Lv5cWZZUXl8ipniUWNQ40MDNqbWz0pxE4wQJNCAEOL8PJ8iWM/YU4Ac1r3FqddrvVq4H/cOT26fZf25lY852ZkYNDVuu317Oc8x7/z+pXaHSjVuOEoNECTQwBDi/DqfLHAAewvTq8lgG+/fvmvJzeX+A6u1wx72Dpdz+IdwDkvfeoIYOvHCGCApoUAhhaApVcXBLgFR4VCnlz4OrU4jyXWMVBS42ng3IrSsef+tbrXKFt+Z3+V/Kzb3uR9D+f43LQ69oClNk8KX9erAABoJGFOiQE0pZz7ZSKpgdRZ1GRLfF6ct/j2iQPP7lmZiCvk8oLKsr52Lo6O9gNftbF1/G9HsKyyopGnt45w7jS2lkcvaOBn3fZvN47yHPnP3PRPO/fTMIKbzObg8z/rWwMANAYCGPTfi7PFbYKbbvf34PM/Z1ze+7Fnn6T3v+D2a3PKS06+erRywCXXk0tD2r0VYNfBysQsIT9zy6PrMzx6f+Mb0oCluMqsCyrLs8qK6tx1rlIokgqyPC3sNIzjJrN5XIA9YIAmhQAG/ffybLHbmFZNs6yNj659e+fU0UFTu6vcdsrKRDy2g4+nmdOzuPz8vrnXsp4XvC5rZ255buj/87Swb9iCGGK4o9BD2nhoHvNRQVZbsYWZoaYDAO0klumlBeXyKhMetwEBAK3Ahw30HCtnX50v7rdO0xlQbYl+mRh++8SFYbNcpdbVh7btL7m35vVHa/w+cvXTyuK62zrdyEqtM4D/zE33smyteRwjA4O2EounRblvtbLVSm0AUCd0wgI9l3WrTOwoErfW+W/Np0U5ky5G/hE8ocb0JSJrL9PKIkXh01qvBq6v7jZON7JT6xztbs6rLlYOdY7mIrVKKcJRaICmgwAGPffibHHb/jo/AVypkI89968vvPv3sWtf60gMtQmWvDzH99GEdeIZwHdyXvlY1XwPLFUdzK2TC3O0URcA8IIABj336lxT9MD6Lv6Mjankk84Bmkdr00/y8nytt+OoLyeJBRGlFudpHk3DTShVuUitUhDAAE0IAQz6TFHFpl0pceyr2wBOyMtYk3B5Q5/RDDGax2wTJHl5pkiLi+5p63xN4+04MkoLNdyEUlUHqVUyLgUGaEIIYNBnWTfLZO1FptaGulsES+zHV/cv9hnURlx3R2uLt0zklayGm0LXV4+6Avh2zktfHsefichFap1ShD1ggKaDAAZ9lnqy0GlQ/e4wVV//TIorkVfO9OjNc3yngeapp7S2E9zT1vmaxsci3cx+0c2m5qcQqulgjj1ggCaFAAZ99vxkkVOIDgM4q6zo7zejfusz2pDh+1FyGqTNAO5u43Qn51WVQlHbCDdfv/DjF8BWJmKWJf5PWAKARkIAg96qKFC8ji9zCNDhM4D/Fndsgpuf5tssq3EKMX95plhRxWqlAJnI1FlieS83rbYR4rJfdLN2qm2oGmdzi2dFuVopDADqhAAGvZV6qsjBX2Jkpqs3+fWs5ydfPlz09qB6TSW2N5K2E2VcK9VWGRr6YaWXFpZUVbSXWvKcVTtzy+fFCGCAJoIABr31/Hih81BdHX+uVMhDL/+xvPsIqcikvtM6D5E+jy7UViU9bZ2v1fJg4KuZz3ratquzb/Z/C5NYYg8YoMkggEFPsfQsuqjdUKmOZv9d/Blnc8u/uvg2YFrnoebPjmvtNHCAfYcL6ck1DrqS+bSPXTv+s3KWWDyv66piANAWBDDop9f3ykRippWbsS5mfj8vfU3C5V/7vN+wyVv3FhckV5RmVmmlGE8Lu5Kqihr3XK9mPu2t4c5c1bQzt3yOPWCApoIABv30LLrIabBOdn+rFIppl/74tutgPhf+1sjAiGkTLNFWX2iGmKDWrmfTHqu1l8or7+am9bDh2wOLiNqZWz5FAAM0FQQw6Kfn0YXOg3VyAnhp/GmpyGT6W70aM5N2Q6RaPAo9wNH95KtHao1XM595W7Y2r88p6vbmVk9xLw6ApoIABj1U9lqefbus7QDt34HyQnryhoex/+w7ln/Pphq1Hyl9frxQXqGdi5GGtfU88fJhpUKu2ngu/Ulwa7d6zcdBLM2vKCuVV2qlKgDQDAEMeujpkYK2A7V/AVJGaeGHF3ZuDfg/R7GskbMS2xtZvGXy6rx2nozkKJa5Sq0vZaSoNp58+bC/Y/0CmCHG2dziKR7JANAkEMCgh5IPFLiMamxGqpGzinHnd0517zm4zVtamaHLe7Lk/QVamRURvd++y+7k28o/UwpzUgpz+tm71nc+7c2tcBoYoGkggEHfVBYqXl0oaT9CywE8O/aAiYHR1z71u+2GBq6jZckHCli5do5Cf+Tqt+/ZvZKqfx893p1ye3T7LkYG9f6Ad5Ba4ZEMAE0DAQz6JuVwgWOg2Fimzff2P+6evpGV+kfwRwZMo079qpJ1MJa0Eb26qJ3HAzuKZQH2HbYmXSeiUnnl2oTLoQ3qJoZ+WABNBgEM+ubxHwWuoxt4gVCNtiRd35p049igqfXqUcyH2weypD352prb1z6Dvos/k1ZSsOxuTE9bZ5/63KFaqYO5VTLOAQM0CQQw6JXKIsXLs0Xt39HaFcB7n8YvuhV9PGSavZn2ryp2H2fxZG++tvpC+1m3ndMpwHXvdzuTb63r3cCbhLhIrVLwUEKAJmEkdAEA2pRyqMChr8TE0lArc9v3ND7s2sHoQaEdZbZamaEaqbPIqpNp6kmt/WL43Dt4WseeEiNjU8MGfrRdpNbYAwZoGtgDBr2SuC3PY4KFVmYVmXJnduyBqEHTulg5aGWGNXIf1yppt9aOQhORtYm4welL3FOBic0p186ZaQDQAAEM+qMotTLrVmmHkVro//zzg0uf3ThyesiMhp1J5c91tOxZVGFloUKnS6kXF6l1Mo5CA+geAhj0x4NNuR3/amFo2qiOyiyxf4s7+uvDq5eGze5s0VpbtdXGzNaoTX/Jo13N6BlELlJrXIkE0AQQwKAnFFXsg825nafzffh8jcrlVRMv7Lmc8fTy8NntzBs1K/68Zlr9+WszCjwXqdWTAuwBA+gcAhj0RPL+glZuxlZepg2eQ3ppYf/oX8vkVacGz7A0NtNibZo5DTSXl7Fpl5rLaVc3qc3jwmyhqwDQfwhg0BP31uR0mW3d4MnjslN7HFk9pM1bkcHjxUYiLRZWN4a8Pra6t6a57HS6y2ySChDAADqHAAZ9kHW7rPB5RYdRDbyYZ0vS9WGnNv/ca9Qin0GNfMxRw3hOskw9XVT0olk8hshNZvMYAQygewhg0Ad3I7K9PrY2MKp3dpbLqz6+un/Fn+cvDPt4lLOXLmrjQyQ1eGu8xb01zeJMcBtxq7yKUqGrANB/CGB44xUkVzyLKmxA96vnxXkBUWszy4pi35nj0cpOF7Xx5zPf5sHGnLJsed2j6pgBw7hIG34wHwB4QgDDG+/Gt1nes6xNLOp396tLGSm9j/48zsXnj+CPZKKGd93SFnMnkcv7sjsRzeLYr7vMRugSAPQfAhjebDn3y58fL/SZX7/AWP3g4uiz27cGjP20cz9BTvrWqPvXdvc35JSkVQldCOno1psAoAoBDG+2q1+kd/3Clv/DB8v/f3v3HtbUmecB/JeE3ENuXAoidwQVrCJe8AZYt45Ta7XdOnYere6Udkq3dqpOfR7r1Grr7K6z+ridVp+ZrmNbq9VinXopsCJax1rForQzVgSlVAFR7gm5kBtJ9o/T5ckGVFTIyUm+nz98zJs38Pvx5ry/nJP3nOPsee7rwo9qz597/JXZMalDGtu9UgwXjnpOU7Ghle1AKI3tA/IAwQAFGDis8ZhJV2PL+FftAPtfN3VOL9lucti/nrs8QTHQV/lS1tqIa4cN7X+3shtGqgp7wABDDgUYuKrH4vpq+c2c94YJRAM6hnykoSq76L3FSeMLZy6Rh4iGOrz7I1YLJv/bQ38raHI7B+cehfcnDYegAYYeCjBw1bnXWyInyuLmKO7a0+CwvlT+11e/OXxo1r+sSJ/hP1/69mv0cxqhnP/dFjZXY4VL5ETUbDGyGANAwEMBDhwbNmxgO4TBMZBEGo+b6j435Gy7+40Cixovjzm4pcfl+vv8VdkR8YMQ34Dd54jw6JEPh//jvzpavmH5ZNwL7Y0UZO8rTkAiAYPndrN5pIujeDx//Lv5Z1T34a6JmBodB7Lr/ml37PBH5Hfo9qOxY1XFkWp965+n/vPM6JTBDvPuHmREfjxk+Hpl8y8qkyXaezu9arAsKm5Njf524/g5nlksKm4tnMvV9VnBs4FwhX8m4suosAcMHGPTO4vm1meujrhD9W22GF87/8WkL/44JTLh4oLfslJ9H1DSAmXKQuXRpxucdtZmqHNt9Wz9aoBgEMJ2AAD3wGFyFc2tj52tGLui/0s13ew2/PHy6Z1Xv1mcPP77Ba9Fy5Q+jnAQTdkUdXRhQ+kvGn+2P3aAC80G14X2GzYn+ycl35Wb3Hq71elyGRw/rR63Onsszp8urK0SSvg8HhFRmFJntwh5fIVQzFaoAJ6CtAB/+umn77zzzuXLl0NCQqZNm7Z+/foJEyawHRTchaWtp/jx+ogs6bTNUX2frWhr2FZ9pvhG9eKk8X+f/9vhcpXvIxxcPD7N3hdb+ovGo083zN4XK5T7+njVOO2wo01XfPxLvbRZTR227kaz/la3odVqarOaO2zmW93GLruly2Htslu77FaDw6oRSfk8nur/biIpFoTIBD/d0kpvt/x0DGHNouTP/t3uclqdPSqRRC2SqkUSjUg2TKYMk8iipcpomVIjksYrNNFSpVYs+6lsAwyZYCzAb7/99vr16xMSEvLz8zs7O/fv319aWlpaWjpz5ky2Q4Pbaj7bXbbkRuoS9eS3Ij1XMd8wdx24fvGD2gqr0/FCavZ72U+qROxfV3KwCES8n+2PPVVw8/NpP/5sf6w61ae7bs8kjdv343c++EU9LleDWVetb6036xrN+hvmrnarucGsu27qlAqE4RL5cJl6mEwZIVFEShWpyoiHpAq1SKoSSVQiqUooGeCI83i8Treb+XVdDmuX3aKzWzpt3be6De02c7PFeEnX3Gnvvm7UNVsMph57jEwVJ1fHKTSJCm2sXB2nUMcrNHFyjUQQjNMmDAV//A58SFVVVY0ZMyY9Pb28vFyhUBDR+fPnc3JyoqOja2pqRKIBnR6KtQNDyisRh8n17X+2X97Rmfd+TOIToURkctjK2+r/dqvuf5pqGs36ebGjl6ZMyI1K8rfziwZxRC79ubPizdbxa8LHLA/zzeHoRcWt7z8aOuKvm9pXvOtu0fU2PsgiLKfbdbPbUK1vvWbqaDR31Zt0102dDSZdi9UUJQ0dpXooIVQTK1fHytVhYlm8QhMv1wzi4eJ7Gg67y9lg0jWY9Q0m3XWTrsGsazDp6826JnNXjFwVL9fEKtQJCk28XJMYqk1TRUZJ7/NWmPchULd0P+HLqILuo9y2bdvcbve6deuY6ktEEydOXLZs2fvvv19SUrJgwQJ2wwNPltaeqv/WXXy3QztXGFMmKBFc/O5M07nWhnqzbnxYTG5U8ruTF0yJjBfwAn8tYUaBdvgsxZlVt77f1jlmuXbUc5p7vfnEfVCLpCtGz3hjSVWPyxXCv4c/covF2GnrbjDrG836BrP+urGz3qxrMOlvWQxhYtko1UNJodo4hWZ2TGq8XBOnUMfIVEI+O4u9b0fEF6Qow1P63JTC4XJeM3U2mPQNZl29SXeyuW7H1W9+MLQbHNY4uSZeoYlTqOPkmgSFJk6hiZOrh8lU2GOG2/HHDyBDKiMjo6qqSqfTqdXq3saDBw8+9dRTL7/88rZt2wbyQ/DBbegYHbbExElvv7Hj+0st1/W6rvGmxqiOLpc1QxP1sGbYWG30lMiE0eqHRH42X/drKEakpcLy/Xsd174wxuTJE+aFDp8pVyYNyVW9mJ3dHpdLuOrp2Y/N2Txh3sPaaKax09bdaetmDtvqbZYWq6nNarrZbWg06zts3deNnQqhWCuWxsk1cQp1rFydoNDGKzRxcnWsXM1ioR3SDcTidNSbdPUmXaNZz/yn9zOHmB8Sr9DEKzQakTRapoyWKsMksmEypUYki5QoIqWKe30zB8aWTv6aiC+j8sf8h47b7RYKhVqttrX1/13v/uLFi2PHjp01a9bx48cH8nPwvhkIu8tp7rETkaXHYXU6iMjksDvcToPd2uWw6rotnQZLp9F8Q2/osJrbus2tdnN7iMlBPap2SYIkPCUqbGRaREZE1FhtdGKo1t8OLw/E0I2ITe+sLzFe/8J483R3j9kVNlYSliFWp4pD44VirUAeLRRrBDw+idWC+/6z9R5t5gkE//Fd2Y6r3zRbDI9rVx7u2CoVCMPEsnCJPFqqVIkkkVJFpEQRLVXGKTRhYlmCQisLEQ5mtoOErQ1Eb7fUm3QNZn2nrbvZYrzVbeiwdTd1d+ntllaLqcVqFPEFKpFUKRSrRFK1SKISSjViqUokUQklKpFUxBcohGIhnx8qFIfwBEqhePKkSbX/uMSsERPw+EqRhIgkghCpwB//7Hfgb1MWw6dRuYOJwWAgorS0NK/2hoYGIsrKyhrID3H53zsGAAAGyxAUn/4F15cTZrOZiPqutBKLxURktQ7oFjQ8opra9sZm/aCH558k/BAJb0Dvk1C+hM/jE5FaIBGpBX1P4uDxSaQS8AS8gd89EB6QTe8k90//ej+lc7IREQwtk8vW8387CT1up9llu6eXO9wus8s+BHFxRkqc7+6TFlwFWCKREJHd7v32stlsRNS7LOuu0lLC0lL6vxAEgF9h1mqJNRz4yhwGRQRJ2Q4BBiq4dkRCQ0OFQqFe773zqtPpiMhzWRYAAMCQCq4CLBAIRowY0draajKZPNvr6uqIKD09naW4AAAg6ARXASaivLw8t9t94sQJz8Zjx44RUU5ODktBAQBA0PHHVeBDqrKycuLEiRkZGadPn1apVET07bffTp8+PSoq6urVqyEhwfWlOAAAsCXoCjARrV69esuWLbGxsQsXLuzq6tq3b5/T6Tx69GheXh7boQEAQLAIxgJMRDt37ty+fXt1dbVSqZw4ceKGDRtwNyQAAPClIC3AAAAA7Aq6RVgAAAD+AAUYAACABSjAAAAALEABBgAAYAEKMAAAAAtQgAEAAFiAAgwAAMACFGAAAAAWoAADAACwAAWYw8rLy5944onw8HC5XP7www9v2bLFZrN59fn000+zs7OVSqVWq503b96FCxdYCXWArly5IpfLp0+f3vcp/0/EbDZv2LAhPT1dKpXK5fIJEyZs377d5XJ5dfP/RPriYswBORyc3kACb74aBG7gpj179ggEgtDQ0Pz8/FdffXXEiBFENH/+fJfL1dvnrbfeIqKEhIQVK1YsXbpUIpEIhcIvv/ySxbDvwOFwMFfknjZtmtdT/p+IwWAYN24cEY0cOfKll17Kz8+Piooioqeeeoq7I8LgYswBORyc3kACb74aFCjAnHT9+nW5XB4TE3P16lWmxWq1zpgxg4iKi4uZlkuXLvF4vIyMDKPRyLRUVFRIJJLExESbzcZO3He0bt065kOh1/zCiUTeeOMNInrmmWccDgfT0tXVlZmZSUR79+5lWjiRiBcuxuwO0OHg7gYSkPPVoEAB5qTVq1cT0WeffebZeOrUqV/96ldFRUXMw4KCAiIqLCz07PPiiy8S0cGDB30X68CUl5cLBIIFCxb0nV84kUhqaioRNTY2ejYeOHCAiBYvXsw85EQiXrgYszsQh4PTG0jgzVeDBQWYk9LS0hQKhd1uv0Of9PR0ItLpdJ6Nn3/+ORG9/PLLQxzgvTGZTCkpKampqVeuXOk7v3AikV27dm3evNmrsaSkhDnsyTzkRCJeuBizO+CGg+sbSIDNV4MIi7C4p7u7u7a2NiMjw2w2r1ixIi4uTiKRjB49+p133uldY+J2u2tqaiIiItRqtedrk5OTiaimpoaFuG9v5cqV165d+/jjj2UymddTXElk6dKlr732mlcjs8uVlZVF3EnEExdjZgTYcHB6Awm8+WoQhbAdANyz5uZml8vF4/Gys7ObmppycnKI6Kuvvlq5cuX58+c/+eQTIjKZTE6nU6vVer1Wo9EQkV6v933Yt1NUVLRjx47f/e53kydPvnHjhtezHErEy5dffvnRRx+p1epf//rXxM1EuBjz7XB3OLi+gQTYfDW4sAfMPUajkYjKy8v5fH5VVVVxcXFxcfHFixeTkpL27t1bWFhIRGazmYhEIpHXa8ViMRFZrVafR92/tra2559/PjMzc/369f124EoiXs6ePfvkk08S0e7du8PDw4mbiXAx5n5xdzgCYAMJpPlq0KEAc49AIGD+s3379ri4OOb/iYmJf/jDH4hoz549RCSRSIjIbrd7vZY58U6hUPgs2jt74YUX9Hr97t27hUJhvx24koinw4cPP/roo93d3bt27Xr88ceZRi4mwsWY++L0cATABhJI89WgwyFov1ZWVrZy5UrPljVr1uTm5hIRn8+fOnWq51PMwx9++IGIQkNDhUJh30M3Op2OiLy+aPGBfhOx2WyHDx/evHkzs/6iX5xIZMmSJb0PN23atHbtWoVCceTIkZ///Oe97f6WyEBwMWYvnB6OnTt3cm4D6UulUhHX5iufQQH2a11dXVVVVZ4t7e3tw4cPDw0NNZlMXp2ZFQ3MMg2BQDBixIjq6mqTyeT5+bGuro6I7rA9D5F+EykqKiKi1atXM2cp9Dpz5gyPx0tPT7906RInEmH+09PT8/zzz+/atSshIeHIkSNjxozx7OZviQwEF2PuFQDDsW/fPuLaBtIXF+crn0EB9mtPP/202+3u256bm1tUVFRWVtZ7SI2IKioqiIi52gAR5eXlXb58+cSJE/Pnz+/tc+zYMSJi1kH4Ur+JGAwGr8+2FoulpKQkLCwsLy+v92iV/ydCRG63e8mSJYWFhVOmTDl8+HBERETfPn6VyABxMWYKlOHIycnh3AbSF4/H49x85TvsnP0ED4Z5X6alpbW0tDAtra2tI0eO5PF4586dY1ouXLjA4/HGjBmj1+uZlsrKSqlUmpiY2Ht5IH/T2NhIfU5z5EQizBda48ePN5vNt+vDiUS8cDFmd+AOB0c3kECdrx4cCjBXvf7660QUFhZWUFBQUFDAXOp27dq1nn2YUyFjY2NXrVqVn58vk8nEYvHJkydZCvnu+p1f3H6fiE6nY46bjRo1KrePVatW9fb080T6xbmYA3g4OLqBuAN0vnpwKMAcdujQoby8vNDQUIVCkZ2d3XuRW09/+ctfMjMzJRJJZGTk3Llzz58/7/s4B+5284vbvxMpLS29w0GmWbNmeXb250Ruh1sxB/BwcHQDYQTefPXgeO7+vtACAACAIYXzgAEAAFiAAgwAAMACFGAAAAAWoAADAACwAAUYAACABSjAAAAALEABBgAAYAEKMAAAAAtQgAEAAFiAAgwAAMACFGAAAAAWoAADAACwAAUYAACABSjAAAAALEABBgAAYAEKMAAAAAtQgAEAAFiAAgwAAMACFGAAAAAWoAADAACwAAUYAACABSjAAHAXs2fPnjZtGttRAAQaFGAAuBOTyXTq1Km5c+eyHQhAoEEBBoA7KSsrs9vtKMAAgw4FGADupLi4OCYmZuzYsWwHAhBoUIABglFFRYVQKJRKpTU1Nb2NH374IY/HGzdunNVqZVrcbndJSYnn7m9DQ8NvfvObUaNGyeVyiUSSmJhYUFBw8+ZNXycAwH08t9vNdgwAwILf//7369atmzJlytdff83n82trazMzM4mosrIyLS2N6VNZWTlhwoQjR47MmzePiGpra6dMmaLT6ebMmTNq1Cij0VhWVnbt2rXk5OSqqiqxWMxmPgCc4waAoOR0OmfMmEFEW7dutdvtWVlZRLRr1y7PPhs2bBCLxWazmXn4y1/+koj+9Kc/9XZwOBzMC48fP+7T6AG4L4Td8g8AbOHz+bt37x47duy6desuXbpUWVm5bNmypUuXevYpLi6eOXOmTCZjHi5dunTSpEn5+fm9HUJCQvLy8iorK9va2nwaPQD34RA0QFDbu3fv4sWLiWjkyJEXLlyQy+W9T7W0tERHR7/77rvLly/3fInFYqmqqqqrq6urq7t48eLx48c7Ojr27NnD/BwAGCDsAQMEtccee0ypVBoMhqysLM/qS0QlJSVut9tzBVZnZ+eaNWv27NljsViISCqVZmVlJSYmdnR04KM8wL3CKmiAoFZQUGAwGMLDwz/55JMDBw54PlVUVDR69OjExMTelnnz5u3YsSM3N/fQoUO1tbVGo/H06dPTp0/3edQAgQAFGCB4ffzxx4WFhTk5OadPn5ZIJC+++GJTUxPzlMPhKCsr89z9ra6uPnv2bFJSUklJyfz581NSUgQCARExJzJhDxjgXqEAAwSpa9euvfLKK3K5/IMPPhg5cuTGjRs7OzuXLVvGlNJTp04ZjUbPAsycZWQ2m5njz4ykThXbAAAA+klEQVT9+/eXlpYSkc1m83kGANyGAgwQjJxO57PPPmswGDZt2pScnExEq1atys7OPnHixNatW4mouLhYpVJ53oMhKSnpkUceaWlpyc7OfuuttzZu3DhnzpxFixZFRkYSUXt7O1u5AHAUVkEDBKONGze++eabubm5J0+e5PF4TGN1dXVmZqbb7a6oqFi4cGFmZmZhYaHnq/R6/caNGw8dOtTU1BQeHp6cnPzss8/m5OSkpaVNnTr1zJkzbKQCwFUowAAAACzAIWgAAAAWoAADAACwAAUYAACABSjAAAAALEABBgAAYAEKMAAAAAtQgAEAAFiAAgwAAMCC/wUd1/dgchDUFgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%gnuplot\n",
"set style data line\n",
"set yrange [0:2]\n",
"set xrange [-75:75]\n",
"set xlabel \"x/a\"\n",
"set ylabel \"Psi^2\"\n",
"set title \"Scattering of a Gaussian wave packet (41 waves), E=V/4\"\n",
"plot '<./packet 0' t \"time = 0\",'<./packet 40' t \"time = 40\",'Vo.dat' t \"V(x)\",\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework: Developing the numerical part of the C code\n",
"\n",
"Start by copying the skeleton of packet.c
code from above, so that you can start making changes, recompile, and re-run. If you continue working in this jupyter notebook, the next two cells will get executed repeatedly. You can also open a separate editor window to modify the packet.c
file and use the make
to re-build the code after every change."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting packet.c\n"
]
}
],
"source": [
"%%file packet.c\n",
"/*\n",
" * packet.c\n",
" * packet - generate a Gaussian wavepacket impacting on an energy barrier.\n",
" *\n",
" * Completed: January.2018 (c) E.Sternin\n",
" * Revisions: \n",
" *\n",
" */\n",
"\n",
"#ifndef VERSION /* date-derived in Makefile */\n",
"#define VERSION \"2018.01\" /* default, that's when we first wrote the program */\n",
"#endif\n",
"\n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"\n",
"#define MAX_STR 256\n",
"\n",
"/* Global variables */\n",
"\n",
"static char whoami[MAX_STR] ; /* argv[0] will end up here */\n",
"static int verbosity; /* Show detailed information */\n",
"static char options[] = \"Vvhp:t:\"; /* command-line options, : means takes a value */\n",
"static char help_msg[] = \"\\\n",
"%s []\\n\\\n",
"\\t-V report version information\\n\\\n",
"\\t-v increase verbosity, may combine several\\n\\\n",
"\\t-p # number of points for the packet\\n\\\n",
"\\t-t # time since the beginning\\n\\\n",
"\\t-h help message\\n\\\n",
"\\n\\\n",
"e.g.\\tpacket -v -p 601 -t 20\\n\" \n",
";\n",
"\n",
"/*************************************************service routines************/\n",
"\n",
"void __attribute__((noreturn)) die(char *msg, ...) {\n",
" va_list args;\n",
" va_start(args, msg);\n",
" vfprintf(stderr, msg, args);\n",
" fputc('\\n', stderr);\n",
" exit(1);\n",
" }\n",
"\n",
"/***************************************************************** main *********/\n",
"int main(int argc, char **argv) {\n",
" int i,p;\n",
" double t;\n",
"\n",
"/*\n",
" * default values, may get changed by the command-line options\n",
" */\n",
" verbosity = 0; /* 0 = quiet by default, 1 = info, 2 = debug */\t\n",
" p = 25; /* default to 25 points in the packet */\n",
" t = 0;\n",
"\n",
" strncpy(whoami, argv[0], MAX_STR);\n",
"\n",
" while ((i = getopt(argc, argv, options)) != -1)\n",
" switch (i) {\n",
" case 'V':\n",
" printf(\" %s: version %s\\n\",whoami,VERSION);\n",
" break;\n",
" case 'v':\n",
" verbosity++;\n",
" if (verbosity > 1) printf(\" %s: verbosity level set to %d\\n\",whoami,verbosity);\n",
" break;\n",
" case 'h':\n",
" die(help_msg,whoami);\n",
" break;\n",
" case 'p':\n",
" if ( ((p = atoi(optarg)) > 10000) || p < 10 )\n",
" die(\" %s: -p %d is not a valid number of points (10..10000)\\n\",whoami,p);\n",
" if (verbosity > 1) printf(\" %s: Number of points = %d\\n\",whoami,p);\n",
" break;\n",
" case 't':\n",
" if ( ((t = atof(optarg)) < 0) )\n",
" die(\" %s: -t %d is not a valid time\\n\",whoami,t);\n",
" if (verbosity > 1) printf(\" %s: time = %f\\n\",whoami,t);\n",
" break;\n",
" default:\n",
" if (verbosity > 0) die(\" try %s -h\\n\",whoami);\t/* otherwise, die quietly */\n",
" return 0;\n",
" }\n",
"\n",
"/*\n",
" * when we get here, we parsed all user input, and are ready to calculate things\n",
" */\n",
"\n",
" for (i=0; i < p; i++) { \n",
" if (verbosity > 0) printf(\"%d\\t%f\\n\",i,i*t);\n",
" }\n",
"\n",
" return 0;\n",
" }\n"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rm -f packet *.o *~ core\n",
"cc -O -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" -c packet.c\n",
"cc -O -DVERSION=\"\\\"`date '+%Y-%b-%d-%H:%M'`\\\"\" -o packet packet.o -lm\n",
" ./packet: verbosity level set to 2\n",
" ./packet: verbosity level set to 3\n",
" ./packet: Number of points = 10\n",
"0\t0.000000\n",
"1\t0.000000\n",
"2\t0.000000\n",
"3\t0.000000\n",
"4\t0.000000\n",
"5\t0.000000\n",
"6\t0.000000\n",
"7\t0.000000\n",
"8\t0.000000\n",
"9\t0.000000\n"
]
}
],
"source": [
"%%bash\n",
"make clean; make\n",
"./packet -vvv -p10"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
},
"toc": {
"nav_menu": {},
"number_sections": false,
"sideBar": true,
"skip_h1_title": false,
"toc_cell": false,
"toc_position": {},
"toc_section_display": "block",
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}