/*
 * Sandpile Model (Bak-Tang-Wiesenfeld)
 * C Code (eggx library) from Y.Okamoto ported to Processing (Java)
 * Added delay function for the destruction process (applied after drawing post-collapse)
 Original code is as follows,
[sand_pl.bas]***********************************************[Ver 0.00]
  S.O.C.
 model for earthquakes(Sand_pile model)
        algo.
 from T.Bak & C.Tang(1989),K.Ito & M.Matsuzaki(1990)
   start from random state(0-3)

                    special thanks to Dr.M.Matsuzaki
        prog.by Y.Okamoto 1997.03/17-06/10
*****************************************************************
 modified and ported by Gemini on 25 Dec.2025 
*/

// --------------------------------------------------------------------
// 1. Global Constants
// --------------------------------------------------------------------
final int T_END = 10000; // Maximum number of iterations
final int T_SKIP = 1000; // Number of iterations to skip for log output

final int X00 = 50;      // Display origin X
final int Y00 = 50;      // Display origin Y (Processing Y-axis is downward)
final int IMAX = 128;    // Grid size X
final int JMAX = 128;    // Grid size Y

final int DX = 10;       // Cell interval X
final int DY = 10;       // Cell interval Y

// [Fix 1] Expand window size: Drawing grid (128*10 + 50*2 = 1380) + Info panel (600) = 1980
final int GRID_DRAW_WIDTH = IMAX * DX + 2 * X00;
final int INFO_WIDTH = 600;      // Expanded for graph drawing
final int WINWIDTH = GRID_DRAW_WIDTH + INFO_WIDTH; 
final int WINHEIGHT = 1200;

final int CRITICAL_H = 4;  // Critical height (topples if h >= 4)
final int CMAX = 4096;     // Guideline for maximum cluster size

final float R = 4.0;       // Particle radius (display size)

// [Added] Graph constants
final float LOG_BIN_WIDTH = 0.1f; // Bin width on Log10 scale (for log-binning)

// --------------------------------------------------------------------
// 2. Global Variables
// --------------------------------------------------------------------
// Sand height array. +2 for boundary conditions (i=0, IMAX+1, etc.)
int[][] h = new int[IMAX + 2][JMAX + 2]; 
// Toppling flag array (records if a cell collapsed)
int[][] f = new int[IMAX + 2][JMAX + 2];  

int t = 0;             // Time/Iteration counter (drop count)
int c, cc;             // Toppled cells (c: calculated, cc: during drawing)
boolean flg;           // Avalanche chain flag
int[] n = new int[CMAX + 1];  // Frequency of avalanche size (c) (Calculation-based)
int[] nn = new int[CMAX + 1]; // Frequency of avalanche size (cc) (Drawing-based)

// --------------------------------------------------------------------
// 3. setup(): Initialization and Initial Drawing
// --------------------------------------------------------------------
void setup() {
  // [Fix 2] Set window size according to global variables
  size(2200, 1400); 
  background(255);         // Set background to white
  smooth();                // Anti-aliasing
  frameRate(60);           // Normal drawing speed
  
  // Font settings
  textAlign(LEFT, TOP);
  textSize(14);

  // Initialize grid values (active area: i, j = 1 to IMAX, JMAX)
  for (int i = 1; i <= IMAX; i++) {
    for (int j = 1; j <= JMAX; j++) {
      // Initialize with integers 0-3 (0 <= h < 4)
      h[i][j] = floor(sqrt(16.0 * random(1))); 
      f[i][j] = 0; // Initialize toppling flags
      
      // Initial drawing
      drawParticle(i, j, h[i][j], false);
    }
  }
  
  // Initial display of statistics panel
  drawStatistics();
  
  println("Simulation Start. Outputting data for t > " + T_SKIP);
  println("t      c      cc  (Time, Fluctuation Count, Drawing Fluctuation Count)");
}

