#define _XOPEN_SOURCE 600
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// Option to change numerical precision.
typedef int64_t int_t;
typedef double real_t;

// Simulation parameters: size, step count, and how often to save the state.
const int_t
    N = 1024,
    max_iteration = 4000,
    snapshot_freq = 10;

// Wave equation parameters, time step is derived from the space step.
const real_t
    c = 1.0,
    dx = 1.0;
real_t
    dt;

// Buffers for three time steps, indexed with 2 ghost points for the boundary.
real_t
    *buffers[3] = { NULL, NULL, NULL };

#define U_prv(i) buffers[0][(i) + 1]
#define U(i) buffers[1][(i) + 1]
#define U_nxt(i) buffers[2][(i) + 1]

// Save the present time step in a numbered file under 'data/'.
void domain_save(int_t step) {
  char filename[256];
  sprintf(filename, "data/%.5ld.dat", step);
  FILE *out = fopen(filename, "wb");
  fwrite(&U(0), sizeof(real_t), N, out);
  fclose(out);
}

// TASK: T1
// Set up our three buffers, fill two with an initial cosine wave,
// and set the time step.
void domain_initialize(void) {
  // allocate memory
  for (int i = 0; i < 3; i++)
    buffers[i] = malloc((N + 2) * sizeof(real_t));

  // apply initial condition
  real_t x = 0.0f;
  for (int i = 0; i < N; i++) {
    U_prv(i) = cos(2 * M_PI * x);
    U(i) = U_prv(i);
    U_nxt(i) = 0.0f;
    x += (real_t)dx / (real_t)N;
  }

  // set the time-step dt according to CFL
  dt = dx / c;
}

// TASK T2:
// Return the memory to the OS.
void domain_finalize(void) {
  for (int i = 0; i < 3; i++) {
    free(buffers[i]);
    buffers[i] = NULL;
  }
}

// TASK: T3
// Rotate the time step buffers.
void rotate_buffers(void) {
  real_t *temp = buffers[0];
  buffers[0] = buffers[1];
  buffers[1] = buffers[2];
  buffers[2] = temp;
}

// TASK: T4
// Derive step t+1 from steps t and t-1.
void time_step(void) {
  for (int i = 0; i < N; i++) {
    U_nxt(i) = -U_prv(i) + 2 * U(i) + dt * dt * c * c / (dx * dx) * (U(i - 1) + U(i + 1) - 2 * U(i));
  }
}

// TASK: T5
// Neumann (reflective) boundary condition.
void boundary(void) {
  U_prv(-1) = U_prv(1);
  U(-1) = U(1);
  U_nxt(-1) = U_nxt(1);
  U_prv(N) = U_prv(N - 2);
  U(N) = U(N - 2);
  U_nxt(N) = U_nxt(N - 2);
}

// TASK: T6
// Main time integration.
void simulate(void) {
  int_t iteration = 0;
  while (iteration < max_iteration) {
    boundary();
    time_step();
    rotate_buffers();
    domain_save(iteration / snapshot_freq);
    iteration++;
  }
}

int main(void) {
  domain_initialize();

  simulate();

  domain_finalize();
  exit(EXIT_SUCCESS);
}