#include "frame.h"
#include "optics.h"

#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(void) {
  const int width = 100, height = 100;
  const PointSpreadFunction psf = {.fwhm_pixels = 2.7, .moffat_beta = 4.5};
  const GeodesicTraceConfig trace = {.coordinate_time_step = 0.25,
                                     .max_steps = 100};
  const ObserverState observer = observer_fixed_at_origin();
  Star star = {
      .direction = {0.0, 0.0, -1.0}, .temperature_K = 7000.0, .amplitude = 1.0};
  StarCatalog catalog = {.stars = &star, .count = 1};
  SpacetimeSource spacetime = {0};
  FrameLensMesh mesh = {0};
  double *hdr = calloc((size_t)width * height * 3, sizeof *hdr);
  int result = 1;
  if (hdr == NULL || spacetime_create_minkowski(&spacetime, 10.0) ||
      frame_lens_mesh_build_coarse(&mesh, width, height, 20, 30.0) ||
      frame_lens_mesh_trace(&mesh, &spacetime, &observer, &trace))
    goto done;
  const size_t images =
      frame_splat_catalog(&mesh, &catalog, hdr, width, height, 100.0, &psf);
  if (images != 1 || hdr[3 * (50 * width + 50)] <= 0.0) {
    fputs("flat-space inverse lens-map regression failed\n", stderr);
    goto done;
  }
  /* Private HDR accumulation must preserve the serial splat result. */
  double *serial_hdr = calloc((size_t)width * height * 3, sizeof *serial_hdr);
  double *parallel_hdr = calloc((size_t)width * height * 3, sizeof *parallel_hdr);
  if (serial_hdr == NULL || parallel_hdr == NULL) {
    free(serial_hdr);
    free(parallel_hdr);
    goto done;
  }
  const int original_threads = omp_get_max_threads();
  omp_set_dynamic(0);
  omp_set_num_threads(1);
  const size_t serial_images = frame_splat_catalog(
      &mesh, &catalog, serial_hdr, width, height, 100.0, &psf);
  omp_set_num_threads(4);
  const size_t parallel_images = frame_splat_catalog(
      &mesh, &catalog, parallel_hdr, width, height, 100.0, &psf);
  omp_set_num_threads(original_threads);
  for (int value = 0; value < width * height * 3; ++value)
    if (fabs(serial_hdr[value] - parallel_hdr[value]) >
        1e-12 * fmax(1.0, fabs(serial_hdr[value]))) {
      fputs("parallel catalog splat regression failed\n", stderr);
      free(serial_hdr);
      free(parallel_hdr);
      goto done;
    }
  free(serial_hdr);
  free(parallel_hdr);
  if (serial_images != 1 || parallel_images != serial_images) {
    fputs("parallel catalog image-count regression failed\n", stderr);
    goto done;
  }
  const LinearRgb cool = blackbody_to_linear_rgb(3000.0);
  const LinearRgb hot = blackbody_to_linear_rgb(10000.0);
  if (!(cool.r > cool.b && hot.b > hot.r &&
        hot.r + hot.g + hot.b > cool.r + cool.g + cool.b)) {
    fputs("blackbody spectral-color regression failed\n", stderr);
    goto done;
  }
  /* The Moffat is flux-normalized and retains a measurable, continuous wing
   * beyond the former Gaussian's 3-sigma raster box. */
  memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr);
  splat_moffat(hdr, width, height, 50.5, 50.5,
               (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf);
  double moffat_flux = 0.0;
  for (int pixel = 0; pixel < width * height; ++pixel)
    moffat_flux += hdr[3 * pixel];
  if (fabs(moffat_flux - 1.0) > 0.01 ||
      hdr[3 * (50 * width + 62)] <= 0.0) {
    fputs("Moffat normalization or wing regression failed\n", stderr);
    goto done;
  }
  frame_draw_mesh(&mesh, hdr, width, height, 0.5, 0.5);
  if (hdr[3 * (10 * width + 20)] != 0.25) {
    fputs("mesh diagnostic overlay regression failed\n", stderr);
    goto done;
  }
  /* A fixed absolute edge tolerance used to make tiny source triangles claim
   * sources far outside their field. */
  FrameLensMesh fine_mesh = {0};
  Star fine_stars[2] = {{.direction = {0.0, 0.0, -1.0},
                         .temperature_K = 7000.0,
                         .amplitude = 1.0},
                        {.direction = {0.01, 0.0, -0.9999499987499375},
                         .temperature_K = 7000.0,
                         .amplitude = 1.0}};
  StarCatalog fine_catalog = {.stars = fine_stars, .count = 2};
  memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr);
  if (frame_lens_mesh_build_coarse(&fine_mesh, width, height, 1, 0.1) ||
      frame_lens_mesh_trace(&fine_mesh, &spacetime, &observer, &trace) ||
      frame_splat_catalog(&fine_mesh, &fine_catalog, hdr, width, height,
                          100.0, &psf) != 1) {
    fputs("fine source-triangle containment regression failed\n", stderr);
    frame_lens_mesh_destroy(&fine_mesh);
    goto done;
  }
  frame_lens_mesh_destroy(&fine_mesh);
  result = 0;
done:
  frame_lens_mesh_destroy(&mesh);
  spacetime_destroy(&spacetime);
  free(hdr);
  return result;
}