// --------------------------------------------------------------------
// 4. draw(): Main Loop (Time Iteration)
// --------------------------------------------------------------------
void draw() {
  if (t >= T_END) {
    // Process after time ends
    drawLogLogPlot(); // Final graph update
    noLoop();         // Stop draw() loop
    
    // Output avalanche frequency data (cc, nn[cc])
    println("\n--- Avalanch Frequency (cc, nn[cc]) ---");
    for (int cc_val = 1; cc_val <= CMAX; cc_val++) {
      if (nn[cc_val] > 0) {
        println(nf(cc_val, 6) + " " + nf(nn[cc_val], 6));
      }
    }
    return;
  }

  t = t + 1;
  flg = true; // Avalanche chain reaction flag
  c = 0;      // Toppled cells (Calculation-based)

  // Clear f[][] flags
  for (int i = 0; i <= IMAX + 1; i++) {
    for (int j = 0; j <= JMAX + 1; j++) {
      f[i][j] = 0;
    }
  }

  // 1. Add one grain of sand at a random position
  int drop_i = floor(random(1, IMAX + 1));
  int drop_j = floor(random(1, JMAX + 1));
  h[drop_i][drop_j] = h[drop_i][drop_j] + 1;

  // 2. Stabilization (Avalanche Chain) loop (Calculation only)
  while (flg) {
    flg = false;  
    boolean chain_reaction = false;

    for (int i = 1; i <= IMAX; i++) {
      for (int j = 1; j <= JMAX; j++) {
        if (h[i][j] >= CRITICAL_H) { // Critical value 4
          chain_reaction = true;
          // Toppling process
          h[i][j] = h[i][j] - 4;
          h[i - 1][j] = h[i - 1][j] + 1;
          h[i + 1][j] = h[i + 1][j] + 1;
          h[i][j - 1] = h[i][j - 1] + 1;
          h[i][j + 1] = h[i][j + 1] + 1;
          f[i][j] = 1; // Record toppled cell
          c = c + 1;   // Count toppled cells (Calculation-based)
        }
      }
    }
    if (chain_reaction) flg = true;
  } 

  // 3. Boundary Conditions (h[boundary] = 0)
  for (int i = 0; i <= IMAX + 1; i++) {
    h[i][JMAX + 1] = 0;
    h[i][0] = 0;     
  }
  for (int j = 0; j <= JMAX + 1; j++) {
    h[IMAX + 1][j] = 0;
    h[0][j] = 0;     
  }
  
  // 4. Final drawing and count cc
  cc = 0; 
  // Clear the drawing grid
  fill(255);
  rect(X00 - R, Y00 - R, IMAX * DX + R * 2, JMAX * DY + R * 2);

  for (int i = 1; i <= IMAX; i++) {
    for (int j = 1; j <= JMAX; j++) {
      boolean isFalling = (f[i][j] == 1);
      drawParticle(i, j, h[i][j], isFalling);
      if (isFalling) {
        cc = cc + 1; // Count toppled cells (Drawing-based)
        f[i][j] = 0;  // Reset drawing flag
      }
    }
  }
  
  // 5. Data recording and delay insertion
  if (cc > 0) {
    // Insert delay only if an avalanche (cc > 0) occurred
    delay(50); // 50ms delay (highlights the result after the chain ends)
    
    if (t > T_SKIP) {
      // Log outputs both calculation-based (c) and drawing-based (cc)
      println(nf(t, 6) + " " + nf(c, 6) + " " + nf(cc, 6));
      // Record only drawing-based (cc) in frequency array
      if (cc <= CMAX) nn[cc]++;
    }
  }
  
  // 6. Update statistics display
  drawStatistics();
}

// --------------------------------------------------------------------
// 5. drawParticle(): Particle Drawing Function
// --------------------------------------------------------------------
void drawParticle(int i, int j, int height, boolean isFalling) {
  float xx = i * DX + X00;
  float yy = j * DY + Y00;

  color c;
  if (isFalling) {
    c = color(55, 0, 255); // Highlight toppling cell
  } else {
    // Color coding based on sand height (h[i][j])
    switch (height) {
      case 0: c = color(255); break; // White
      case 1: c = color(180, 255,255); break; // Cyan/Green
      case 2: c = color(255, 205, 75); break; // Yellow
      case 3: c = color(255, 125, 75); break; // Orange
      case 4: c = color(0, 0, 255); break; // Purple/Blue
      default: c = color(255); break; // White (for h > 4 immediately after calculation)
    }
  }
  
  fill(c);
  noStroke();
  ellipse(xx, yy, R * 2, R * 2);
}


