// Created by Andrew Wild - akohdr/2021
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

// needs gamma correction

#define AA                       // uncomment to turn on anti aliasing
#define RADIAL_AA
#define FOUR_ROOK_AA

#define ITS 7                      // iteration count
#define SPEED 7e-2                 // global speed
#define ZOOM 100.                  // zoom level
//#define ORIGIN                   // highlight origin
//#define DEBUG                    // flash any isnan red

// performs two 2d offset rotations moving field origin in trefoil orbit
#define ROTATE(a) mat2(cos(a),-sin(a),sin(a),cos(a))
//#define ROTATE(a) 1.   // no rotation
//#define ROTATE(a) 2.   // see the bigger picture

#define iR iResolution

void calcImage( out vec4 k, in vec2 p )
{
    p*=ZOOM;

    float wt = iTime*SPEED;        // wall time

         wt += 2.5;                // starting time offset

          p *= .05;                // unzoom compensation
          p -= iR.xy/2.;           // center coords
          p += vec2(-1e3);         // offset
          p *= ROTATE(wt);         // outer rotation

     vec2 c  = p-iR.xy/2.;         // centered coords
          c /= iR.y;               // aspect corrected/scaled
          c *= ROTATE(-wt*3.);     // inner rotation
          c += vec2(1);            // movement

     float t = wt+.2;              // local time

     vec4  P = vec4(2,3,5,7),      // primes
          vx =  (P*c.x),           // x space component
          vy =  (P*c.y);           // y space component

           k = vec4(0,0,0,1);      // initial condition

    for(int i=ITS; i>0; i--)      // iterate mat4 pipeline
    {
      // main transform
      mat4 M = mat4(vx,vy,
                    t*P,
                    k*(c.xyyx-1.5));

      vec4 l = k;
      k *= M;
      k = sin(k/9.);

      #ifdef DEBUG
          if(any(isnan(k))) k = l;
      #endif
    }

//    k = vec4(.5/length(k));    // grayscale
//    k = normalize(k);

    // show origin
    #ifdef ORIGIN
        k.g += smoothstep(.1,.0,length(c));
    #endif

    // flash any isnan values red
    #ifdef DEBUG
       if(any(isnan(k))) k = vec4(mod(iTime,1.),0,0,0);
    #endif
}

void mainImage( out vec4 k, in vec2 p )
{
// shader is pretty much torture test for aliasing
#ifdef AA
    vec4 s = k = vec4(0);

#ifdef RADIAL_AA
    for (float a=6.18; a>0.; a-=1.2){
        calcImage(s,p+(a/9.)*vec2(cos(a), sin(a)));
        k = mix(k,s,.5);
    }
#endif

#ifdef FOUR_ROOK_AA
    // apply four rook anti-aliasing
    //   improves but doesn't eliminate issue in high contrast/freq. areas
    float d1 = .125, d2 = .375;
    calcImage(s,p+vec2( d1, d2));  k = mix(k,s,.5);
    calcImage(s,p+vec2( d2,-d1));  k = mix(k,s,.5);
    calcImage(s,p+vec2(-d2, d1));  k = mix(k,s,.5);
    calcImage(s,p+vec2(-d1,-d2));  k = mix(k,s,.5);
#endif

#else
    calcImage(k,p);
#endif
}