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

typedef struct {
    double x, y, z;
} vector;

static inline vector mult_sv(double scalar, vector vec) {
    vector r = { scalar * vec.x, scalar * vec.y, scalar * vec.z };
    return r;
}

static inline vector add_vvv(vector v1, vector v2, vector v3) {
     vector r = {
         v1.x + v2.x + v3.x,
         v1.y + v2.y + v3.y,
         v1.z + v2.z + v3.z };
     return r;
}

static inline double dot_vv(vector v1, vector v2) {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
}

static inline int cube(int n) { return n*n*n; }

vector a = {1.0, 0.0, 0.0};
vector b = {0.0, 1.0, 0.0};
vector c = {0.5, 0.0, 0.8};
double meshStepLen = 10.0/64.0;
int nPerSide = 12;
int meshSideN = 100;

/*
 * The Fourier transform of δ(x - a, y - b, z - c) is
 *     exp(i(aX + bY + cZ))
 *     --------------------
 *          (2√2)(√π)³
 * using Mathematica's default definition of FT.
 */
double complex FTDD(vector realDDPos, vector reciprPos) {
    return cexp(I * dot_vv(realDDPos, reciprPos)) /
        (M_SQRT2 * M_2_SQRTPI * M_PI * M_PI);
}

int main() {
    // Make progress bar work.
    setvbuf (stdout, NULL, _IONBF, BUFSIZ);

    vector *directLattice = (vector*) malloc(cube(nPerSide)*sizeof(vector));
    {
        vector *directLatticeTmp = directLattice;
        for (int i = 0; i < nPerSide; i++) {
            for (int j = 0; j < nPerSide; j++) {
                for (int k = 0; k < nPerSide; k++) {
                    directLatticeTmp[0] = add_vvv(mult_sv(i, a), mult_sv(j, b), mult_sv(k, c));
                    directLatticeTmp++;
                }
            }
        }
    }

    double complex *reciprocalLattice = (double complex*) malloc(cube(meshSideN)*sizeof(double complex));
    for (int i = 0; i < cube(meshSideN); i++) {
        reciprocalLattice[i] = 0.0;
    }

    {
        double complex *reciprocalLatticeTmp = reciprocalLattice;

        for (int i = 0; i < meshSideN; i++) {
            putchar('.');
            for (int j = 0; j < meshSideN; j++) {
                for (int k = 0; k < meshSideN; k++) {
                    vector imagPoint = { i*meshStepLen, j*meshStepLen, k*meshStepLen };
                    for (int l = 0; l < cube(nPerSide); l++) {
                        reciprocalLatticeTmp[0] += FTDD(directLattice[l], imagPoint);
                    }
                    reciprocalLatticeTmp++;
                }
            }
        }
    }
    printf(" Created complex lattice\n");

    double *magData = (double*) malloc(cube(meshSideN)*sizeof(double));
    for (int i = 0; i < cube(meshSideN); i++) {
        magData[i] = cabs(reciprocalLattice[i]);
    }
    FILE *datafile = fopen("rlattice.bin", "wb");
    if (!datafile) {
        perror(0);
        goto free_rlattice;
    }
    fwrite(magData, sizeof(double)*cube(meshSideN), 1, datafile);
    fclose(datafile);

    printf("\nDone.\n");

free_rlattice:
    free(reciprocalLattice);
    free(directLattice);
    return 0;
}