// --------------------------------------------------------------------
// 6. drawStatistics(): Statistics Info Drawing Function
// --------------------------------------------------------------------
void drawStatistics() {
  // Clear information area / Gray background
  fill(240);
  rect(GRID_DRAW_WIDTH, 0, INFO_WIDTH, WINHEIGHT);
  
  int startX = GRID_DRAW_WIDTH + 20;
  int currentY = 50;
  
  fill(0); // Set text color to black
  textSize(16);
  text("--- SOC Sandpile Stats ---", startX, currentY);
  currentY += 30;
  
  textSize(14);
  // Drop count (Time t)
  text("Time (t):", startX, currentY);
  text(nf(t, 6) + " / " + T_END, startX, currentY + 20);
  currentY += 50;

  // Avalanche Size (c: Calculation-based)
  text("Avalanche Size (C):", startX, currentY);
  if (c > 0) fill(255, 0, 0); // Highlight if avalanche occurs
  text(nf(c, 4), startX, currentY + 20);
  currentY += 50;

  // Drawn Avalanche Size (cc: Drawing-based)
  fill(0);
  text("Drawn Size (CC):", startX, currentY);
  if (cc > 0) fill(55, 0, 255); // Highlight with purple
  text(nf(cc, 4), startX, currentY + 20);
  currentY += 50;

  // Log Output Status
  fill(0);
  text("Log Output > T_SKIP:", startX, currentY);
  fill(t > T_SKIP ? color(0, 150, 0) : color(150, 0, 0));
  text(t > T_SKIP ? "ON" : "OFF (t=" + T_SKIP + ")", startX, currentY + 20);
  currentY += 50;

  // Draw plot if simulation is finished
  if (t >= T_END) drawLogLogPlot();
}

