#define PI 3.1415926
/*
Analytical 2-order motion blur
Version: 0.2
Author: Fedor Petrov

Inspired  by https://www.shadertoy.com/view/MdSGDm

The blur achieved in https://www.shadertoy.com/view/MdSGDm is always
linear. This is OK for relatively slow motion, even fast motion 
looks nice but doesn't feel fully right. At high velocities 
the blured dots start to look like rigid sticks.

The blur in this shader follows the parabola, which is constructed 
taking into account local coordinate, local velocity  and local acceleration.
To find the distance from a point to a parabola and thus the area covered
with color one needs to solve the third order equation for time.
Here I was too lazy and made the assumtion, that t**3 term is negligible.
It works quite fine, however there are some artifacts (work in progress). 

I know that the code is not clean enough, there are some parts to
be refactored.

This shader, however, doesn't simulate a lot of things that become
complicated when the trajectory is not linear (e.g., color overlapping
for a particle that rapidly changes direction of motion etc). 
*/
mat2 rotate(float alpha) {
return mat2(sin(alpha),cos(alpha),cos(alpha),-sin(alpha));
}

float trace(vec2 p, vec2 r0, vec2 v0, vec2 a0, float R, float dt)
{
	float rat = 1.5* PI*R*R/(length(v0)*dt*2.0*R+PI*R*R);
    float a = 1.5*(a0.y*v0.y+a0.x*v0.x);
    float b = (v0.x*v0.x+v0.y*v0.y) + a0.y*(r0.y-p.y) + a0.x*(r0.x-p.x);
    float c = v0.x*(r0.x-p.x)+v0.y*(r0.y-p.y);
    float D = b*b-4.0*a*c;
    float t1 = dt*10.0;
    float t2 = dt*10.0;
    float T=dt*10.0;
    float L1 = 2.0*R;
    float L2 = 2.0*R;
    float LM1 = 2.0*R;
    float LM2 = 2.0*R;
    float coef = 0.0;
    vec2 V1 = v0-a0*dt;
    float RT1 = R/length(V1);
    vec2 V2 = v0+a0*dt;
    float RT2 = R/length(V2);
    vec2 R1 = vec2(2.0*R);
    vec2 R2 = vec2(2.0*R);
    vec2 RM1 = 0.5*a0*dt*dt - v0*dt+r0-p;
    vec2 RM2 = 0.5*a0*dt*dt + v0*dt+r0-p;
    LM2 = length(RM2);
    LM1 = length(RM1);
    if(D>=0.0)
    {
    t1 = (-b + sqrt(D))*0.5/a;
    t2 = (-b - sqrt(D))*0.5/a;
    R1 = 0.5*a0*t1*t1 + v0*t1+r0-p;
    R2 = 0.5*a0*t2*t2 + v0*t2+r0-p;
    L1 = length(R1);
    L2 = length(R2); 
	if((abs(t1)<dt && L1<R) 
      || (abs(t2)<dt && L2<R) 
      || (min(LM1,LM2)<R) 
      )
    {
    coef = 1.0;
    }   
    }
    else if(min(LM1,LM2)<R) {
        coef = 1.0;
        }
         
    if(coef!=0.0) {
        float mt = min(L1,L2);
        if(mt>min(LM1,LM2))
        {
        mt = min(LM1,LM2);
        }   
        if(LM1<R && LM2<R)
        {
        rat = 1.0;
        }
        else if(LM1<R) {
            float RMdir = dot(RM1,V1)/length(V1);
            float RMperp = sqrt(LM1*LM1-RMdir*RMdir);
            float fact = 0.5-0.5*RMdir/sqrt(R*R-RMperp*RMperp);
            coef *=sqrt(1.0-pow(clamp(RMperp,0.0,R)/R,2.0))*fact;
        }
        else if(LM2<R) {
            float RMdir = dot(RM2,V2)/length(V2);
            float RMperp = sqrt(LM2*LM2-RMdir*RMdir);
            float fact = 0.5+0.5*RMdir/sqrt(R*R-RMperp*RMperp);
            coef *=sqrt(1.0-pow(clamp(RMperp,0.0,R)/R,2.0))*fact;
        }
        else {
        coef *= sqrt(1.0-pow(mt/R,2.0));
        }
        }    
	return coef*rat;
}


float randomoff(float x)
{
return fract(271828.459045+3.1415926*sqrt(abs(x)));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 cent = vec2(0.5,0.5);
    float alpha = 2.0;
    float beta = 2.0;
    float ampxx = 0.3;
    float ampyy = 0.3;
    float ratio = iResolution.x/iResolution.y;
    vec2 uv = fragCoord.xy/iResolution.xy;
    fragColor = texture(iChannel0,uv);
    fragColor.rgb*=0.5;
    float maxx=30.0;
    for(float k=0.0;k<maxx;k++)
    {
    alpha = alpha + 0.111111;
    beta = beta + randomoff(beta);
    float tt = iTime*0.5;
    float phi = randomoff(k*2.71828182845*10.0);
    float ampx = ampxx*(0.9+0.2*randomoff(k*alpha));
    float ampy = ampyy*(0.9+0.2*randomoff(k*beta));
    vec2 r0 = vec2(ampx*sin(alpha*tt+phi),
                   ampy*cos(beta*tt+phi));
    vec2 v0 = vec2(alpha*ampx*cos(alpha*tt+phi),
                   -beta*ampy*sin(beta*tt+phi));
    vec2 a0 = vec2(-alpha*alpha*ampx*sin(alpha*tt+phi),
                   -beta*beta*ampy*cos(beta*tt+phi));
    a0 = a0*rotate(k*0.1);
    v0 = v0*rotate(k*0.1);
    r0 = r0*rotate(k*0.1)+vec2(ratio*0.5,0.5);
      
    float R = 0.02+k/maxx*0.03;
    float DT=0.02;
    // Normalized pixel coordinates (from 0 to 1)
    uv = fragCoord/iResolution.y;

    // Time varying pixel color
    vec3 col = 0.5 + 0.5*cos(iTime+uv.xyx+vec3(0,2,4)+phi*2.0);
	float mult = trace(uv,r0,v0,a0,R,DT);
    // Output to screen
    fragColor = fragColor*(1.0-mult) + vec4(col*1.5*vec3(k/maxx),1.0)*mult;
    
    }
}