{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 3. C programming" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code is needed only if you are using the notebook, because make command treats spaces and TAB characters differently. We need to be able to reconfigure Jupyter to treat TABS as characters, not replace them with 4 spaces, as it does by default." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "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", "\n", "IPython.tab_as_tab_everywhere()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gaussian packet evolution: a C project" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the first homework to be handed in, 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", "

\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": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "## start with an empty directory for the project\n", "if [ ! -d ~/5P10 ]; then\n", " mkdir ~/5P10\n", "fi\n", "cd ~/5P10\n", "rm -rf packet\n", "mkdir packet" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%cd ~/5P10/packet" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "pwd\n", "ls -l" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Physics of tunnelling\n", "\n", "Review these notes from a QM course (PDF).\n", "\n", "Last session, we saw a working Fortran program. The relevant numerical part is this, and the notation is a straightforward translation of the expressions from the PDF notes. FORTRAN = \"FORmula TRANslation\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%file add_wave.f\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": "markdown", "metadata": {}, "source": [ "## The beginnings of a C program\n", "\n", "Now we are ready to dive into C. We'll start with a simple skeleton C program, slightly enhanced compared to the version discussed the last time. The numeric heart of the program is still missing, and the output consists of meaningless numbers.\n", "\n", "We will need a Makefile and a gnuplot kernel loaded, for inline graphs." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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": null, "metadata": {}, "outputs": [], "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": null, "metadata": {}, "outputs": [], "source": [ "%load_ext gnuplot_kernel" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%gnuplot inline pngcairo font \"Arial,16\" size 800,600" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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", "#include \n", "#include \n", "\n", "#define MAX_STR 256\n", "#define NPTS 601 /* these cannot be changed via command-line switches */\n", "#define X_MIN -75.0\n", "#define X_MAX 75.0\n", "#define NS 4 /* +/- NS * sigma is the extent of the Gaussian packet we calculate */\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[] = \"Vvhn:t:U:p:s:x:\"; /* 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-n # number of waves making up the packet\\n\\\n", "\\t-t # time since the beginning\\n\\\n", "\\t-U # U[=16] height of the barrier U(x), in epsilon\\n\\\n", "\\t-p # p0[=0.5] is the mean momentum of the packet\\n\\\n", "\\t-s # s[=0.1] is the width of the packet in p-space\\n\\\n", "\\t-x # x0[=-25] is the initial position of the packet\\n\\\n", "\\t-h help message\\n\\\n", "\\n\\\n", "e.g.\\tpacket -v -n 20 -t 20 -p 0.5 -s 0.1\\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,j,n;\n", " double t,U,p0,s,x0,dx,x[NPTS],A0,p_step,p_now,E;\n", " double complex Psi[NPTS];\n", "\n", " double Pi = 2.*acos(0.);\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", " n = 20; /* default to 41=(2*n+1) waves in the packet */\n", " p0 = 0.5; /* default to 0.5 in p-space */\n", " s = 0.1; /* keep to less than 1/4 of p0, so that (p0+/-4s) > 0, all */\n", " /* constituent waves travel to the right initially */\n", " U = 1; /* barrier height, in epsilon */\n", " x0 = -25; /* initial position of 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 'n':\n", " if ( ((n = atoi(optarg)) > 10000) || n < 10 )\n", " die(\" %s: -n %d is not a valid number of points (10..10000)\\n\",whoami,n);\n", " if (verbosity > 1) printf(\" %s: Number of points = %d\\n\",whoami,n);\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", " case 'U':\n", " U = atof(optarg);\n", " if (verbosity > 1) printf(\" %s: barrier height is = %f epsilon\\n\",whoami,U);\n", " break;\n", " case 'p':\n", " p0 = atof(optarg);\n", " if (verbosity > 1) printf(\" %s: mean momentum of packet is p0 = %f\\n\",whoami,p0);\n", " break;\n", " case 's':\n", " if ( (s = atof(optarg)) < 1e-2 )\n", " die(\" %s: -s %f is not a valid packet width, >= 1e-2\\n\",whoami,s);\n", " if (verbosity > 1) printf(\" %s: width of packet, s = %f\\n\",whoami,s);\n", " break;\n", " case 'x':\n", " if ( ((x0 = atof(optarg)) < X_MIN) || x0 > 0)\n", " die(\" %s: -x %f is not a valid packet position, %d..0\\n\",whoami,x0,X_MIN);\n", " if (verbosity > 1) printf(\" %s: initial position of packet, x0 = %f\\n\",whoami,x0);\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", " if ((p0-NS*s) <= 0) /* is the smallest value of p still positive ? */\n", " die(\" %s: some waves in the packet p=%f+/-%d*%f are not travelling to the right at t=0\\n\",whoami,p0,NS,s);\n", "\n", " if ((p0+NS*s) >= sqrt(U)) /* is the largest value of p still below U of the barrier ? */\n", " die(\" %s: some waves in the packet p=%f+/-%d*%f have energy >= U=%f\\n\",whoami,p0,NS,s,U);\n", "\n", " dx = (X_MAX-X_MIN)/(double)(NPTS-1);\n", " for (j=0; j < NPTS; j++) { \n", " x[j]=X_MIN + j*dx;\n", " Psi[j] = (double complex)0.0; /* start with a blank Psi(x) */\n", " }\n", "\n", "\n", " /*\n", " * here we compute Psi(x,t), blank for now\n", " */\n", "\n", "\n", " /* output the total Psi as is, without normalization */\n", " for (j=0; j < NPTS; j++) {\n", " if (verbosity > 0) printf(\"\\t%f\\t%f\\t%f\\t%f\\n\",x[j],cabs(Psi[j]),creal(Psi[j]),cimag(Psi[j]));\n", " }\n", "\n", " return 0;\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "make\n", "./packet -v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can even test the gnuplot script, but of course for all t values, the output is zero everywhere." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%gnuplot\n", "set style data line\n", "set yrange [-0.1:1.1]\n", "set xrange [-60:60]\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 \\\n", " 'Vo.dat' t \"V(x)\",\\\n", " '<./packet -v -t 0' title \"t=0\",\\\n", " '<./packet -v -t 20' title \"t=20\",\\\n", " '<./packet -v -t 50' title \"t=50\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Notes on Fortran-to-C conversion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fortran:\n", "

\n",
    "    real*4      x(601), ...\n",
    "\tcomplex*8   P(601), ...\n",
    "\t...\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", "\n", "C:\n", "
\n",
    "...\n",
    "#include <math.h>\n",
    "#include <complex.h>\n",
    "\n",
    "#define NPTS    601       /* these cannot be changed via command-line switches */\n",
    "#define X_MIN   -75.0\n",
    "#define X_MAX    75.0\n",
    "...\n",
    "  double dx,x[NPTS], ...;\n",
    "  double complex Psi[NPTS];\n",
    "...\n",
    "  dx = (X_MAX-X_MIN)/(double)(NPTS-1);\n",
    "  for (j=0; j < NPTS; j++) { \n",
    "     x[j]=X_MIN + j*dx;\n",
    "     Psi[j] = (double complex)0.0;  /* start with a blank Psi(x) */\n",
    "     }\n",
    "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fortran:\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", "\n", "C:\n", "
\n",
    "#include <complex.h>\n",
    "...\n",
    "  for (i=-n; i < n; i++) {\n",
    "    p_now = p0 + (double)i * p_step;\n",
    "    E = p_now*p_now;\n",
    "    A = A0*exp(-pow((p_now-p0)/(2*s),2))*cexp(-I*(p_now*x0 + E*t));\n",
    "\n",
    "    k1 = csqrt((double complex)(E));\n",
    "    k2 = csqrt((double complex)(E-U));\n",
    "    ce4 = cexp(4*I*k2);\n",
    "    C1 = -2*k1*k2-k1*k1-2*k1*k2*ce4+ce4*(k1*k1+k2*k2)-k2*k2;\n",
    "\n",
    "    b = (k1*k1-k2*k2)/ C1 *(cexp(2*I*(2*k2-k1))-cexp(-2*I*k1));\n",
    "    c = -2*k1*(k2+k1)/ C1 *cexp(I*(k2-k1));\n",
    "    d = (-2*k2+2*k1)*k1/ C1 *cexp(I*(-k1+3*k2));\n",
    "    f = -4*k1*k2/ C1 *cexp(2*I*(k2-k1));\n",
    "\n",
    "    ...\n",
    "  }\n",
    "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fortran:\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",
    "
\n", "\n", "C:\n", "
\n",
    "    for (j=0; j < NPTS; j++) {\n",
    "       if (x[j] < -1.) {         /* to the left of the barrier, x < -1 */\n",
    "          Psi[j] += A*( cexp(I*k1*x[j]) + b*cexp(-I*k1*x[j]) );\n",
    "          }\n",
    "       else if (x[j] > 1.) {     /* to the right of the barrier, x > 1 */\n",
    "          Psi[j] += A*f*cexp(I*k1*x[j]);\n",
    "          }\n",
    "       else {                    /* inside the barrier, -1 < x < 1 */\n",
    "          Psi[j] += A*( c*cexp(I*k2*x[j]) + d*cexp(-I*k2*x[j]) );\n",
    "          }\n",
    "       }\n",
    "
" ] }, { "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": { "height": "1033.6px", "left": "0px", "right": "1440.97px", "top": "79.4px", "width": "156.033px" }, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }