Simple wavefunction example with potential wall

wavefunc-my.c

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

const double m = 1;
const double h = 1;
const double U = 1.5;

complex double A(complex double E, double a) {
	complex double SE = csqrt(E);
	complex double z = csqrt(E-U);
	double k = sqrt(2.*m) * a / h;

	return (U/2.) / ((E-U/2) + I*SE*z*1./ctan(k*z));
}

complex double G(complex double E, double a) {
	complex double SE = csqrt(E);
	complex double z = csqrt(E-U);
	complex double k1, _k1;

	k1 = csqrt(2.*m*E) / h;
	_k1 = csqrt(2.*m*(E-U)) / h;


	return 2.*I*k1*_k1 /
		((k1*k1 + _k1*_k1) * csin(_k1*a) + 2.*I*k1*_k1*ccos(_k1*a));
}

complex double B(complex double E, double a) {
	complex double k = csqrt(E / (E - U));
	return 0.5 * (1 + k + A(E, a) * (1 - k));
}

complex double C(complex double E, double a) {
	complex double k = csqrt(E / (E - U));
	return 0.5 * (1 - k + A(E, a) * (1 + k));
}

double before(double x, double E, double a) {
	complex double k = csqrt(2 * m * E /h/h);
	complex double f = cexp(I*k*x) + A(E, a) *cexp(-I*k*x);
	return pow(cabs(f), 2);
}

double inside(double x, double E, double a) {
	complex double _k = csqrt(2 * m * (E - U) / h / h);
	complex double f = B(E, a) * cexp(I*_k*x) + C(E, a) * cexp(-I*_k*x);
	return pow(cabs(f), 2);
}

double after(double x, double E, double a) {
	complex double k = csqrt(2 * m * E /h/h);
	complex double f = G(E, a) * cexp(I*k*(x-a));
	return pow(cabs(f), 2);
}

int main() {
	double xinit, xend;

	double a = 1.5;
	int i; int N = 500;
	xinit = -2;
	xend = a + 2;

	double step = (xend - xinit) / N;
	double x = xinit;
	double y;

	FILE *f = fopen("plt.txt", "w");
	double E = 5;
	for (i=0; i<N; i++) {
		if (x < 0) y = before(x, E, a);
		if (x >= 0 && x < a) y = inside(x, E, a);
		if (x >= a) y = after(x, E, a);

		fprintf(f, "%lf %lf\n", x, y);
		x += step;
	}
	fclose(f);

	return 0;
}