2. Shell and other kinds of programming

Shell programming: advanced parts of the bash tutorial

As you go through the tutorial, copy the example commands into your jupyter notebook and test them to make sure you understand how they work. Feel free to do that more than once, changing the given examples as you wish. An example from the "Using Shell Arrays" section of the tutorial would look like this:

%%bash

NAME[0]="Zara"
NAME[1]="Qadir"
NAME[2]="Mahnaz"
NAME[3]="Ayan"
NAME[4]="Daisy"
echo "First Method: ${NAME[*]}"
echo "Second Method: ${NAME[@]}"

First Method: Zara Qadir Mahnaz Ayan Daisy
Second Method: Zara Qadir Mahnaz Ayan Daisy

You can also save the commands into a shell script and execute it as a "command" in your own filespace:

%%bash
mkdir -p ~/5P10/Lab2

%cd ~/5P10/Lab2

/home/esternin/5P10/Lab2

%%file list.sh
#!/bin/bash

NAME[0]="Zara"
NAME[1]="Qadir"
NAME[2]="Mahnaz"
NAME[3]="Ayan"
NAME[4]="Daisy"
echo "First Method: ${NAME[*]}"
echo "Second Method: ${NAME[@]}"

Overwriting list.sh

%%bash
chmod 755 list.sh
./list.sh

First Method: Zara Qadir Mahnaz Ayan Daisy
Second Method: Zara Qadir Mahnaz Ayan Daisy

%%bash
echo "Hello, world!"

### These are all of the bash environment variables (LOTS OF LINES!):
#set

### These are the variables that have "TRIUMF" in their name:
set | grep TRIUMF

### Select a specific variable and parse the NAME=VALUE for VALUE:
set | grep TRIUMF_FONTS | cut -d= -f2

### instead of typing it, assign it to a variable:
FDIR=`set | grep TRIUMF_FONTS | cut -d= -f2`

echo "Fonts are in \"$FDIR\", look there"

Hello, world!
TRIUMF_FONTS=/usr/local/lib
TRIUMF_PLOTTER_TYPE=PSA
TRIUMF_TERMINAL_TYPE=X
/usr/local/lib
Fonts are in "/usr/local/lib", look there

%%bash
VARS=`set | grep TRIUMF`

for v in $VARS
do
  name=`echo $v | cut -d= -f1`
  value=`echo $v | cut -d= -f2`
  echo "Environment variable \"$name\" has a value \"$value\" "
done

Environment variable "TRIUMF_FONTS" has a value "/usr/local/lib" 
Environment variable "TRIUMF_PLOTTER_TYPE" has a value "PSA" 
Environment variable "TRIUMF_TERMINAL_TYPE" has a value "X" 

Writing and compiling programs

Let's try a very simple C program. "Running" a program requires that you first

• compile it, i.e. convert the human-readable source into machine codes, or object;
• link object to the system libraries to make a complete executable;
• ensure the resulting file has executable flag set.

Usually, the compiler can take care of all of these steps at once.

%%file hello.c
#include <stdio.h>
#include <math.h>

int main() {
  double f;
    
  f=atan(3.);
  printf("Hello, world!  atan(3)=%f\n",f);
    
  return 0; /* a non-zero return reports an error to OS */
  }

Overwriting hello.c

%%bash
cc hello.c
./a.out

Hello, world!  atan(3)=1.249046

The above command cc (that stands for "C Compiler") has many switches. Without them, it defaults to a full compile-and-link cycle and places the executable into an executable called a.out, which can get confusing if there is more than one program. This is better:

%%bash
cc -o hello hello.c
./hello

Hello, world!  atan(3)=1.249046

The link stage of the cc command found all the libraries it needed, including the system library that provided the atan() function.

Homework

Use bash commands ar, nm, grep, find, cut, and bash loops as needed to find which of the *.a or *.so files in /usr/lib64 and/or in /lib and/or /lib64 and all of their subdirectories contains a subroutine to calculate arctan(). Ignore symbolic link files, just use the actual library file, to remove duplication.

Solution

%%bash
# executables in this directory become personal "commands"
if [ ! -d ~/bin ]; then
  mkdir ~/bin
fi

%%file ~/bin/search_libs
#!/bin/sh

# could make it into a parameter
LIBDIR=/usr/lib

