{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. C programming\n",
    "\n",
    "<div align=\"right\">Edward Sternin</br>September, 2025</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "WORKDIR=~/5P10/Lab2\n",
    "# if the directory exists, remove it and all its contents\n",
    "if [ -d $WORKDIR ]; then\n",
    "  rm -rf $WORKDIR\n",
    "fi\n",
    "# and re-create it as a blank work directory\n",
    "mkdir -p $WORKDIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/edik/5P10/Lab2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'/home/edik/5P10/Lab2'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now tell jupyter where to work\n",
    "%cd ~/5P10/Lab2\n",
    "%pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code is needed only if you are using the notebook, because <tt>make</tt> 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": [
    "## Physics of tunnelling\n",
    "\n",
    "Review <a href=\"https://www.physics.brocku.ca/Courses/5P10/Lectures/wavepackets.pdf\" target=\"_blank\">these notes from a QM course (PDF)</a>.\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",
    "Last time, we also built a simple skeleton C program. 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": "markdown",
   "metadata": {},
   "source": [
    "The code below is our current draft, already slightly enhanced compared to last session.  We have added a couple of switches and some error chacking and parsing of parameters.  We have converted the problem into a dimensionless form by choosing apropriate energy and spatial units, and we have control over output of complex wavefunction Psi, which is all zeros to begin with. "
   ]
  },
  {
   "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 <stdarg.h>\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <string.h>\n",
    "#include <ctype.h>\n",
    "#include <unistd.h>\n",
    "#include <sys/types.h>\n",
    "#include <sys/stat.h>\n",
    "#include <fcntl.h>\n",
    "#include <math.h>\n",
    "#include <complex.h>\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 [<options>]\\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",
    "<pre>\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",
    "</pre>\n",
    "\n",
    "C:\n",
    "<pre>\n",
    "...\n",
    "#include &lt;math.h&gt;\n",
    "#include &lt;complex.h&gt;\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",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortran:\n",
    "<pre>\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",
    "</pre>\n",
    "\n",
    "C:\n",
    "<pre>\n",
    "#include &lt;complex.h&gt;\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",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortran:\n",
    "<pre>\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",
    "</pre>\n",
    "\n",
    "C:\n",
    "<pre>\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",
    "</pre>"
   ]
  },
  {
   "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.12.3"
  },
  "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
}