// color ---------------------------------------------------------------------------------------------------------------------


color dkBrownColor = color(77,47,31);
color medBrownColor = color(125,88,64);
color ltBrownColor = color(185,142,92);
color medGreenColor = color(66,90,33);
color ltGreenColor = color(122,140,66);
color ltltGreenColor = color(169,190,99);
color dkBlueColor = color(14,30,60);
color medBlueColor = color(34,65,105);
color ltBlueColor = color(54,100,150);
color ltltBlueColor = color(74,135,195);
color dkPurpleColor = color(57,24,51);
color medPurpleColor = color(80,53,73);
color dkRedColor = color(112,15,0);
color medRedColor = color(146,54,38);
color ltRedColor = color(180,90,76);
color ltltRedColor = color(214,126,114);
color orangeColor = color(186,92,38);
color yellowColor = color(191,146,54);
color offWhiteColor = color(255,255,243);
color dkGrayColor = color(0);
color medGrayColor = color(102);
color ltGrayColor = color(204);


color shift(color col, float r) {
  if (r>0) {
    return lerpColor(col,color(255,255,255,alpha(col)),r);
  } else {
    return lerpColor(col,color(0,0,0,alpha(col)),-r);
  }
}

color lighten(color col) {return shift(col,0.25);}

color darken(color col) {return shift(col,-0.25);}

color randomshift(color col) {
  return shift(col, random(-0.5,0.5));
}

color transparency(color col, float f) {
  return color(red(col),green(col),blue(col),f*alpha(col));
}

color desaturate(color col) {
  return color(brightness(col),brightness(col),brightness(col),alpha(col));
}

color[] colormapcat(int N, color col1, color col2) {return colormapcat(N, new color[] {col1, col2});}
color[] colormapcat(int N, color col1, color col2, color col3) {return colormapcat(N, new color[] {col1, col2, col3});}
color[] colormapcat(int N, color col1, color col2, color col3, color col4) {return colormapcat(N, new color[] {col1, col2, col3, col4});}
color[] colormapcat(int N, color col1, color col2, color col3, color col4, color col5) {return colormapcat(N, new color[] {col1, col2, col3, col4, col5});}
color[] colormapcat(int N, color[] cols) {
  int N0 = cols.length;
  float[] R0 = new float[N0], G0 = new float[N0], B0 = new float[N0], ind0 = new float[N0];
  for (int i=0; i<N0; i++) {
    R0[i] = red(cols[i]);
    G0[i] = green(cols[i]);
    B0[i] = blue(cols[i]);
    ind0[i] = (float)i / (N0-1);
  }
  color[] cmap = new color[N];
  for (int i=0; i<N; i++) {
    float ind = (float)i / (N-1);
    cmap[i] = color(interp1(ind0,R0,ind), interp1(ind0,G0,ind), interp1(ind0,B0,ind));
  }
  return cmap;
}

color[] warmrainbow() {return warmrainbow(20);}
color[] warmrainbow(int n) {return warmrainbow(n, 0.9, 0.3);}
color[] warmrainbow(int n, float bright, float contrast) {
  // nice colorscale from blue to yellow to red
  float Hblue = 0.55; // hues for end & middle colors
  float Hyellow = 0.16667;
  float Hred = 0;
  float dipWidth = 1.5; // width of the dip in saturation over green
  float[] H = new float[n];
  float[] S = new float[n];
  float[] B = new float[n];
  int N = n-1;
  int iy = floor(n/(float)2); // index of yellow, the middle color
  // hue
  for (int i=0; i<=iy; i++) H[i] = Hblue - (Hblue-Hyellow) *sq(i/(float)iy);
  for (int i=iy+1; i<n; i++) H[i] = lerp(H[iy],Hred,(i-iy)/((float)n-iy));
  //saturation
  // find greenest color
  int ig = 0;
  for (int i=1; i<n; i++) if (abs(H[i]-0.3333) < abs(H[ig]-0.3333)) ig = i;
  // gaussian dip in saturation
  for (int i=0; i<n; i++) S[i] = 1 - 0.5*exp(-dipWidth*sq(i/(float)ig-1));
  // brightness
  float b = 4*contrast/N;
  float a = -b/N;
  for (int i=0; i<iy; i++) B[i] = bright - lerp(contrast,0,i/((float)iy-1));
  for (int i=iy; i<n; i++) B[i] = a*sq(i) + b*i + bright - contrast;
  colorMode(HSB);
  color[] map = new color[n];
  for (int i=0; i<map.length; i++) map[i] = color(H[i]*255,S[i]*255,B[i]*255);
  colorMode(RGB);
  return map;
}


// math ---------------------------------------------------------------------------------------------------------------------


float Inf = 1./0.;

boolean isnan(float x) {
  return !((x<=0) || (x>0));
}

boolean isnan(float[] x) {
  boolean found = false;
  for (int i=0; i<x.length && !found; i++) found = isnan(x[i]);
  return found;
}

float sum(float[] a) {
  float s = 0;
  for (int i=0; i<a.length; i++) s += a[i];
  return s;
}

float mean(float[] a) {
  return sum(a) / a.length;
}

int[] round(float[] a) {
  int[] r = new int[a.length];
  for (int i=0; i<a.length; i++) r[i] = round(a[i]);
  return r;
}

float randn() {
  // returns a normally distributed random number with mean 0 and std dev 1
  return sqrt(-2*log(random(1)))*cos(TWO_PI*random(1));
}

