#define too_small 1e-8

//common math objects everyone knows about and can use freely

vec3 angle_axis_rotation( vec3 v, vec3 k, float a ) {

    //rotates vector:v around axis: k (normalized) by angle:a
    //https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
    //Euler already knew it though..
    
    //Please use it, Abuse it!
    float cc=cos(a), ss=sin(a);
    
    return v*cc + cross(k,v)*ss + k*dot(k,v)*(1.-cc);
}


mat3 inv3d( mat3 a ) {

    //good ol' fashioned inverse of 3x3 matrix
    //Please use it, Abuse it!

    //determinant
    float det = 
        a[0][0] * a[1][1] * a[2][2] 
      + a[1][0] * a[2][1] * a[0][2]
      + a[2][0] * a[0][1] * a[1][2]
      - a[2][0] * a[1][1] * a[0][2]
      - a[1][0] * a[0][1] * a[2][2]
      - a[0][0] * a[2][1] * a[1][2];
 
    if ( abs(det) < too_small ) 
    { det = sign(det) * 1./too_small; }
    else { det = 1./det; }

    return det * 
    mat3( a[1][1]*a[2][2]-a[2][1]*a[1][2], a[2][0]*a[1][2]-a[1][0]*a[2][2], a[1][0]*a[2][1]-a[2][0]*a[1][1],
          a[2][1]*a[0][2]-a[0][1]*a[2][2], a[0][0]*a[2][2]-a[2][0]*a[0][2], a[2][0]*a[0][1]-a[0][0]*a[2][1],
          a[0][1]*a[1][2]-a[1][1]*a[0][2], a[1][0]*a[0][2]-a[0][0]*a[1][2], a[0][0]*a[1][1]-a[1][0]*a[0][1] );

}

#define newton_iter 7.
float timefac;

vec3 newton_color( vec3 p ) {

    p = mod(p-2.,4.)-2.;
    vec3 q, color=vec3(0.); mat3 dq;
    for (float i=0.;i<newton_iter;i++) {
    
        //you can plug and play with your own 
        //vector functions and derivatives here
        float pxx=p.x*p.x,pyy=p.y*p.y,pzz=p.z*p.z;
        //vector function
        q = vec3(p.y*pyy, p.z*pzz, p.x*pxx) - vec3(-3.1);
        //matrix of partial derivatives
        dq = mat3( 0.,0.,3.*pxx,
                   3.*pyy,0.,0.,
                   0.,3.*pzz,0.);
                   
        p -= (1.-timefac)*inv3d(dq) * q;  //the "Newton Step"
        
        if ( p.x*p.y > 0. ) color.x++;
        if ( p.y*p.z > 0. ) color.y++;
        if ( p.z*p.x > 0. ) color.z++;
    }
    
    return color/newton_iter;

}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{

    float mt = mod(iTime,30.);
    vec3 ro = vec3(0.,0.,0.);

    if (mt > 12. ) {
        float t1 = 3.;
        float tt = mod( (mt-12.)/3.,2.*t1);
    
        timefac = tt;
        if ( tt>t1 ) timefac = t1 - mod(tt,t1);
    }
    else if (mt>2.) {
    
        float xx = sin( (mt-2.)/2.*1.2);
        ro = vec3(-xx,-xx,-xx );
    
    }
    
    
    vec2 uv = (2.*fragCoord - iResolution.xy)/iResolution.y/3.*(1.+timefac);
    
    vec3 rd = normalize( vec3(uv,-1.4) );
    
    vec3 p = ro;
    
    //get oriented to the good place using angle axis rotation
    //overkill here but I wanted to code it
    rd = angle_axis_rotation(rd, vec3(1.,0.,0.), -.65);
    rd = angle_axis_rotation(rd, vec3(0.,1.,0.), .8);
    
    vec3 color = vec3(0.);
    float totdist = 0.;
    float stepsize = .015*(1.+timefac/5.);
    
    //doing raymarch not measuring any distances but accumulating
    //color along the path, more of a volumetric density kind of thing
    for (float i=0.; i<80.; i++) {
    
        totdist += stepsize;
        float tt = totdist*totdist; tt*=tt;
        color += newton_color(p) / (1.+tt);
        p += stepsize * rd;
    }
    
    color = cos(color/32.);
    color *= color;
    
    fragColor = vec4(color,1.);

}