class gumowskiMira {
  PImage pi;
  float pa, pb, newx, newy, oldx, oldy, logmaxd;
  int N = width;
  int maxdense = 0;
  int[][] density = new int[N][N];

  void construct() {
    //Produces the four variables to pass to the attractor
    pa = map(mouseX, 0, width, -1, 1);
    pb = map(mouseY, 0, height, 0.90, 1.001);
    oldx = width/2;
    oldy = height/2;
  }

  void populate(int s, boolean c) {
    //Populate array with density info with s number of samples
    int samples = s;
    boolean clear = c;
    if (clear) {
      for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
          density[i][j] = 0;
        }
      }
    }
    for (int i = 0; i < samples; i++) {
      oldx = random(0, width);
      oldy = random(0, height);
      for (int j = 0; j < 20000; j++) {
        newx = (pb * oldy + gumowskimira(oldx));
        newy = (gumowskimira(newx) - oldx);
        //Smoothie
        newx += random(-0.001, 0.001);
        newy += random(-0.001, 0.001);
        //If coordinates are within range, up density count at its position
        float plotx = (newx * 50) + width / 2;
        float ploty = (newy * 50) + height / 2;
        if ((plotx > 0) && (plotx < N) && (ploty > 0) && (ploty < N) ) {
          density[int(plotx)][int(ploty)] += 1;
        }
        oldx = newx;
        oldy = newy;
      }
    }
    //Put maximum density and its log()-value into variables
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
        if (density[i][j] > maxdense) {
          maxdense = density[i][j];
          logmaxd = log(maxdense);
        }
      }
    }
  }

  float gumowskimira(float x) {
    return pa * x + 2 * (1 - pa) * sq(x) * pow((1 + sq(x)), -2);
  }


  void updateloop() {
    stop = false;
    stepCounter = 0;
    pi = createImage(N, N, RGB);
  }

  void incrementalupdate() {
    //Loops the non-clearing update and plotting to produce low-noise render
    populate(16, false);
    plot(0, false);
    redraw();
  }

  void reparam() {
    //Fast reparametrization of variables
    construct();
    populate(1, true);
    plot(255, true);
  }

  PImage plot(int f, boolean c) {
    int factor = f;
    boolean clear = c;
    //Plot image from density array
    if (clear) {
      pi = createImage(N, N, RGB);
    }
    pi.loadPixels();
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
        if (density[i][j] > 0) {
          float myhue = map(logX(density[i][j]), 0, logX(maxdense), 0, 255);
          float mysat = map(logX(density[i][j]), 0, logX(maxdense), 0, 255);
          float mybright = map(logX(density[i][j]), 0, logX(maxdense), 0, 255) + factor;
          color newc = color(myhue, mysat, mybright);
          color oldc = pi.pixels[i * N + j];
          newc = blendColor(newc, oldc, LIGHTEST);
          pi.pixels[i * N + j] = newc;
        }
      }
    }
    pi.updatePixels();
    return pi;
  }
  
  float logX (float x) {
  return (log(x) / log(2));
}

}