float randn(float a, float b) {
  // for repeatable randn values
  float r = random(1);
  randomSeed(round((a+b)*1e6)); // temporarily take control of the random sequence
  float result = randn();
  randomSeed(round(1e6*r)); // lose control of it again
  return result;
}

float log10(float a) {
  return log(a)/log(10);
}

float[] series(float min, float stride, float max) {
  float[] S = new float[floor((max-min)/stride)+1];
  for (int i=0; i<S.length; i++) S[i] = min + i*stride;
  return S;
}

float[] rotatePolyX(float[] x, float[] y, float ang) { // rotates a polygon around the center of a box (0..1, 0..1)
  ang = -ang*PI/180;
  float[] xr = new float[x.length];
  for (int i=0; i<x.length; i++) {
    xr[i] = (x[i]-0.5)*cos(ang) - (y[i]-0.5)*sin(ang) + 0.5;
  }
  return xr;
}
  
float[] rotatePolyY(float[] x, float[] y, float ang) { // rotates a polygon around the center of a box (0..1, 0..1)
  ang = -ang*PI/180;
  float[] yr = new float[x.length];
  for (int i=0; i<x.length; i++) {
    yr[i] = (x[i]-0.5)*sin(ang) + (y[i]-0.5)*cos(ang) + 0.5;
  }
  return yr;
}


int[] findBeforeAfter(float[] X, float xi) {
  // returns the indices into X that bracket xi.
  // could replace this with a non-integer find().
  boolean failed = false;
  int Nx = X.length;
  int dir;
  // is X increasing or decreasing?
  if (X[Nx-1] > X[0]) {
    dir = 1;
    failed = ((xi<X[0]) || (xi>X[Nx-1]));
  } else {
    dir = -1;
    failed = ((xi>X[0]) || (xi<X[Nx-1]));
  }
  if (failed) {
    return new int[] {-1,-1};
  } else {
    // find the indices that bracket xi
    int nbefore = 0;
    int nafter = Nx-1;
    int nmid = (nbefore+nafter)/2;
    while (nafter-nbefore > 1) {
      if (dir==1) {
        if (xi>X[nmid]) {
          nbefore = nmid;
        } else {
          nafter = nmid;
        }
      } else {
        if (xi<X[nmid]) {
          nbefore = nmid;
        } else {
          nafter = nmid;
        }
      }
      nmid = (nbefore+nafter)/2;
    }
    return new int[] {nbefore, nafter};
  }
}


float interp1(float[] X, float[] F, float xi) {
  int[] iba = findBeforeAfter(X,xi);
  int i0 = iba[0];
  int i1 = iba[1];
  if (i0==-1) return 0./0.;
  float Fi = ( (X[i1]-xi) * F[i0] + (xi-X[i0]) * F[i1] ) / (X[i1] - X[i0]);
  return Fi;
}

float interp2(float[] X, float[] Y, float[][] F, float xi, float yi) {
  int[] iba = findBeforeAfter(X,xi);
  int i0 = iba[0];
  int i1 = iba[1];
  if (i0==-1) return 0./0.;
  int[] jba = findBeforeAfter(Y,yi);
  int j0 = jba[0];
  int j1 = jba[1];
  if (j0==-1) return 0./0.;
  // interp in x
  float F0 = ( (X[i1]-xi) * F[i0][j0] + (xi-X[i0]) * F[i1][j0] ) / (X[i1] - X[i0]);
  float F1 = ( (X[i1]-xi) * F[i0][j1] + (xi-X[i0]) * F[i1][j1] ) / (X[i1] - X[i0]);
  // interp in y
  float Fi = ( (Y[j1]-yi) * F0        + (yi-Y[j0]) * F1        ) / (Y[j1] - Y[j0]);
  return Fi;
}


// input --------------------------------------------------------------------------------------



float[] interpFromFile(String filename, float[] xi) {
  // reads a file with columns x,y and interpolates to xi
  String[] lines = loadStrings(filename);
  float[] x = new float[lines.length];
  float[] y = new float[lines.length];
  float xmin = Inf, xmax = -Inf;
  for (int i=0; i<lines.length; i++) {
    String[] a = splitTokens(lines[i]);
    x[i] = float(a[0]);
    y[i] = float(a[1]);
    xmin = min(xmin,x[i]);
    xmax = max(xmax,x[i]);
  }
  float[] yi = new float[xi.length];
  for (int i=0; i<xi.length; i++) {
    if (xi[i] < xmin) {
      yi[i] = interp1(x,y,xmin);
    } else if (xi[i] > xmax) {
      yi[i] = interp1(x,y,xmax);
    } else {
      yi[i] = interp1(x,y,xi[i]);
    }
  }
  return yi;
}



// text -------------------------------------------------------------------------------------

String nicenum(float a) {
  if (a==round(a)) {
    return "" + round(a);
  } else {
    return "" + a;
  }
}

void textVAlign(String S, float x, float y) {textVAlign(S,x,y,"baseline");}
void textVAlign(String S, float x, float y, String align) {
  if (align.equals("baseline")) {
    text(S,x,y);
  } else if (align.equals("middle")) {
    text(S,x,y+0.5*(textAscent()));
  } else if (align.equals("bottom")) {
    text(S,x,y-textDescent());
  } else if (align.equals("top")) {
    text(S,x,y+textAscent());
  }
}

int countLines(String S) {
  int c = 1;
  for (int i=0; i<S.length(); i++) {
    if (S.charAt(i) == '\n') c++;
  }
  return c;
}