Back (Current repo: kronig-penney)

kp.c (1786B)

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

void solve(const double E, double complex k[2], double b, double V_0, double Delta) {
    double complex q = csqrt(2.0 * E);
    double complex kappa = csqrt(2.0 * (E - V_0));

    double complex psi_11, psi_12, psi_21, psi_22;

    psi_11 = cexp(I * q * (b - Delta)) / (4.0 * q * kappa)
        * ( cexp(I * kappa * Delta)   * (q + kappa) * (q + kappa)
          - cexp(-I * kappa * Delta)  * (q - kappa) * (q - kappa) );

    psi_22 = conj(psi_11);

    psi_12 = -I * cexp(I * q * (b - Delta)) / (2.0 * q * kappa)
        * (q * q - kappa * kappa) * csin(kappa * Delta);

    psi_21 = conj(psi_12);

    // quadratic determinantal eqn
    double complex u = -(psi_11 + psi_22);
    double complex v = (psi_11 * psi_22 - psi_12 * psi_21);
    double complex lam[2];

    lam[0] = (-u + csqrt(u * u - 4.0 * v)) / 2.0;
    lam[1] = (-u - csqrt(u * u - 4.0 * v)) / 2.0;

    for (int x = 0; x < 2; x++) {
        k[x] = clog(lam[x]) / (I * b);
    }
}

void gen_data(double b, double V_0, double Delta) {
    const double pi = 4.0 * atan(1.0);
    const double dE = 0.01;
    double E = dE;
    const int N = 3000;

    FILE *fp = fopen("data.csv", "w");
    if (!fp) {
        fp = fopen("/dev/null", "w");
    } 

    for (int y = 0; y < N; y++) {
        double complex q[2];
        solve(E, q, b, V_0, Delta);

        double rq = creal(q[0]);

        if (rq > 0 && rq < pi / b) {
            fprintf(fp, "%.15g,%.15g\n", rq, E);
            fprintf(fp, "%.15g,%.15g\n", -rq, E);
        }

        rq = creal(q[1]);

        if (rq > 0 && rq < pi / b) {
            fprintf(fp, "%.15g,%.15g\n", rq, E);
            fprintf(fp, "%.15g,%.15g\n", -rq, E);
        }

        E += dE;
    }

    fclose(fp);
}