#define t iTime*.5
#define PI 3.14159265359
const float radius=0.15;
const float M=1.;
const float c=1.;
const float freq=2.;


const int iterations=20; //try 15 iterations for cool clitch :)
vec2 p10=radius*vec2(1,0);
vec2 p20=radius*vec2(-1,0);


mat2 rot(float a){
    return mat2(cos(a), -sin(a),sin(a),cos(a));
}

vec2 getPos(vec2 p0, float time){
    return rot(freq*time)*p0;
}


vec2 retardedPos(vec2 uv, vec2 p0, float time){
    float upper =time;
    float lower =0.;
    float tr = time*.5;
    int n=0;
    vec2 retardedP=getPos(p0, tr);
    while(n<iterations){
        retardedP=getPos(p0, tr);
        float dist= length(uv-retardedP);
        if(dist/c>(time-tr)){
            upper = tr;
            tr=(upper+lower)*.5;
        }else{
            lower= tr;
            tr=(upper +lower)*.5;
        }if(tr<.01){
           tr=0.;
           break;
        }
        retardedP=getPos(p0, tr);
        n++;
    }

    return retardedP;
}


vec2 force1(vec2 uv, vec2 p){
//this is the good old newtonian force:
    vec2 separation=uv-p;
    float dist=length(separation);
    float scl=0.01;
    return -scl*M*pow(dist,-2.)*normalize(separation);
}

//these are relativistic extra components:
vec2 force2(vec2 uv, vec2 p){
    vec2 separation=uv-p;
    float dist=length(separation);
    vec2 a = -freq*freq*p;
    float scl=.05;

    return scl*4.*M/(c*c*dist)*a;;
}

vec2 force3( vec2 uv, vec2 p){
    vec2 separation=uv-p;
    float dist=length(separation);
    vec2 v = freq*rot(PI*.5)*p;
    vec2 a = -freq*freq*p;
    float scl=.1;

    return scl*4.*M*dot(a,separation)/(c*c*c*dist*dist)*v; // +1./(c*c*dist)*a);
}


#define FORCE force3(p.xz, retardedP1)+force3(p.xz, retardedP2)
//#define FORCE force2(p.xz, retardedP1)+force2(p.xz, retardedP2)
//#define FORCE force1(p.xz, retardedP1)+force1(p.xz, retardedP2)

#define MAX_ITER 200.
#define MAX_DIST 10.
#define SURF .005

#define BODY 1
#define SURFACE -1

vec3 col=vec3(0);

vec2 force;
vec2 retardedP1,retardedP2;


vec3 getRayDir(vec2 uv, vec3 ro,vec3 lookAt, float zoom){

    vec3 f= normalize(lookAt-ro);
    vec3 r= normalize(cross(vec3(0,1,0),f));
    vec3 u= cross(f,r);
    vec3 i= ro+f*zoom+uv.x*r+uv.y*u;

    return normalize( i-ro);
}

float sdSphere(vec3 p, vec3 center, float rad){
return length(p-center)-rad;
}


float sdPlane(vec3 p){
    vec3 s= vec3(9,2.,9);
    p.y+=s.y*.5;
    s*=.5;

    p= abs(p)-s;
    return length(max(p,0.))+ min(max(p.x,max(p.y,p.z)),0.);

}

int getMaterial(vec3 p){
    if(sdSphere(p, vec3(retardedP1.x,0.02,retardedP1.y), .06)<7.*SURF){
        return BODY; //1 ;
    }else if(sdSphere(p, vec3(retardedP2.x,0.02,retardedP2.y), .06)<7.*SURF){
        return BODY; //-1;
    }
    else {
        return SURFACE;
    }

}

float getDist(vec3 p){

   float dist =sdPlane(p);

   vec2 p1=getPos(p10, t);
   vec2 p2=getPos(p20, t);

   retardedP1=retardedPos(p.xz, p1, t);
   retardedP2=retardedPos(p.xz, p2, t);


    if(dist<MAX_DIST){
        //bumpmap:
        force=FORCE;

        float displacement= length(force);  //*sign(p.y);
        dist+=displacement;

    }else
        return MAX_DIST;

    //the mass bodies:
    float material = float(getMaterial(p));
    float sd1= -material* sdSphere(p, vec3(retardedP1.x,0.02,retardedP1.y), .06);
    float sd2= -material* sdSphere(p, vec3(retardedP2.x,0.02,retardedP2.y), .06);

    return min(sd2,min(sd1,dist));
}


vec3 getNormal(vec3 p){
  vec2 e= vec2(.01,0);
   float d=getDist(p);
   vec3 n = d-vec3(getDist(p- e.xyy),getDist(p- e.yxy),getDist(p- e.yyx));

   return normalize(n);
}

float RayMarch(vec3 ro, vec3 rd, float inside){
    float dO=0.;
    float i=0.;
   while(i<MAX_ITER){
      vec3 p= ro+dO*rd;

      float dS=inside*getDist(p);

        //conservative stepsize closer to bodies:
      float d1 =length(p-vec3(retardedP1.x,.02,retardedP1.y));
      float d2 =length(p-vec3(retardedP2.x,.02,retardedP2.y));

      if(d1<.9)
        dS*=(.1+d1);
      else if(d2<.9)
        dS*=(.1+d2);

      dO+=dS;

      if(dO>MAX_DIST){
          float halo= 3.*i/100.;
          col+=halo*halo*vec3(0.4,0.2,1);
          break;
          }
      else if(dS<SURF){
          float iter= i+2.*dS/SURF; //smoothiteration
          float halo= 11.*iter/MAX_ITER;
          col+=halo*vec3(0.4,0.2,1);
          break;
      }
      i++;
    }



      return dO;
}




void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord-iResolution.xy*.5)/iResolution.y;
    vec2 m = (iMouse.xy-.5)/iResolution.xy;


    //camera:
    float zoom= 1.;//smoothstep(-1.,3.,iTime);

    vec3 ro= vec3(-1.,.9,-0.);
    ro.xz*=rot(PI/2.);



    if(sign(iMouse.z)>0.){
        ro.yz*=rot(-(m.y-.5)*PI);
        ro.xz*=rot(-(m.x-.5)*PI);
    }else{
        ro.yz*=rot(sin((t*.5-.5)*PI)*.5-PI*.2);
    }


    vec3 lookAt=vec3(0,.1,0);

    vec3 rd= getRayDir(uv, ro, lookAt,zoom);

    float d= RayMarch(ro,rd,sign(ro.y));

    vec3 p=ro;

     if(d<MAX_DIST){//if we hit the object:

          p= p+ d*rd;

          //next march if we begin below surface:
          if(sign(ro.y)<=0.){
             vec3 n=getNormal(p);

             //float material = float(getMaterial(p));
             d= RayMarch(p+5.*SURF*n, rd,1.);
             p= p+ d*rd;
             /*
             if(getMaterial(p)==BODY){
              col=vec3(0);
              }
             */

         }
         /*
          else if(getMaterial(p)==BODY){
              col=vec3(1);
          }
          */

    }

    fragColor = vec4(col,1.);
}