// --------------------------------------------------------------------
// 7. drawLogLogPlot(): Log-Log Plot Function (Log-binning & bin width correction)
// --------------------------------------------------------------------
void drawLogLogPlot() {
  fill(255); 
  rect(GRID_DRAW_WIDTH, 0, INFO_WIDTH, WINHEIGHT);

  // [Step 1] Determine maximum value
  int max_size_cc = 0;
  for (int size = 1; size <= CMAX; size++) {
    if (nn[size] > 0) {
      if (size > max_size_cc) max_size_cc = size;
    }
  }

  if (max_size_cc <= 1) {
    fill(0);
    text("No data to plot (or too small).", GRID_DRAW_WIDTH + 50, 50);
    return;
  }
  
  // [Step 2] Calculate Log-bin array N_bin
  float log10_max_size = (float) (Math.log10(max_size_cc) + LOG_BIN_WIDTH);
  int num_bins = (int) Math.ceil(log10_max_size / LOG_BIN_WIDTH);
  
  float[] N_bin = new float[num_bins]; // Total frequency per bin
  float min_log_density = Float.MAX_VALUE;

  for (int S = 1; S <= max_size_cc; S++) {
      if (nn[S] > 0) {
          int bin_index = (int) Math.floor(Math.log10(S) / LOG_BIN_WIDTH);
          if (bin_index < num_bins) N_bin[bin_index] += nn[S];
      }
  }
  
  // [Step 3] Calculate Cumulative Frequency N_cum_bin
  float[] N_cum_bin = new float[num_bins];
  float cumulative_sum = 0;
  
  // Accumulate in reverse from the largest bin (CCDF)
  for (int i = num_bins - 1; i >= 0; i--) {
      cumulative_sum += N_bin[i];
      N_cum_bin[i] = cumulative_sum;
  }
  
  float total_count_float = cumulative_sum; // Sum of N_bin = Y-axis max guideline

  // [Step 4] Calculate Minimum Log(Density) for discrete values
  for (int bin_index = 0; bin_index < num_bins; bin_index++) {
    if (N_bin[bin_index] > 0) {
      float count = N_bin[bin_index];
      float S_min = (float) Math.pow(10, bin_index * LOG_BIN_WIDTH);
      float S_max = (float) Math.pow(10, (bin_index + 1) * LOG_BIN_WIDTH);
      float delta_S = S_max - S_min;
      
      // Bin width correction (Density = Frequency / Bin width)
      float density = count / delta_S;
      float log_density = (float) log(density);
      if (log_density < min_log_density) min_log_density = log_density;
    }
  }
  
  // Plot area definition
  final int MARGIN = 40;
  final int PLOT_W = (INFO_WIDTH - MARGIN * 2); // Approximately double width
  final int PLOT_H = 600;
  final int PLOT_X = GRID_DRAW_WIDTH + MARGIN;
  final int PLOT_Y = WINHEIGHT - PLOT_H - MARGIN;

  // Log scale calculations
  float log_max_size = log(max_size_cc);
  float log_max_y = log(total_count_float); // Max value of CCDF (total events)
  float log_min_size = log(1);
  // Set Y-axis min (use smaller of log(1) or min log-density with padding)
  float log_min_y = min(log(1), min_log_density); 
  log_min_y -= 0.5f;

  // Draw Axes and Labels
  stroke(0);
  fill(0);
  textSize(10);
  noFill();
  rect(PLOT_X, PLOT_Y, PLOT_W, PLOT_H);
  
  textAlign(CENTER, BOTTOM);
  text("Cluster Size S (log_e)", PLOT_X + PLOT_W / 2, PLOT_Y + PLOT_H + 20);
  
  textAlign(CENTER, CENTER);
  pushMatrix();
  translate(PLOT_X - 15, PLOT_Y + PLOT_H / 2);
  rotate(-HALF_PI);
  text("Density / Cumulative Count (log_e)", 0, 0);
  popMatrix();

  fill(0);
  textAlign(LEFT, TOP);
  textSize(14);
  text("Size Distribution (Log Binned Discrete Density & CCDF)", PLOT_X, PLOT_Y - 25);

  // X-axis ticks (10^N)
  for (int i = 0; i < 4; i++) {
      float val = (float)Math.pow(10, i);
      if (log(val) > log_max_size * 1.05f) break; 
      float x_tick = map(log(val), log_min_size, log_max_size, PLOT_X, PLOT_X + PLOT_W);
      line(x_tick, PLOT_Y + PLOT_H, x_tick, PLOT_Y + PLOT_H + 5);
      textAlign(CENTER, TOP);
      text("10^" + i, x_tick, PLOT_Y + PLOT_H + 5);
  }
  
  // Y-axis ticks (10^N)
  float log10_max_y = (float) (Math.log10(total_count_float) + 1.0);
  float log10_min_y = (float) Math.log10(Math.exp(log_min_y));
  
  for (int i = (int)Math.ceil(log10_min_y); i <= (int)Math.floor(log10_max_y); i++) {
      float val = (float)Math.pow(10, i);
      float log_val = log(val); 
      if (log_val < log_min_y - 0.1f || log_val > log_max_y + 0.1f) continue;
      
      float y_tick = map(log_val, log_min_y, log_max_y, PLOT_Y + PLOT_H, PLOT_Y);
      line(PLOT_X, y_tick, PLOT_X - 5, y_tick);
      textAlign(RIGHT, CENTER);
      text("10^" + i, PLOT_X - 10, y_tick);
  }

  // Data Plotting
  // 1. [Discrete Plot] - Blue (Log-binned density, bin-width corrected)
  for (int bin_index = 0; bin_index < num_bins; bin_index++) {
    if (N_bin[bin_index] > 0) {
      float count = N_bin[bin_index];
      float log10_S_center = (bin_index * LOG_BIN_WIDTH) + (LOG_BIN_WIDTH / 2.0f);
      float S_center = (float) Math.pow(10, log10_S_center);
      float S_min = (float) Math.pow(10, bin_index * LOG_BIN_WIDTH);
      float S_max = (float) Math.pow(10, (bin_index + 1) * LOG_BIN_WIDTH);
      float density = count / (S_max - S_min); // Bin width correction
      
      float x_plot = map((float)log(S_center), log_min_size, log_max_size, PLOT_X, PLOT_X + PLOT_W);
      float y_plot = map((float)log(density), log_min_y, log_max_y, PLOT_Y + PLOT_H, PLOT_Y);
      fill(0, 0, 255); // Blue
      noStroke();
      ellipse(x_plot, y_plot, 4, 4);
    }
  }
  
  // 2. [Cumulative Plot] - Red (Using log-binned cumulative values)
  for (int bin_index = 0; bin_index < num_bins; bin_index++) {
      if (N_cum_bin[bin_index] > 0) {
          float log10_S_center = (bin_index * LOG_BIN_WIDTH) + (LOG_BIN_WIDTH / 2.0f);
          float S_center = (float) Math.pow(10, log10_S_center);
          float x_plot = map(log(S_center), log_min_size, log_max_size, PLOT_X, PLOT_X + PLOT_W);
          float y_plot_cum = map(log(N_cum_bin[bin_index]), log_min_y, log_max_y, PLOT_Y + PLOT_H, PLOT_Y);
          fill(255, 0, 0); // Red
          noStroke();
          rect(x_plot - 2, y_plot_cum - 2, 4, 4); // Plot with square
      }
  }

  // Legend
  textSize(12);
  fill(0, 0, 255);
  rect(PLOT_X + PLOT_W - 150, PLOT_Y + 10, 8, 8);
  fill(0);
  textAlign(LEFT, TOP);
  text("Discrete (Density, Bin Corrected)", PLOT_X + PLOT_W - 140, PLOT_Y + 8);
  fill(255, 0, 0);
  rect(PLOT_X + PLOT_W - 150, PLOT_Y + 30, 8, 8);
  fill(0);
  text("Cumulative (CCDF)", PLOT_X + PLOT_W - 140, PLOT_Y + 28);
}