3. C programming

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.

Gaussian packet evolution: a C project

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:

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.

Physics of tunnelling

Review these notes from a QM course (PDF).

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".

The beginnings of a C program

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.

We will need a Makefile and a gnuplot kernel loaded, for inline graphs.

We can even test the gnuplot script, but of course for all t values, the output is zero everywhere.

Notes on Fortran-to-C conversion

Fortran:

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

C:

...
#include <math.h>
#include <complex.h>

#define NPTS    601       /* these cannot be changed via command-line switches */
#define X_MIN   -75.0
#define X_MAX    75.0
...
  double dx,x[NPTS], ...;
  double complex Psi[NPTS];
...
  dx = (X_MAX-X_MIN)/(double)(NPTS-1);
  for (j=0; j < NPTS; j++) { 
     x[j]=X_MIN + j*dx;
     Psi[j] = (double complex)0.0;  /* start with a blank Psi(x) */
     }

Fortran:

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

C:

#include <complex.h>
...
  for (i=-n; i < n; i++) {
    p_now = p0 + (double)i * p_step;
    E = p_now*p_now;
    A = A0*exp(-pow((p_now-p0)/(2*s),2))*cexp(-I*(p_now*x0 + E*t));

    k1 = csqrt((double complex)(E));
    k2 = csqrt((double complex)(E-U));
    ce4 = cexp(4*I*k2);
    C1 = -2*k1*k2-k1*k1-2*k1*k2*ce4+ce4*(k1*k1+k2*k2)-k2*k2;

    b = (k1*k1-k2*k2)/ C1 *(cexp(2*I*(2*k2-k1))-cexp(-2*I*k1));
    c = -2*k1*(k2+k1)/ C1 *cexp(I*(k2-k1));
    d = (-2*k2+2*k1)*k1/ C1 *cexp(I*(-k1+3*k2));
    f = -4*k1*k2/ C1 *cexp(2*I*(k2-k1));

    ...
  }

Fortran:

    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

C:

    for (j=0; j < NPTS; j++) {
       if (x[j] < -1.) {         /* to the left of the barrier, x < -1 */
          Psi[j] += A*( cexp(I*k1*x[j]) + b*cexp(-I*k1*x[j]) );
          }
       else if (x[j] > 1.) {     /* to the right of the barrier, x > 1 */
          Psi[j] += A*f*cexp(I*k1*x[j]);
          }
       else {                    /* inside the barrier, -1 < x < 1 */
          Psi[j] += A*( c*cexp(I*k2*x[j]) + d*cexp(-I*k2*x[j]) );
          }
       }