import java.util.*;

int INNER_RADIUS = 7;
int OUTER_RADIUS = 3 * INNER_RADIUS;
double B1 = 0.278;
double B2 = 0.365;
double D1 = 0.267;
double D2 = 0.445;
double ALPHA_N = 0.028;
double ALPHA_M = 0.147;
int LOG_RES = 9;
int TRUE_RES = (1<<LOG_RES);
int[] COLOR_SHIFT = new int[]{0, 0, 0};
int[] COLOR_SCALE = new int[]{256, 256, 256};
int NR_OF_FIELDS = 2;

int current_field;
List<double[][]> fields = new ArrayList<double[][]>();
double[][] imaginary_field, M_re_buffer, M_im_buffer, N_re_buffer, N_im_buffer;
double[][] M_re, M_im, N_re, N_im;

void setup() {
  size(800, 800);

  current_field = 0;
  for (int i = 0; i < NR_OF_FIELDS; i++) {
    fields.add(zeromatrix());
  }
  imaginary_field = zeromatrix();
  M_re_buffer = zeromatrix();
  M_im_buffer = zeromatrix();
  N_re_buffer = zeromatrix();
  N_im_buffer = zeromatrix();


  //Precalculate multipliers for m,n
  BesselJ inner_bessel = BesselJ(INNER_RADIUS);
  BesselJ outer_bessel = BesselJ(OUTER_RADIUS);
  double inner_w = 1.0 / inner_bessel.w;
  double outer_w = 1.0 / (outer_bessel.w - inner_bessel.w);
  M_re = inner_bessel.re;
  M_im = inner_bessel.im;
  N_re = outer_bessel.re;
  N_im = outer_bessel.im;
  for (int i=0; i < TRUE_RES; ++i) {
    for (int j=0; j < TRUE_RES; ++j) {
      N_re[i][j] = outer_w * (N_re[i][j] - M_re[i][j]);
      N_im[i][j] = outer_w * (N_im[i][j] - M_im[i][j]);
      M_re[i][j] *= inner_w;
      M_im[i][j] *= inner_w;
    }
  }

  background(0);
  frameRate(30);
  init_field(0);
  speckle_field(2000, 1);
}

void draw() {
  perform_calculation_step();

  double[][] cur_field = fields.get(current_field);

  if (mousePressed){
    int v = (int) (((double) mouseX / (double) width) * (double) TRUE_RES);
    int u = (int) (((double) mouseY / (double) height) * (double) TRUE_RES);
    for (int x = -OUTER_RADIUS/2; x < OUTER_RADIUS/2; x++) {
      for (int y = -OUTER_RADIUS/2; y < OUTER_RADIUS/2; y++) {
        cur_field[u+x][v+y] = 1;
      }
    }
    for (int x = -INNER_RADIUS/2; x < INNER_RADIUS/2; x++) {
      for (int y = -INNER_RADIUS/2; y < INNER_RADIUS/2; y++) {
        cur_field[u+x][v+y] = 0;
      }
    }
  }

  int image_ptr = 0;
  PImage img = createImage(TRUE_RES, TRUE_RES, RGB);
  img.loadPixels();
  for (int i = 0; i < TRUE_RES; i++) {
    for (int j = 0; j < TRUE_RES; j++) {
      double s = cur_field[i][j];
      int[] clr = new int[]{255, 255, 255};
      for (int k = 0; k < 3; k++) {
        clr[k] = (int) Math.max(0, Math.min(255, Math.floor(COLOR_SHIFT[k] + COLOR_SCALE[k]*s)));
      }
      img.pixels[image_ptr] = color(clr[0], clr[1], clr[2]);
      image_ptr++;
    }
  }
  img.updatePixels();

  //Blit to screen
  img.resize(width, height);
  image(img, 0, 0);
}

void init_field(double x) {
  double[][] cur_field = fields.get(current_field);
  for (int i = 0; i < TRUE_RES; i++) {
    for (int j = 0; j < TRUE_RES; j++) {
      cur_field[i][j] = x;
    }
  }
}

void speckle_field(int count, double intensity) {
  double[][] cur_field = fields.get(current_field);
  for (int i = 0; i < count; i++) {
    int u = (int) Math.floor(Math.random() * (TRUE_RES - INNER_RADIUS));
    int v = (int) Math.floor(Math.random() * (TRUE_RES - INNER_RADIUS));
    for (int x = 0; x < INNER_RADIUS; x++) {
      for (int y = 0; y < INNER_RADIUS; y++) {
        cur_field[u+x][v+y] = intensity;
      }
    }
  }
}

void perform_calculation_step() {
  //Read in fields
  double[][] cur_field = fields.get(current_field);
  current_field = (current_field + 1) % fields.size();
  double[][] next_field = fields.get(current_field);
  //Clear extra imaginary field
  for (int i = 0; i < TRUE_RES; i++) {
    for (int j = 0; j < TRUE_RES; j++) {
      imaginary_field[i][j] = 0.0;
    }
  }
  //Compute m,n fields
  fft2(true, LOG_RES, cur_field, imaginary_field);
  field_multiply(cur_field, imaginary_field, M_re, M_im, M_re_buffer, M_im_buffer);
  fft2(false, LOG_RES, M_re_buffer, M_im_buffer);
  field_multiply(cur_field, imaginary_field, N_re, N_im, N_re_buffer, N_im_buffer);
  fft2(false, LOG_RES, N_re_buffer, N_im_buffer);
  //Step s
  for (int i=0; i<next_field.length; ++i) {
    for (int j=0; j<next_field.length; ++j) {
      next_field[i][j] = S(N_re_buffer[i][j], M_re_buffer[i][j]);
    }
  }
}

