/* Copyright ©2015-2024 Evan Pretti.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 *     1.  Redistributions of source code must retain the above copyright
 *         notice, this list of conditions and the following disclaimer.
 *     2.  Redistributions in binary form must reproduce the above copyright
 *         notice, this list of conditions and the following disclaimer in the
 *         documentation and/or other materials provided with the distribution.
 *
 * This software is provided by the author and contributors "as is" and any
 * express or implied warranties, including, but not limited to, the implied
 * warranties of merchantability and fitness for a particular purpose are
 * disclaimed.  In no event shall the author or contributors be liable for any
 * direct, indirect, incidental, special, exemplary, or consequential damages
 * (including, but not limited to, procurement of substitute goods or services;
 * loss of use, data, or profits; or business interruption) however caused and
 * on any theory of liability, whether in contract, strict liability, or tort
 * (including negligence or otherwise) arising in any way out of the use of this
 * software, even if advised of the possibility of such damage. */

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define π 3.1415926535897932
typedef struct xyi_t {
    size_t x;
    size_t y;
} xyi_t;
typedef struct xyd_t {
    double x;
    double y;
} xyd_t;
double sq(double x);
xyi_t xyi(size_t x, size_t y);
xyd_t xyd(double x, double y);
uint32_t color(double z);
double evaluate(xyd_t r, double λ, xyd_t* sources, size_t source_count);
void render(uint32_t* data, xyi_t size, xyd_t r1, xyd_t r2, double λ, xyd_t* sources, size_t source_count);
int main();
double sq(double x) {
    return x * x;
}
xyi_t xyi(size_t x, size_t y) {
    xyi_t result;
    result.x = x;
    result.y = y;
    return result;
}
xyd_t xyd(double x, double y) {
    xyd_t result;
    result.x = x;
    result.y = y;
    return result;
}
uint32_t color(double z) {
    uint8_t i;
    if(z < 0.0) {
        z = 0.0;
    }
    if(z > 1.0) {
        z = 1.0;
    }
    if(z < 0.25) {
        z *= 4.0;
        i = (int)(z * 255.0);
        return 0 | i << 8 | 255 << 16;
    }
    if(z < 0.5) {
        z -= 0.25;
        z *= 4.0;
        i = (int)(z * 255.0);
        return 0 | 255 << 8 | (255 - i) << 16;
    }
    if(z < 0.75) {
        z -= 0.5;
        z *= 4.0;
        i = (int)(z * 255.0);
        return i | 255 << 8 | 0 << 16;
    } else {
        z -= 0.75;
        z *= 4.0;
        i = (int)(z * 255.0);
        return 255 | (255 - i) << 8 | 0 << 16;
    }
}
double evaluate(xyd_t r, double λ, xyd_t* sources, size_t source_count) {
    size_t source_index;
    xyd_t source;
    double phase;
    double real = 0.0;
    double imaginary = 0.0;
    for(source_index = 0; source_index < source_count; source_index++) {
        source = sources[source_index];
        phase = 2 * π * sqrt(sq(r.x - source.x) + sq(r.y - source.y)) / λ;
        real += cos(phase);
        imaginary += sin(phase);
    }
    return sqrt(sq(real / source_count) + sq(imaginary / source_count));
}
void render(uint32_t* data, xyi_t size, xyd_t r1, xyd_t r2, double λ, xyd_t* sources, size_t source_count) {
    xyi_t image_point;
    xyd_t field_point;
    size_t pixel_count = size.x * size.y;
    for(size_t pixel_index = 0; pixel_index < pixel_count; pixel_index++) {
        if(pixel_index % 10000 == 0) {
            printf("%zu\n", pixel_index);
        }
        image_point = xyi(pixel_index % size.x, pixel_index / size.x);
        field_point = xyd((((double)image_point.x / size.x) * (r2.x - r1.x)) + r1.x, ((1.0 - ((double)image_point.y / size.y)) * (r2.y - r1.y)) + r1.y);
        data[pixel_index] = color(evaluate(field_point, λ, sources, source_count)) | 0xFF000000;
    }
}
int main() {
    size_t pixel_count;
    uint32_t* data;
    FILE* file;
    uint64_t width;
    uint64_t height;
#define SLITS 16
#define SAMPLES 128
#define SPACING 4.0
#define WIDTH 2.0
#define IMAGE 2048.0
#define PATWIDTH ((SPACING * (SLITS - 1)) + WIDTH)
    xyi_t size = xyi(1024, 1024);
    xyd_t r1 = xyd((PATWIDTH - IMAGE) / 2, IMAGE * 0);
    xyd_t r2 = xyd((IMAGE + PATWIDTH) / 2, IMAGE * 1);
    double λ = 1.0;
    xyd_t sources[SLITS * SAMPLES];
    size_t source_count = SLITS * SAMPLES;
    size_t source_index;
    char* path = "interference.dat";
    for(size_t slit = 0; slit < SLITS; slit++) {
        for(size_t sample = 0; sample < SAMPLES; sample++) {
            sources[(slit * SAMPLES) + sample] = xyd((SPACING * slit) + sample * (WIDTH / SAMPLES), 0.0);
        }
    }
    pixel_count = size.x * size.y;
    data = calloc(pixel_count, sizeof(uint32_t));
    render(data, size, r1, r2, λ, sources, source_count);
    file = fopen(path, "wb");
    width = (uint64_t)(size.x);
    height = (uint64_t)(size.y);
    fwrite(&width, sizeof(uint64_t), 1, file);
    fwrite(&height, sizeof(uint64_t), 1, file);
    fwrite(data, sizeof(uint32_t), pixel_count, file);
    fclose(file);
    free(data);
}