#include "fracture.h"
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_laguerre.h>
#include <gsl/gsl_vector.h>

void get_mean(uint_t len, uint32_t *data, double result[2]) {
  uint_t total = 0;
  double mean = 0;
  double mean_err = 0;

  for (uint_t i = 0; i < len; i++) {
    uint32_t datai = data[i];

    total += datai;
    mean += i * datai;
    mean_err += gsl_pow_2(i) * datai;
  }

  double meanf = mean / total;
  double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);

  result[0] = meanf;
  result[1] = meanf_err;
}

void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2],
             double result[2]) {
  uint_t total = 0;
  double mom = 0;
  double mom_err = 0;

  for (uint_t i = 0; i < len; i++) {
    uint32_t datai = data[i];
    double in = pow(fabs(((double)i) - mean[0]), n);

    total += datai;
    mom += in * datai;
    mom_err += gsl_pow_2(in) *
               datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0])));
  }

  double momf = mom / total;
  double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);

  result[0] = momf;
  result[1] = momf_err;
}

int main(int argc, char *argv[]) {
  uint_t nc = argc - 1;
  uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
  double *betas_c = (double *)malloc(nc * sizeof(double));
  double *vals_c1 = (double *)malloc(nc * sizeof(double));
  double *errors_c1 = (double *)malloc(nc * sizeof(double));
  double *vals_c2 = (double *)malloc(nc * sizeof(double));
  double *errors_c2 = (double *)malloc(nc * sizeof(double));
  double *vals_c3 = (double *)malloc(nc * sizeof(double));
  double *errors_c3 = (double *)malloc(nc * sizeof(double));
  double *vals_c4 = (double *)malloc(nc * sizeof(double));
  double *errors_c4 = (double *)malloc(nc * sizeof(double));

  char *out_filename = (char *)malloc(100 * sizeof(char));
  uint_t out_filename_len = 0;

  for (uint_t i = 0; i < nc; i++) {
    char *fn = argv[1 + i];
    char *buff = (char *)malloc(20 * sizeof(char));
    uint_t pos = 0;
    uint_t c = 0;
    while (c < 5) {
      if (fn[pos] == '_') {
        c++;
      }
      if (i == 0) {
        out_filename[pos] = fn[pos];
        out_filename_len++;
      }
      pos++;
    }
    c = 0;
    while (fn[pos] != '_') {
      buff[c] = fn[pos];
      pos++;
      c++;
    }
    buff[c] = '\0';
    Ls_c[i] = atoi(buff);
    c = 0;
    pos++;
    while (fn[pos] != '_') {
      buff[c] = fn[pos];
      pos++;
      c++;
    }
    buff[c] = '\0';
    betas_c[i] = atof(buff);
    free(buff);

    uint_t dist_len = 4 * pow(Ls_c[i], 2);
    uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
    FILE *dist_file = fopen(fn, "rb");
    fread(dist, sizeof(uint32_t), dist_len, dist_file);
    fclose(dist_file);
    {
      double mom1[2];
      get_mean(dist_len, dist, mom1);
      vals_c1[i] = mom1[0];
      errors_c1[i] = mom1[1];

      double mom2[2];
      get_mom(dist_len, 2, dist, mom1, mom2);
      vals_c2[i] = mom2[0];
      errors_c2[i] = mom2[1];

      double mom3[2];
      get_mom(dist_len, 3, dist, mom1, mom3);
      vals_c3[i] = mom3[0];
      errors_c3[i] = mom3[1];

      double mom4[2];
      get_mom(dist_len, 4, dist, mom1, mom4);
      vals_c4[i] = mom4[0];
      errors_c4[i] = mom4[1];
    }
    free(dist);
  }

  out_filename[out_filename_len - 1] = '.';
  out_filename[out_filename_len] = 't';
  out_filename[out_filename_len + 1] = 'x';
  out_filename[out_filename_len + 2] = 't';
  out_filename[out_filename_len + 3] = '\0';

  FILE *out_file = fopen(out_filename, "w");

  for (uint_t i = 0; i < nc; i++) {
    fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
            vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
            errors_c3[i], vals_c4[i], errors_c4[i]);
  }

  fclose(out_file);

  free(Ls_c);
  free(betas_c);
  free(vals_c1);
  free(errors_c1);
  free(vals_c2);
  free(errors_c2);
  free(vals_c3);
  free(errors_c3);

  return 0;
}