class BesselJ {
  double[][] re;
  double[][] im;
  double w;
  BesselJ(double[][] re, double[][] im, double w) {
    this.re = re;
    this.im = im;
    this.w = w;
  }
}

BesselJ BesselJ(double radius) {
  double[][] field = zeromatrix();
  double weight = 0.0;
  for (int i = 0; i < field.length; i++) {
    for (int j = 0; j < field.length; j++) {
      double ii = ((i + field.length/2) % field.length) - field.length/2;
      double jj = ((j + field.length/2) % field.length) - field.length/2;
      double r = Math.sqrt(ii*ii + jj*jj) - radius;
      double v = 1.0 / (1.0 + Math.exp(LOG_RES * r));
      weight += v;
      field[i][j] = v;
    }
  }
  double[][] imag_field = zeromatrix();
  fft2(true, LOG_RES, field, imag_field);
  return new BesselJ(field, imag_field, weight);
}

double sigma(double x, double a, double alpha) {
  return 1.0 / (1.0 + Math.exp(-4.0/alpha * (x - a)));
}

double sigma_2(double x, double a, double b) {
  return sigma(x, a, ALPHA_N) * (1.0 - sigma(x, b, ALPHA_N));
}

double lerp(double a, double b, double t) {
  return (1.0-t)*a + t*b;
}

double S(double n, double m) {
  double alive = sigma(m, 0.5, ALPHA_M);
  return sigma_2(n, lerp(B1, D1, alive), lerp(B2, D2, alive));
}

void field_multiply(double[][] a_r, double[][] a_i,
  double[][] b_r, double[][] b_i,
  double[][] c_r, double[][] c_i) {
  for (int i = 0; i < TRUE_RES; i++) {
    double[] Ar = a_r[i];
    double[] Ai = a_i[i];
    double[] Br = b_r[i];
    double[] Bi = b_i[i];
    double[] Cr = c_r[i];
    double[] Ci = c_i[i];
    for (int j = 0; j < TRUE_RES; j++) {
      double a = Ar[j];
      double b = Ai[j];
      double c = Br[j];
      double d = Bi[j];
      double t = a * (c + d);
      Cr[j] = t - d*(a+b);
      Ci[j] = t + c*(b-a);
    }
  }
}

void fft(boolean dir, int m, double[] x, double[] y) {
  int nn, i, i1, j, k, i2, l, l1, l2;
  double c1, c2, tx, ty, t1, t2, u1, u2, z;
  nn = x.length;
  /* Do the bit reversal */
  i2 = nn >> 1;
  j = 0;
  for (i=0; i<nn-1; i++) {
    if (i < j) {
      tx = x[i];
      ty = y[i];
      x[i] = x[j];
      y[i] = y[j];
      x[j] = tx;
      y[j] = ty;
    }
    k = i2;
    while (k <= j) {
      j -= k;
      k >>= 1;
    }
    j += k;
  }
  /* Compute the FFT */
  c1 = -1.0;
  c2 = 0.0;
  l2 = 1;
  for (l = 0; l < m; l++) {
    l1 = l2;
    l2 <<= 1;
    u1 = 1.0;
    u2 = 0.0;
    for (j=0; j<l1; j++) {
      for (i=j; i<nn; i+=l2) {
        i1 = i + l1;
        t1 = u1 * x[i1] - u2 * y[i1];
        t2 = u1 * y[i1] + u2 * x[i1];
        x[i1] = x[i] - t1;
        y[i1] = y[i] - t2;
        x[i] += t1;
        y[i] += t2;
      }
      z =  u1 * c1 - u2 * c2;
      u2 = u1 * c2 + u2 * c1;
      u1 = z;
    }
    c2 = Math.sqrt((1.0 - c1) / 2.0);
    if (dir)
      c2 = -c2;
    c1 = Math.sqrt((1.0 + c1) / 2.0);
  }
  /* Scaling for forward transform */
  if (!dir) {
    double scale_f = 1.0 / nn;
    for (i=0; i<nn; i++) {
      x[i] *= scale_f;
      y[i] *= scale_f;
    }
  }
}

void fft2(boolean dir, int m, double[][] x, double[][] y) {
  /* In place 2D FFT */
  for (int i = 0; i < x.length; ++i) {
    fft(dir, m, x[i], y[i]);
  }
  for (int i=0; i<x.length; ++i) {
    for (int j=0; j<i; ++j) {
      double t = x[i][j];
      x[i][j] = x[j][i];
      x[j][i] = t;
    }
  }
  for (int i=0; i<y.length; ++i) {
    for (int j=0; j<i; ++j) {
      double t = y[i][j];
      y[i][j] = y[j][i];
      y[j][i] = t;
    }
  }
  for (int i=0; i < x.length; ++i) {
    fft(dir, m, x[i], y[i]);
  }
}

double[][] zeromatrix() {
  return matrix(TRUE_RES, TRUE_RES, 0.0);
}

static double[][] matrix(int rows, int cols, double value) {
  double[][] matrix = new double[rows][cols];
  for (int row = 0; row < rows; row++)
    for (int col = 0; col < cols; col++)
      matrix[row][col] = value;
  return matrix;
}