# check input line for parameters, at least one is mandatory
if [ $# -lt 1 ]; then
  echo "  usage: $0 subroutine_name"
  exit 1
fi

# if more than one, warn the user we are ignoring the extra input
if [ $# -gt 1 ]; then
  echo "  $0 warning: extra input \"$2 ...\" ignored"
fi

FILES=`find -P $LIBDIR \( -name "*.a" -o -name "*.so" \)`

for f in $FILES
do
   # figure out the file type, apply appropriate command
   ELF=`file $f | grep "shared object"`
   AR=`file $f | grep "ar archive"`
    
   # if both of the above strings are empty, this is not a file for us 
   if [ "$ELF" != "" ]; then
     nm -D --defined-only --print-file $f | grep "$1\$"
   elif [ "$AR" != "" ]; then
     FOUND=`ar -t $f | grep "$1"`
     if [ ! -z "$FOUND" ]; then # another way of checking if string is empty
        echo "$f:$FOUND"
     fi
   fi
done
    
exit 0

Overwriting /home/esternin/bin/search_libs

%%bash
chmod 755 ~/bin/search_libs

%%bash 
search_libs atan

/usr/lib/jvm/java-1.7.0-openjdk-1.7.0.261-2.6.22.2.el7_8.x86_64/jre/lib/amd64/libjava.so:0000000000011620 T Java_java_lang_StrictMath_atan
/usr/lib/jvm/java-1.8.0-openjdk-1.8.0.382.b05-1.el7_9.x86_64/jre/lib/amd64/libjava.so:00000000000127f0 T Java_java_lang_StrictMath_atan
/usr/lib/i686/nosegneg/libm-2.17.so:0000ab20 W atan
/usr/lib/i686/nosegneg/libm-2.17.so:000153e0 W catan
/usr/lib/libm-2.17.so:0000a8e0 W atan
/usr/lib/libm-2.17.so:00015190 W catan
/usr/lib/jalbum/jre64/lib/libjava.so:00000000000146d0 T Java_java_lang_StrictMath_atan

%%bash 
search_libs start_main

/usr/lib/i686/nosegneg/libc-2.17.so:0001a200 T __libc_start_main
/usr/lib/libc-2.17.so:0001a1e0 T __libc_start_main

We need to link to the libm*.so library to get the atan() function, and to libc*.so to initiate the running of a main(). It so happens that these are on the default list, so specifying any libraries is not necessary. However, the following should fail to link, and should produce a few errors:

%%bash
cc -nodefaultlibs -o hello hello.c

/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../lib64/crt1.o: In function `_start':
(.text+0x12): undefined reference to `__libc_csu_fini'
/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../lib64/crt1.o: In function `_start':
(.text+0x19): undefined reference to `__libc_csu_init'
/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../lib64/crt1.o: In function `_start':
(.text+0x25): undefined reference to `__libc_start_main'
/tmp/cc9zQUWU.o: In function `main':
hello.c:(.text+0x2e): undefined reference to `printf'
collect2: error: ld returned 1 exit status

---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
<ipython-input-16-3360b2085b27> in <module>
----> 1 get_ipython().run_cell_magic('bash', '', 'cc -nodefaultlibs -o hello hello.c\n')

/usr/local/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
   2369             with self.builtin_trap:
   2370                 args = (magic_arg_s, cell)
-> 2371                 result = fn(*args, **kwargs)
   2372             return result
   2373 

/usr/local/lib/python3.6/site-packages/IPython/core/magics/script.py in named_script_magic(line, cell)
    140             else:
    141                 line = script
--> 142             return self.shebang(line, cell)
    143 
    144         # write a basic docstring:

/usr/local/lib/python3.6/site-packages/decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

/usr/local/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188 
    189         if callable(arg):

/usr/local/lib/python3.6/site-packages/IPython/core/magics/script.py in shebang(self, line, cell)
    243             sys.stderr.flush()
    244         if args.raise_error and p.returncode!=0:
--> 245             raise CalledProcessError(p.returncode, cell, output=out, stderr=err)
    246 
    247     def _run_script(self, p, cell, to_close):

CalledProcessError: Command 'b'cc -nodefaultlibs -o hello hello.c\n'' returned non-zero exit status 1.

To fix the errors, we need to tell the linker where to look. The compiler/linker switch -lxxx forces it to look for a file libxxx.so.\* in the usual directories. The directory to look in can be also set, using the -L/some_directory/subdir switch. Strictly speaking, in our case we should use -L/lib -lm-2.17 but

• by default the compiler is already looking in /lib, and
• the system has a symbolic link, /lib/libm.so.6 -> /lib/libm-2.17.so

so this (simplified) version without specific version numbers should be sufficient:

%%bash
cc -nodefaultlibs -o hello hello.c -lm -lc
./hello

Hello, world!  atan(3)=1.249046

Using make

The make command treats spaces and TAB characters differently than the jupyter notebook. The following javascript commands should reconfigure jupyter to treat TABS as characters, not replace them with 4 spaces, as it does by default. This is necessary because the syntax of Makefiles mandates the use of TABs.

%%javascript
IPython.tab_as_tab_everywhere = function(use_tabs) {
    if (use_tabs === undefined) {
        use_tabs = true; 
    }
    // apply setting to all current CodeMirror instances
    IPython.notebook.get_cells().map(
        function(c) {  return c.code_mirror.options.indentWithTabs=use_tabs;  }
    );
    // make sure new CodeMirror instances created in the future also use this setting
    CodeMirror.defaults.indentWithTabs=use_tabs;
    };
IPython.tab_as_tab_everywhere()

%%file Makefile
hello: hello.c
	cc -o hello hello.c

Overwriting Makefile

%%bash
rm -f hello
make
./hello

cc -o hello hello.c
Hello, world!  atan(3)=1.249046

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.

For details, consult this makefile tutorial.

%%file Makefile

OBJ = hello.o

# these are pre-defined, but we can also change the default behaviour, if we have multiple compilers
CC = gcc
#CC = icc
CFLAGS = -O
LIBS = -lm

# special abbreviations: $@ is the file to be made; $? is the changed dependent(s).

hello: $(OBJ)
#	@echo -n " 2.link:    "
	$(CC) $(CFLAGS) -o $@ $(OBJ) $(LIBS)

# this is the default implicit rule to convert xxx.c into xxx.o 
.c.o:
#	@echo -n " 1.compile: "
	$(CC) $(CFLAGS) -c $<

clean:
#	@echo -n " 0.cleanup: "
	rm -f hello a.out *.o *~ core

Overwriting Makefile

%%bash
make clean; make
./hello

rm -f hello a.out *.o *~ core
gcc -O -c hello.c
gcc -O -o hello hello.o -lm
Hello, world!  atan(3)=1.249046

Gaussian packet evolution: a C project

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:

• make a project subdirectory
• create a basic skeleton of a program that can read and parse user input (though does not compute anything at first)
• have a basic Makefile in place

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.

%%bash
rm -rf packet
mkdir packet

%cd packet

/home/esternin/5P10/Lab2/packet

%%file packet.c
/*
 * packet.c
 *   packet - generate a Gaussian wavepacket impacting on an energy barrier.
 *
 * Completed: January.2018 (c) E.Sternin
 * Revisions: 
 *
 */

#ifndef VERSION           /* date-derived in Makefile */
#define VERSION "2018.01" /* default, that's when we first wrote the program */
#endif

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

#define MAX_STR 256

/* Global variables */

static char whoami[MAX_STR] ;        /* argv[0] will end up here */
static int  verbosity;               /* Show detailed information */
static char options[] = "Vvhp:t:";   /* command-line options, : means takes a value */
static char help_msg[] = "\
%s [<options>]\n\
\t-V        report version information\n\
\t-v        increase verbosity, may combine several\n\
\t-p #      number of points for the packet\n\
\t-t #      time since the beginning\n\
\t-h        help message\n\
\n\
e.g.\tpacket -v -p 601 -t 20\n" 
;

/*************************************************service routines************/

void __attribute__((noreturn)) die(char *msg, ...) {
  va_list args;
  va_start(args, msg);
  vfprintf(stderr, msg, args);
  fputc('\n', stderr);
  exit(1);
  }

/***************************************************************** main *********/
int main(int argc, char **argv) {
  int    i,p;
  double t;

/*
 *  default values, may get changed by the command-line options
 */
  verbosity = 0;    /* 0 = quiet by default, 1 = info, 2 = debug */	
  p = 25;           /* default to 25 points in the packet */
  t = 0;

  strncpy(whoami, argv[0], MAX_STR);

  while ((i = getopt(argc, argv, options)) != -1)
    switch (i) {
      case 'V':
         printf(" %s: version %s\n",whoami,VERSION);
         break;
      case 'v':
         verbosity++;
         if (verbosity > 1) printf(" %s: verbosity level set to %d\n",whoami,verbosity);
         break;
      case 'h':
         die(help_msg,whoami);
         break;
      case 'p':
         if ( ((p = atoi(optarg)) > 10000) || p < 10 )
           die(" %s: -p %d is not a valid number of points (10..10000)\n",whoami,p);
         if (verbosity > 1) printf(" %s: Number of points = %d\n",whoami,p);
         break;
      case 't':
         if ( ((t = atof(optarg)) < 0) )
           die(" %s: -t %d is not a valid time\n",whoami,t);
         if (verbosity > 1) printf(" %s: time = %f\n",whoami,t);
         break;
      default:
         if (verbosity > 0) die(" try %s -h\n",whoami);	/* otherwise, die quietly */
         return 0;
      }

/*
 * when we get here, we parsed all user input, and are ready to calculate things
 */

  for (i=0; i < p; i++) { 
     if (verbosity > 0) printf("%d\t%f\n",i,i*t);
     }

  return 0;
  }

Writing packet.c

%%bash
cc -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\"" -o packet packet.c
./packet -vvvV -p10

 ./packet: verbosity level set to 2
 ./packet: verbosity level set to 3
 ./packet: version 2023-Sep-12-12:52
 ./packet: Number of points = 10
0	0.000000
1	0.000000
2	0.000000
3	0.000000
4	0.000000
5	0.000000
6	0.000000
7	0.000000
8	0.000000
9	0.000000

%%file Makefile

OBJ = packet.o

# these are pre-defined, but we can also change the default behaviour, if we have multiple compilers
CC = cc
#CC = icc
CFLAGS = -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\"" 
LIBS = -lm

# special abbreviations: $@ is the file to be made; $? is the changed dependent(s).

packet: $(OBJ)
	$(CC) $(CFLAGS) -o $@ $(OBJ) $(LIBS)

# this is the default implicit rule to convert xxx.c into xxx.o 
.c.o:
	$(CC) $(CFLAGS) -c $<

clean:
	rm -f packet *.o *~ core

Writing Makefile

%%bash
make clean; make
./packet -Vvvv -p 11 -t 25

rm -f packet *.o *~ core
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -c packet.c
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -o packet packet.o -lm
 ./packet: version 2023-Sep-12-12:52
 ./packet: verbosity level set to 2
 ./packet: verbosity level set to 3
 ./packet: Number of points = 11
 ./packet: time = 25.000000
0	0.000000
1	25.000000
2	50.000000
3	75.000000
4	100.000000
5	125.000000
6	150.000000
7	175.000000
8	200.000000
9	225.000000
10	250.000000

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.

%%file packet.f
C packet.f -- propagation of a Gaussian wave  packet through a barrier in 1D
C	x = coordinate (in units of a, the barrier is from -a to +a)
C	P = wavefunction  (not normalized)
C	E = energy of the wavefunction (in units of V, the barrier height)
C
C	Written by:	E.Sternin
C	Completed:	30.03.93
C=======================================================================	
	real*4		x(601),x0,p0,sp,p_now,p_step,Pi,A0,t
	complex*8	P(601),I,A
	integer		j,Nargs,Nwaves
	character*80 	buffer

	Nwaves = 20	! add up (2*Nwaves+1) eigenwaves to form a packet
	x0 = -25.	! start the packet far away from the barrier
	p0 = 0.5	! keep incident energy low, E=0.5**2 < 1 (units of V)
	sp = 0.1	! p0 +/- 4 sigma_p > 0,	all eigenwaves move to the right
	 
	Nargs = iargc() 	! arguments supplied on the command line ?
	if (Nargs .lt. 1) then 
	  write(*,'(A,$)') '  Error: need time t: '
	  read(*,*) t
	else
	  call getarg(1,buffer)
	  read(buffer,*) t
	end if

	I = cmplx(0.,1.)		!complex i
	Pi = 2.*acos(0.)
	A0 = 1./sqrt( sqrt(2.*Pi) * sp ) /float (2*Nwaves+1)

	do j = 1,601
	  x(j) = -75 + (j-1)*0.25
	  P(j) = cmplx(0.,0.)
	end do

	p_step = 4.*sp/float(Nwaves)	! calculate over +/- 4*sigma_p
	do j = -Nwaves,Nwaves		! in Nwaves steps
	  p_now = p0 + j * p_step 
	  E = p_now**2
	  A = A0*exp( - ( (p_now-p0)/(2*sp) )**2 - I*(p_now*x0 + E*t) )
	  call add_wave(x,P,A,E)
	end do

C..output the saved values as is, without normalization
  10	do j=1,	601
	  write(*,*) x(j),abs(P(j))**2,real(P(j)),aimag(P(j))
	end do
	end

C=======================================================================	
	subroutine add_wave(x,P,A,E)
C..adds the specified wave to the the supplied complex storage array P
C  the values of Psi are only calculated at points given in real array x 
	real*4 		x(601),E
	complex*8	P(601),I,A,k1,K2,C1,b,c,d,f

	I = cmplx(0.,1.)	!complex i

	k1 = sqrt(cmplx(E,0.))
	K2 = sqrt(cmplx(E-1.,0.))

	C1 = -2*k1*K2-k1**2-2*k1*K2*exp(4*I*K2)+exp(4*I*K2)*k1**2
     *	-K2**2+K2**2*exp(4*I*K2)

	b = (k1**2-K2**2)/ C1 *(exp(2*I*(2*K2-k1))-exp(-2*I*k1))
	c = -2*k1*(K2+k1)/ C1 *exp(I*(K2-k1))
	d = (-2*K2+2*k1)*k1/ C1 *exp(I*(-k1+3*K2))
	f = -4*k1*K2/ C1 *exp(2*I*(K2-k1))

	do j = 1,296		! left of the barrier, x < -1
	  P(j) = P(j) + A * ( exp(I*k1*x(j)) + b*exp(-I*k1*x(j)) )
	end do
	do j=297,304		! inside the barrier, -1 < x < 1
	  P(j) = P(j) + A * ( c*exp(I*K2*x(j)) + d*exp(-I*K2*x(j)) )
	end do
	do j=305,601		! to the right of the barrier, x > 1
	  P(j) = P(j) + A * f*exp(I*k1*x(j))
	end do
	  
 	return
	end

Writing packet.f

%%bash
gfortran -o packet packet.f
./packet 40

  -75.0000000       8.52197530E-08   1.16693496E-04   2.67586205E-04
  -74.7500000       3.20557263E-07   3.57286714E-04   4.39207768E-04
  -74.5000000       7.47516651E-07   6.27484871E-04   5.94793586E-04
  -74.2500000       1.36120696E-06   9.13982512E-04   7.25150225E-04
  -74.0000000       2.12284544E-06   1.20234303E-03   8.22931761E-04
  -73.7500000       2.96266194E-06   1.47765456E-03   8.82722496E-04
  -73.5000000       3.78836876E-06   1.72579114E-03   9.00007668E-04
  -73.2500000       4.50573361E-06   1.93375722E-03   8.75395373E-04
  -73.0000000       5.02329522E-06   2.09036935E-03   8.08486657E-04
  -72.7500000       5.27878183E-06   2.18690187E-03   7.04444654E-04
  -72.5000000       5.24243342E-06   2.21807300E-03   5.67965733E-04
  -72.2500000       4.92594882E-06   2.18175026E-03   4.07325802E-04
  -72.0000000       4.37627432E-06   2.07901793E-03   2.32291524E-04
  -71.7500000       3.67039229E-06   1.91512029E-03   5.20246103E-05
  -71.5000000       2.90051889E-06   1.69875333E-03  -1.21474615E-04
  -71.2500000       2.15335467E-06   1.44088222E-03  -2.77872314E-04
  -71.0000000       1.50178641E-06   1.15568750E-03  -4.07642918E-04
  -70.7500000       9.90534545E-07   8.59129708E-04  -5.02424780E-04
  -70.5000000       6.29514943E-07   5.67211304E-04  -5.54784900E-04
  -70.2500000       4.03137477E-07   2.97086081E-04  -5.61139314E-04
  -70.0000000       2.74815335E-07   6.43452368E-05  -5.20264381E-04
  -69.7500000       2.00807236E-07  -1.17278614E-04  -4.32496192E-04
  -69.5000000       1.46606908E-07  -2.35196581E-04  -3.02141474E-04
  -69.2500000       9.81526753E-08  -2.82348512E-04  -1.35764451E-04
  -69.0000000       6.73033540E-08  -2.52886268E-04   5.78955514E-05
  -68.7500000       9.34146485E-08  -1.47164857E-04   2.67875264E-04
  -68.5000000       2.34040129E-07   3.21054831E-05   4.82710428E-04
  -68.2500000       5.51280777E-07   2.76413339E-04   6.89112814E-04
  -68.0000000       1.09518010E-06   5.73823112E-04   8.75161262E-04
  -67.7500000       1.88719878E-06   9.09983180E-04   1.02914008E-03
  -67.5000000       2.90695652E-06   1.26707926E-03   1.14081835E-03
  -67.2500000       4.09146287E-06   1.62594335E-03   1.20323373E-03
  -67.0000000       5.33215916E-06   1.96641800E-03   1.21052016E-03
  -66.7500000       6.49898357E-06   2.26926315E-03   1.16164889E-03
  -66.5000000       7.44904264E-06   2.51625036E-03   1.05713122E-03
  -66.2500000       8.05988930E-06   2.69161747E-03   9.02820029E-04
  -66.0000000       8.24550898E-06   2.78319465E-03   7.06637395E-04
  -65.7500000       7.98025758E-06   2.78387824E-03   4.79874376E-04
  -65.5000000       7.28935311E-06   2.68957182E-03   2.35704472E-04
  -65.2500000       6.26574547E-06   2.50312570E-03  -1.03411730E-05
  -65.0000000       5.03992305E-06   2.23173853E-03  -2.43447372E-04
  -64.7500000       3.76202502E-06   1.88760250E-03  -4.46073827E-04
  -64.5000000       2.57725833E-06   1.48712832E-03  -6.04737783E-04
  -64.2500000       1.60647710E-06   1.05198694E-03  -7.06965802E-04
  -64.0000000       9.16446481E-07   6.04671310E-04  -7.42171891E-04
  -63.7500000       5.25519908E-07   1.70712476E-04  -7.04540405E-04
  -63.5000000       4.00786774E-07  -2.23876239E-04  -5.92170749E-04
  -63.2500000       4.74573596E-07  -5.55795734E-04  -4.07019281E-04
  -63.0000000       6.67948825E-07  -8.02147144E-04  -1.56552880E-04
  -62.7500000       9.16697616E-07  -9.45987180E-04   1.47668005E-04
  -62.5000000       1.18827541E-06  -9.73457354E-04   4.90567181E-04
  -62.2500000       1.49675054E-06  -8.77579034E-04   8.52411729E-04
  -62.0000000       1.90388926E-06  -6.57910947E-04   1.21286535E-03
  -61.7500000       2.50078847E-06  -3.20466526E-04   1.54857663E-03
  -61.5000000       3.39026064E-06   1.21302844E-04   1.83726603E-03
  -61.2500000       4.65240237E-06   6.47406210E-04   2.05749064E-03
  -61.0000000       6.31526882E-06   1.23210717E-03   2.19024671E-03
  -60.7500000       8.33706690E-06   1.84477249E-03   2.22123414E-03
  -60.5000000       1.05908302E-05   2.45194626E-03   2.13981071E-03
  -60.2500000       1.28799593E-05   3.01780133E-03   1.94237835E-03
  -60.0000000       1.49629741E-05   3.50721832E-03   1.63168437E-03
  -59.7500000       1.65833972E-05   3.88620817E-03   1.21687388E-03
  -59.5000000       1.75303558E-05   4.12542885E-03   7.14978145E-04
  -59.2500000       1.76757167E-05   4.20163572E-03   1.48241408E-04
  -59.0000000       1.70022122E-05   4.09815973E-03  -4.55302768E-04
  -58.7500000       1.56223559E-05   3.80731258E-03  -1.06147316E-03
  -58.5000000       1.37751495E-05   3.33144492E-03  -1.63603912E-03
  -58.2500000       1.17840091E-05   2.68272613E-03  -2.14172597E-03
  -58.0000000       1.00079678E-05   1.88314915E-03  -2.54199072E-03
  -57.7500000       8.79676918E-06   9.65245708E-04  -2.80447328E-03
  -57.5000000       8.41458677E-06  -3.08033777E-05  -2.90062721E-03
  -57.2500000       9.01049771E-06  -1.05794892E-03  -2.80913548E-03
  -57.0000000       1.05928020E-05  -2.06299732E-03  -2.51730881E-03
  -56.7500000       1.30348853E-05  -2.99088843E-03  -2.02224450E-03
  -56.5000000       1.61138578E-05  -3.78680765E-03  -1.33189501E-03
  -56.2500000       1.95610373E-05  -4.39817924E-03  -4.65892372E-04
  -56.0000000       2.31367121E-05  -4.77899890E-03   5.45784831E-04
  -55.7500000       2.66794759E-05  -4.89073480E-03   1.66138215E-03
  -55.5000000       3.01715318E-05  -4.70670173E-03   2.83169374E-03
  -55.2500000       3.37367819E-05  -4.21285350E-03   3.99858039E-03
  -55.0000000       3.76395037E-05  -3.41027766E-03   5.09995222E-03
  -54.7500000       4.22198063E-05  -2.31518969E-03   6.07121922E-03
  -54.5000000       4.78313013E-05  -9.60440026E-04   6.84900396E-03
  -54.2500000       5.47439013E-05   6.05673355E-04   7.37408036E-03
  -54.0000000       6.30851209E-05   2.32096715E-03   7.59593491E-03
  -53.7500000       7.27457809E-05   4.11053933E-03   7.47323502E-03
  -53.5000000       8.34101229E-05   5.89031586E-03   6.97956327E-03
  -53.2500000       9.45698048E-05   7.57003017E-03   6.10446092E-03
  -53.0000000       1.05608735E-04   9.05702077E-03   4.85583209E-03
  -52.7500000       1.15908770E-04   1.02605270E-02   3.26042413E-03
  -52.5000000       1.24996470E-04   1.10966908E-02   1.36378640E-03
  -52.2500000       1.32657777E-04   1.14920577E-02  -7.68357306E-04
  -52.0000000       1.39047828E-04   1.13893524E-02  -3.05458345E-03
  -51.7500000       1.44682970E-04   1.07490271E-02  -5.39827580E-03
  -51.5000000       1.50450665E-04   9.55420081E-03  -7.69206742E-03
  -51.2500000       1.57490605E-04   7.81247579E-03  -9.82119329E-03
  -51.0000000       1.67064078E-04   5.55775221E-03  -1.16694244E-02
  -50.7500000       1.80347270E-04   2.84943869E-03  -1.31235654E-02
  -50.5000000       1.98269292E-04  -2.27141310E-04  -1.40789812E-02
  -50.2500000       2.21380440E-04  -3.56273819E-03  -1.44460145E-02
  -50.0000000       2.49718258E-04  -7.02941464E-03  -1.41529366E-02
  -49.7500000       2.82896624E-04  -1.04816696E-02  -1.31541332E-02
  -49.5000000       3.20124294E-04  -1.37642613E-02  -1.14310738E-02
  -49.2500000       3.60413804E-04  -1.67173911E-02  -8.99681263E-03
  -49.0000000       4.02785983E-04  -1.91835742E-02  -5.89715503E-03
  -48.7500000       4.46478196E-04  -2.10139658E-02  -2.21166899E-03
  -48.5000000       4.91187500E-04  -2.20769364E-02   1.94842578E-03
  -48.2500000       5.37160086E-04  -2.22638827E-02   6.44046534E-03
  -48.0000000       5.85187867E-04  -2.14950033E-02   1.10974200E-02
  -47.7500000       6.36558922E-04  -1.97261460E-02   1.57301631E-02
  -47.5000000       6.92870934E-04  -1.69525705E-02   2.01365650E-02
  -47.2500000       7.55716639E-04  -1.32109318E-02   2.41078399E-02
  -47.0000000       8.26478878E-04  -8.58242996E-03   2.74375789E-02
  -46.7500000       9.05986526E-04  -3.19017912E-03   2.99300738E-02
  -46.5000000       9.94412228E-04   2.80074077E-03   3.14096808E-02
  -46.2500000       1.09125383E-03   9.18914936E-03   3.17303203E-02
  -46.0000000       1.19533867E-03   1.57422144E-02   3.07818335E-02
  -45.7500000       1.30525301E-03   2.22040173E-02   2.84997337E-02
  -45.5000000       1.41946936E-03   2.83027254E-02   2.48681568E-02
  -45.2500000       1.53699284E-03   3.37631181E-02   1.99259799E-02
  -45.0000000       1.65757514E-03   3.83147635E-02   1.37678627E-02
  -44.7500000       1.78213743E-03   4.17049080E-02   6.54508593E-03
  -44.5000000       1.91278907E-03   4.37084027E-02  -1.53771020E-03
  -44.2500000       2.05288385E-03   4.41394821E-02  -1.02269184E-02
  -44.0000000       2.20666802E-03   4.28607017E-02  -1.92257129E-02
  -43.7500000       2.37896340E-03   3.97920758E-02  -2.82055680E-02
  -43.5000000       2.57440354E-03   3.49168628E-02  -3.68132629E-02
  -43.2500000       2.79705203E-03   2.82865763E-02  -4.46869284E-02
  -43.0000000       3.04991286E-03   2.00244263E-02  -5.14678098E-02
  -42.7500000       3.33440048E-03   1.03225214E-02  -5.68141378E-02
  -42.5000000       3.65031790E-03  -5.59103151E-04  -6.04152754E-02
  -42.2500000       3.99616593E-03  -1.22994892E-02  -6.20071627E-02
  -42.0000000       4.36929567E-03  -2.45239753E-02  -6.13829829E-02
  -41.7500000       4.76662349E-03  -3.68160494E-02  -5.84054962E-02
  -41.5000000       5.18542528E-03  -4.87304702E-02  -5.30166626E-02
  -41.2500000       5.62385516E-03  -5.98065779E-02  -4.52441014E-02
  -41.0000000       6.08159089E-03  -6.95862025E-02  -3.52044180E-02
  -40.7500000       6.56005880E-03  -7.76285455E-02  -2.31055636E-02
  -40.5000000       7.06258509E-03  -8.35293457E-02  -9.24300961E-03
  -40.2500000       7.59381475E-03  -8.69352892E-02   6.00590743E-03
  -40.0000000       8.15933105E-03  -8.75608996E-02   2.21905336E-02
  -39.7500000       8.76477081E-03  -8.52020681E-02   3.87991965E-02
  -39.5000000       9.41498019E-03  -7.97474012E-02   5.52750565E-02
  -39.2500000       1.01132030E-02  -7.11869895E-02   7.10324943E-02
  -39.0000000       1.08607309E-02  -5.96186668E-02   8.54771659

limit_output extension: Maximum message size of 10000 exceeded with 41469 characters

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.

%%file Vo.dat
-75	0
-1	0
-1	1
1	1
1	0
75	0

Writing Vo.dat

%%bash
gnuplot --persist
plot 'Vo.dat' with lines

%%bash
gnuplot --persist
set style data line
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"

We can even create an "animation" of sorts, by using a macro file.

%%file packet.gnu
set style data line
set yrange [0:2]
set xrange [-75:75]
set xlabel "x/a"
set ylabel "Psi^2"
set title "Scattering of a Gaussian wave packet (41 waves), E=V/4"
plot '<./packet   0' t "time =  0",'Vo.dat' t "V(x)"
pause 2 "initially, the packet is at x = -25 a"
plot '<./packet   5' t "time =  5",'Vo.dat' t "V(x)"   
pause 2 "as time goes on..."
plot '<./packet  10' t "time = 10",'Vo.dat' t "V(x)"   
pause 2 "...the packet approaches the barrier"
plot '<./packet  15' t "time = 15",'Vo.dat' t "V(x)"   
pause 2 "the interference is observed"
plot '<./packet  20' t "time = 20",'Vo.dat' t "V(x)"   
pause 2 "the interference is observed"
plot '<./packet  25' t "time = 25",'Vo.dat' t "V(x)"   
pause 2 "the interference is observed"
plot '<./packet  30' t "time = 30",'Vo.dat' t "V(x)"   
pause 2 "the interference is observed"
plot '<./packet  35' t "time = 35",'Vo.dat' t "V(x)"   
pause 2 "transmitted and reflected packets form"
plot '<./packet  40' t "time = 40",'Vo.dat' t "V(x)"   
pause 2 "transmitted and reflected packets form"
plot '<./packet  45' t "time = 45",'Vo.dat' t "V(x)"   
pause 2 "and move away from the barrier"
plot '<./packet  50' t "time = 50",'Vo.dat' t "V(x)"   
pause 2 "and move away from the barrier"
plot '<./packet  55' t "time = 55",'Vo.dat' t "V(x)"   
pause 2 "note how they spread out with time"
plot '<./packet  60' t "time = 60",'Vo.dat' t "V(x)"   
pause 2 "note how they spread out with time"
plot '<./packet  65' t "time = 65",'Vo.dat' t "V(x)"   
pause 2 "note how they spread out with time"
plot '<./packet  70' t "time = 70",'Vo.dat' t "V(x)"   
pause 2 "note how they spread out with time"
plot '<./packet  80' t "time = 80",'Vo.dat' t "V(x)"   
pause 2
plot '<./packet  90' t "time = 90",'Vo.dat' t "V(x)"
pause 2

Writing packet.gnu

%%bash
gnuplot
load 'packet.gnu'

initially, the packet is at x = -25 a
as time goes on...
...the packet approaches the barrier
the interference is observed
the interference is observed
the interference is observed
the interference is observed
transmitted and reflected packets form
transmitted and reflected packets form
and move away from the barrier
and move away from the barrier
note how they spread out with time
note how they spread out with time
note how they spread out with time
note how they spread out with time

Another way of using gnuplot is to load it right into the notebook itself.

%load_ext gnuplot_kernel

%gnuplot inline pngcairo font "Arial,16" size 640,480

%%gnuplot
set style data line
set yrange [0:2]
set xrange [-75:75]
set xlabel "x/a"
set ylabel "Psi^2"
set title "Scattering of a Gaussian wave packet (41 waves), E=V/4"
plot '<./packet   0' t "time =  0",'<./packet   40' t "time =  40",'Vo.dat' t "V(x)",

Homework: Developing the numerical part of the C code

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.

%%file packet.c
/*
 * packet.c
 *   packet - generate a Gaussian wavepacket impacting on an energy barrier.
 *
 * Completed: January.2018 (c) E.Sternin
 * Revisions: 
 *
 */

#ifndef VERSION           /* date-derived in Makefile */
#define VERSION "2018.01" /* default, that's when we first wrote the program */
#endif

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

#define MAX_STR 256

/* Global variables */

static char whoami[MAX_STR] ;        /* argv[0] will end up here */
static int  verbosity;               /* Show detailed information */
static char options[] = "Vvhp:t:";   /* command-line options, : means takes a value */
static char help_msg[] = "\
%s [<options>]\n\
\t-V        report version information\n\
\t-v        increase verbosity, may combine several\n\
\t-p #      number of points for the packet\n\
\t-t #      time since the beginning\n\
\t-h        help message\n\
\n\
e.g.\tpacket -v -p 601 -t 20\n" 
;

/*************************************************service routines************/

void __attribute__((noreturn)) die(char *msg, ...) {
  va_list args;
  va_start(args, msg);
  vfprintf(stderr, msg, args);
  fputc('\n', stderr);
  exit(1);
  }

/***************************************************************** main *********/
int main(int argc, char **argv) {
  int    i,p;
  double t;

/*
 *  default values, may get changed by the command-line options
 */
  verbosity = 0;    /* 0 = quiet by default, 1 = info, 2 = debug */	
  p = 25;           /* default to 25 points in the packet */
  t = 0;

  strncpy(whoami, argv[0], MAX_STR);

  while ((i = getopt(argc, argv, options)) != -1)
    switch (i) {
      case 'V':
         printf(" %s: version %s\n",whoami,VERSION);
         break;
      case 'v':
         verbosity++;
         if (verbosity > 1) printf(" %s: verbosity level set to %d\n",whoami,verbosity);
         break;
      case 'h':
         die(help_msg,whoami);
         break;
      case 'p':
         if ( ((p = atoi(optarg)) > 10000) || p < 10 )
           die(" %s: -p %d is not a valid number of points (10..10000)\n",whoami,p);
         if (verbosity > 1) printf(" %s: Number of points = %d\n",whoami,p);
         break;
      case 't':
         if ( ((t = atof(optarg)) < 0) )
           die(" %s: -t %d is not a valid time\n",whoami,t);
         if (verbosity > 1) printf(" %s: time = %f\n",whoami,t);
         break;
      default:
         if (verbosity > 0) die(" try %s -h\n",whoami);	/* otherwise, die quietly */
         return 0;
      }

/*
 * when we get here, we parsed all user input, and are ready to calculate things
 */

  for (i=0; i < p; i++) { 
     if (verbosity > 0) printf("%d\t%f\n",i,i*t);
     }

  return 0;
  }

Overwriting packet.c

%%bash
make clean; make
./packet -vvv -p10

rm -f packet *.o *~ core
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -c packet.c
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -o packet packet.o -lm
 ./packet: verbosity level set to 2
 ./packet: verbosity level set to 3
 ./packet: Number of points = 10
0	0.000000
1	0.000000
2	0.000000
3	0.000000
4	0.000000
5	0.000000
6	0.000000
7	0.000000
8	0.000000
9	0.000000