varying vec2 destcoords; // destination location
uniform sampler2D src; // previous attractor

/*---- convert plane-to-texture -------*/
uniform float texscale,texscalei; 
float T(float x) { // plane to texture
  x*=texscale; 
  return x/sqrt(1.0+x*x)*0.5+0.5; 
} 
vec2 T(vec2 p) { return vec2(T(p.x),T(p.y)); }  
float P(float s) { // texture to plane
  float u=2.0*s-1.0; 
  return texscalei*u/sqrt(1.0-u*u); 
}
vec2 P(vec2 s) { return vec2(P(s.x),P(s.y)); }

/*---------------- w0 -----------------*/
uniform vec4 color0;
uniform vec2 w0x,w0y,w0o; // w0's parameters
uniform float w0j,w0v; 
float jacobian0(vec2 t) {return 1;}
vec4 f0(vec2 inv) { // runs at each inverse
  float area=1e-2+abs(w0j*jacobian0(inv));
  vec2 t=T(w0x*inv.x+w0y*inv.y+w0o);
  return (exp2(texture2D(src,t)*20)-1)/area;
}
vec4 nonlinear_inverse0(vec2 p) {
  p=p*w0v; // Draves' variation coefficient
  return f0(p);
}

/*---------------- w1 -----------------*/
uniform vec4 color1;
uniform vec2 w1x,w1y,w1o; // w1's parameters
uniform float w1j,w1v;
float jacobian1(vec2 t) {
  float r2=dot(t,t); 
  return (1-2*t.y*t.y/r2)/r2;
}
vec4 f1(vec2 inv) { // runs at each inverse
  float area=1e-2+abs(w1j*jacobian1(inv));
  vec2 t=T(w1x*inv.x+w1y*inv.y+w1o);
  return (exp2(texture2D(src,t)*20)-1)/area;
}
vec4 nonlinear_inverse1(vec2 p) {
  p=p*w1v; // Draves' variation coefficient
  float ix = 0.5/p.x; 
  float det=1 - 4*p.x*p.x*p.y*p.y; 
  if (det>=0) { 
    float sq = sqrt(det); 
    return f1(vec2(ix*(1 - sq),p.y)) 
          +f1(vec2(ix*(1 + sq),p.y)); 
  } else { return vec4(0); }
}

/* Combined inverse-sampling function */
vec4 sum_inverses(vec2 p) {
  vec4 sum=vec4(0);
  sum+=color0*nonlinear_inverse0(p);
  sum+=color1*nonlinear_inverse1(p);
  return log2(sum+1)*(1.0/20);
}

void main(void) { 
  gl_FragColor = sum_inverses(P(destcoords));
}