2. C programming

Edward Sternin
September, 2025

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.

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

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.

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

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.

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]) );
          }
       }