#version 430

// layout (location=0) in vec3 vertexPosition;

out vec2 uv;

out float no;

out vec3 norm;
out vec3 tan;
out vec3 bit;

uniform mat4 projectionMatrix;
uniform mat4 modelViewMatrix;
uniform mat4 viewMatrix;
uniform mat4 projectionInvMatrix;
uniform float g_windowWidth;
uniform float g_windowHeight;

uniform float g_time;

vec2 rot2(vec2 p, float a) {
  vec2 r = p;
  r.x = cos(a)*p.x - sin(a)*p.y;
  r.y = sin(a)*p.x + cos(a)*p.y;
  return r;
}

vec4 rotateXZ(vec4 p, float a) {
    vec4 r = p;
    r.x = cos(a)*p.x - sin(a)*p.z;
    r.z = sin(a)*p.x + cos(a)*p.z;
    return r;
}
vec4 rotateYZ(vec4 p, float a) {
    vec4 r = p;
    r.y = cos(a)*p.y - sin(a)*p.z;
    r.z = sin(a)*p.y + cos(a)*p.z;
    return r;
}

// from pouet raymarching thread by las of mercury
float perlin(vec3 p) {
    vec3 i = floor(p);
    vec4 a = dot(i, vec3(1., 57., 21.)) + vec4(0., 57., 21., 78.);
    vec3 f = cos((p-i)*3.14159265)*(-.5)+.5;
    a = mix(sin(cos(a)*a),sin(cos(1.+a)*(1.+a)), f.x);
    a.xy = mix(a.xz, a.yw, f.y);
    return mix(a.x, a.y, f.z);
}
float turb(vec3 c) {
        float r=0.0;
        float s=0.5;
        c*=1.0;
        for (int i=0;i<4;i++) {
                //if (i<2) continue;
                r+=s*(perlin(c)*1.0+0.0);
                c*=3.0;
                s*=0.65;
        }
        return r;
}

vec3 posd(float pw) {
    return vec3(sin(pw)*17.0, cos(pw*1.7)*13.0, sin(pw*0.78+cos(pw*0.4))*11.0)*32.0;
}

float getSplineValue(float t, float p0, float pt0, float p1, float pt1) {
    float t2 = t*t;
    float t3 = t*t*t;
    float p = (2.0*t3-3.0*t2+1.0)*p0+(t3-2.0*t2+t)*pt0+(-2.0*t3+3.0*t2)*p1+(t3-t2)*pt1;
    return p;
}

uniform float splineEnabled = 0.0;


float perlin4(vec4 p) {
//    vec4 pf1 = floor(p);
//    vec4 pf2 = floor(p+vec4(1.0));

    float p1 = perlin(p.xyz+posd(floor(p.w)+1.0));
    float p2 = perlin(p.xyz+posd(floor(p.w)+2.0));

    if (splineEnabled > 0.5) {
        float p0 = perlin(p.xyz+posd(floor(p.w)));
        float p3 = perlin(p.xyz+posd(floor(p.w)+3.0));
        float ki = fract(p.w);
        return getSplineValue(ki, p1, p2-p0, p2, p3-p1);
    } else {
        float ki = smoothstep(0.0, 1.0, fract(p.w));
        return p1*(1.0-ki)+p2*(ki);
    }
}

uniform float noiseOct = 4.0;
uniform float noiseOctMul = 3.0;
uniform float noiseOctAmp = 0.75;

float turb4(vec4 c) {
        float r=0.0;
        float s=0.5;
        for (int i=0;i<int(noiseOct);i++) {
                r+=s*(perlin4(c));
                c*=noiseOctMul;
                s*=noiseOctAmp;
        }
        return r;
}



float zn = 0.0;
float zf = 1.0;

vec4 CalcEyeFromWindow(in vec3 windowSpace) {
    vec3 ndcPos;
    vec4 viewport = vec4(0.0, 0.0, g_windowWidth, g_windowHeight);
    ndcPos.xy = ((2.0 * windowSpace.xy) - (2.0 * viewport.xy)) / (viewport.zw) - 1;
    ndcPos.z = (2.0 * windowSpace.z - zn - zf) / (zf - zn);
    vec4 clipPos;
    clipPos.w = projectionMatrix[3][2]/(ndcPos.z-(projectionMatrix[2][2]/projectionMatrix[2][3]));
    clipPos.xyz = ndcPos * clipPos.w;
    return projectionInvMatrix * clipPos;
}

void calcIdo(inout float id, inout float ido, inout float idoSub, inout float objMirror, float numVert, float genObs) {
    //float mirror = objParams.x;
    float mirror = 0;
    float numVertMirror = numVert*(floor(mirror));
    float idom = floor(id/numVertMirror); // object id mirrored
    float idm = id-idom*numVertMirror; // vertex id mirrored

    ido = floor(id/numVert);
    id = id-ido*numVert;

    objMirror = floor(idm/numVert);

    float idoC = floor(idom/genObs); // object id combining all sub objects a single object consists of
    idoSub = idom-idoC*genObs;
    ido = idoC;
}


float getPlaneNV(float rows, float cols) {
    return (rows*cols*6.0);
}
// id: vertex id number
vec4 getPlane(float id, float rows, float cols, inout vec4 n, inout vec2 tc, inout vec4 ta, inout vec4 qc) {
    float id3 = fract(id/3.0)*3.0; // index of the vertex in a triangle 0,1,2
    float idv = floor((id/3.0-fract(id/3.0))+0.5); // triangle index 0,1,2,3,...
    float idvq = floor((id/6.0-fract(id/6.0))+0.5); // quad index 0,1,2,3,...
    float idv2 = fract(idv/2.0)*2.0; // every second triangle 0 and 1

    vec4 vp = vec4(0.0, 0.0, 0.0, 0.0);

    float triScale = 1.0;
    float gap = 1.0;

    if (idv2 < 0.5) {
        if (id3 < 0.50) {
            vp = vec4(0.0, 0.0, 0.0, 1.0);
        } else if (id3 < 1.5) {
            vp = vec4(1.0, 0.0, 0.0, 1.0);
        } else if (id3 < 2.5) {
            vp = vec4(0.0, 1.0, 0.0, 1.0);
        }
    } else {
        if (id3 < 0.50) {
            vp = vec4(1.0, 0.0, 0.0, 1.0);
        } else if (id3 < 1.5) {
            vp = vec4(1.0, 1.0, 0.0, 1.0);
        } else if (id3 < 2.5) {
            vp = vec4(0.0, 1.0, 0.0, 1.0);
        }
    }
    qc = vec4(0.5, 0.5, 0.0, 1.0);
    vp *= triScale;
    qc *= triScale;

    n = vec4(0.0, 0.0, -1.0, 0.0);
    ta = vec4(1.0, 0.0, 0.0, 0.0);


    float col = fract(idvq/cols)*cols;
    float row = floor(idvq/cols);

    vec2 vpo = vec2(col*triScale*gap-cols*0.5*triScale, row*triScale*gap-rows*0.5*triScale);

    vp.xy += vpo.xy;
    vp.xy /= vec2(cols, rows);
//    vp.xy /= vec2(cols, rows);

    qc.xy += vpo.xy;
    qc.xy /= vec2(cols, rows);

//    vp.x += col*triScale*gap-cols*0.5*triScale;
//    vp.y += row*triScale*gap-rows*0.5*triScale;

    //vp.x /= cols;
    //vp.y /= rows;


    tc = vec2(vp.x+0.5, vp.y+0.5);

    if (row > rows) vp *= 0.0;
    if (row > rows) qc *= 0.0;
    //  vp.xyz += vec3(sin(2.3*idv), cos(3.4*idv), sin(cos(idv*2.1)*1.5+idv*1.7));

    vp.w = 1.0;
    qc.w = 1.0;

    return vp;
}



float getCubeNV(float rows, float cols) {
    return (rows*cols*6.0*6.0);
}
vec4 getCube(float id, float rows, float cols, inout vec4 n, inout vec2 tc, inout vec4 ta, inout vec4 qc) {
    float idc = floor(id/(rows*cols*6.0)-fract(id/(rows*cols*6.0))+0.5);
    float ids = id-idc*rows*cols*6.0; // vertex in cube side
    vec4 vp = getPlane(ids, rows, cols, n, tc, ta, qc);
    vp.z = -0.5;
    qc.z = -0.5;

    if (idc < 0.5) {
    } else if (idc < 1.5) {
        vp = rotateXZ(vp, 0.5*3.141592);
        n = rotateXZ(n, 0.5*3.141592);
        ta = rotateXZ(ta, 0.5*3.141592);
        qc = rotateXZ(qc, 0.5*3.141592);
    } else if (idc < 2.5) {
        vp = rotateXZ(vp, 3.141592);
        n = rotateXZ(n, 3.141592);
        ta = rotateXZ(ta, 3.141592);
        qc = rotateXZ(qc, 3.141592);
    } else if (idc < 3.5) {
        vp = rotateXZ(vp, 1.5*3.141592);
        n = rotateXZ(n, 1.5*3.141592);
        ta = rotateXZ(ta, 1.5*3.141592);
        qc = rotateXZ(qc, 1.5*3.141592);
    } else if (idc < 4.5) {
        vp = rotateYZ(vp, 0.5*3.141592);
        n = rotateYZ(n, 0.5*3.141592);
        ta = rotateYZ(ta, 0.5*3.141592);
        qc = rotateYZ(qc, 0.5*3.141592);
    } else if (idc < 5.5) {
        vp = rotateYZ(vp, -0.5*3.141592);
        n = rotateYZ(n, -0.5*3.141592);
        ta = rotateYZ(ta, -0.5*3.141592);
        qc = rotateYZ(qc, -0.5*3.141592);
    } else {
        vp = vec4(0.0, 0.0, 0.0, 1.0);
    }

    vp.w = 1.0;

    return vp;
}



float getBallNV(float rows, float cols) {
    return (rows*cols*6.0);
}

vec4 getSpherical(vec3 v, float rad) {
   v.x += 0.5;
   v.y *= 2.0;
   v.z = 0.0;

   float yr = sqrt(rad*rad-v.y*v.y);

   vec4 vpr = vec4(yr*cos(v.x*2.0*3.141592), v.y, yr*sin(v.x*2.0*3.141592), 1.0);
   vpr.xyz *= 0.5;

   return vpr;
}

vec4 getBall(float id, float rows, float cols, inout vec4 n, inout vec2 tc, inout vec4 ta, inout vec4 qc) {
    vec4 vp = getPlane(id, rows, cols, n, tc, ta, qc);
    vp.x += 0.5; // 0 .. 1 range
    vp.y *= 2.0; // -1 .. 1 range
    vp.z = 0.0;

    float rad = 1.0;

    // y^2+yr^2 = rad^2
    // yr^2 = rad^2-y^2;
    float yr = sqrt(rad*rad-vp.y*vp.y);

    vec4 vpr = vec4(yr*cos(vp.x*2.0*3.141592), vp.y, yr*sin(vp.x*2.0*3.141592), 1.0);
    n.xyz = vpr.xyz;

    ta = rotateXZ(n, tc.x*2.0*3.141592);
    ta.y = 0.0;

    vpr.xyz *= 0.5;


    qc = getSpherical(qc.xyz, rad);

    return vpr;
}




uniform float div = 4.0;
uniform float noiseAmp = 1.0;
uniform float noiseFreq = 1.0;


uniform float noisePos = 0.0;
uniform float noiseOfsX = 0.0;
uniform float noiseOfsY = 0.0;
uniform float noiseOfsZ = 0.0;

uniform float minLim = 0.02;
uniform float fadeArea = 0.20;
uniform float fadeAmp = 160.0;

uniform float explo = 0.0;

uniform float noiseExp = 1.0;
uniform float noiseSub = 0.0;

uniform float subFreq = 4.0;
uniform float subSinFreq = 1.0;
uniform float subAmp = 0.250;
uniform float subEps = 0.050;
uniform float subPos = 0.0;

uniform float subBaseNorm = 0.0;

uniform float uvScale = 2.0;

uniform float scaleX = 1.0;
uniform float scaleY = 1.0;
uniform float scaleZ = 1.0;

float t4(vec3 p, float fr, float pos) {
    return turb4((vec4(p.xyz, 0.0)+vec4(noiseOfsX, noiseOfsY, noiseOfsZ, noisePos+pos))*fr);
}

void tn4(vec3 p, float eps, inout float vbase, inout vec3 von, float fr, float pos) {
    vbase = t4(p.xyz, fr, pos);
    float vox = t4(p.xyz+vec3(eps, 0.0, 0.0), fr, pos);
    float voy = t4(p.xyz+vec3(0.0, eps, 0.0), fr, pos);
    float voz = t4(p.xyz+vec3(0.0, 0.0, eps), fr, pos);
    von = vec3(vox, voy, voz)-vec3(vbase);
    von = normalize(von);
}


void tn4ntb(vec3 p, float eps, inout float vbase, inout vec3 von, float fr, float pos, vec3 n, inout vec3 ta, inout vec3 tb, float noiAmp, float noiSub, float noiExp) {
    vbase = t4(p.xyz, fr, pos);
    float vox = t4(p.xyz+ta*eps, fr, pos);
    float voy = t4(p.xyz+tb*eps, fr, pos);
    float voz = t4(p.xyz+n*eps, fr, pos);
    vec3 vob = vec3(vox, voy, voz)-vec3(vbase);

 //   float noa = noiAmp-noiseSub*clamp(-vbase, 0.0, 1024.0);
 //   noa = clamp(noa, 0.0, 1024.0);

    float no = (vbase-noiSub)*noiAmp;
    no = sign(no)*pow(abs(no), noiExp);
    no = clamp(no, minLim, 1024.0);
    vbase = no;

    float sur = noiAmp*(1.0-pow(clamp((-vbase+noiSub)*noiAmp, 0.0, 1.0), noiExp));

    // vob.xy *= ;
    // vob.xy *= noiseAmp*(0.0+clamp(no, 0.0, 1.0));
    von = -(vob.x*ta*sur+vob.y*tb*sur-n*eps);
    von = normalize(von);


    vec3 nta = -(-ta*eps+vob.y*tb*sur+n*vob.z*sur);
//    vec3 nta = ta;

    nta.xyz = nta.xyz-dot(-von.xyz, nta.xyz)*von.xyz;

    ta = normalize(nta);

    tb = cross(von, ta);

  //  no *= noiAmp;
  //  vbase = no;

}



void main() {
  vec4 vp = vec4(0.0);

  float id = float(gl_VertexID);

  float id3 = fract(id/3.0)*3.0; // index of the vertex in a triangle 0,1,2
  float idv = floor((id/3.0-fract(id/3.0))+0.5); // triangle index 0,1,2,3,...
  float idv2 = fract(idv/2.0)*2.0; // every second triangle 0 and 1
  float idq = floor((id/6.0-fract(id/6.0))+0.5); // quad index 0,1,2,3,...


  vec2 tc;
  vec4 ta;
  vec4 qc;
  vec4 vpN = vec4(0.0);

  float numVert = 0.0;
  float ido = 0.0;
  float idoSub = 0.0;
  float objMirror = 0.0;
  float genObs = 1.0;


  float rows = div;
  float cols = div;



//  numVert = getCubeNV(rows, cols);
//  calcIdo(id, ido, idoSub, objMirror, numVert, genObs);
//  vp = getCube(id, rows, cols, vpN, tc, ta, qc);

  numVert = getBallNV(rows, cols);
  calcIdo(id, ido, idoSub, objMirror, numVert, genObs);
  vp = getBall(id, rows, cols, vpN, tc, ta, qc);


  vp.xyz /= length(vp.xyz);
  qc.xyz /= length(qc.xyz);

  vpN.xyz = vp.xyz;

  vp.xyz *= vec3(scaleX, scaleY, scaleZ);
  vpN.xyz /= vec3(scaleX, scaleY, scaleZ);

  vpN.xyz = normalize(vpN.xyz);

//  ta.xyz = normalize(ta.xyz);

  ta.xyz = ta.xyz-dot(-vpN.xyz, ta.xyz)*vpN.xyz;

  ta.xyz = normalize(ta.xyz);

  vec3 tb = cross(vpN.xyz, ta.xyz);

  float vbase;
  vec3 von;

//  no = turb4((vec4(vp.xyz, 0.0)+vec4(noiseOfsX, noiseOfsY, noiseOfsZ, noisePos))*noiseFreq)*noiseAmp;
  //tn4(vp.xyz, subEps, vbase, von, noiseFreq, 0.0);
  tn4ntb(vp.xyz, subEps, vbase, von, noiseFreq, 0.0, vpN.xyz, ta.xyz, tb, noiseAmp, noiseSub, noiseExp);
//  no = vbase*noiseAmp;
//  no = sign(no)*pow(abs(no), noiseExp)-noiseSub;

  no = vbase;

  float noo = no;
  no = clamp(no, minLim, 1024.0);

  norm = von;
//  norm = vpN.xyz+von*no;
//  norm = vpN.xyz;

  vp.xyz += vpN.xyz*no;

//  float vbase = t4(vp.xyz);
//  float vox = t4(vp.xyz+vec3(eps, 0.0, 0.0));
//  float voy = t4(vp.xyz+vec3(0.0, eps, 0.0));
//  float voz = t4(vp.xyz+vec3(0.0, 0.0, eps));
//  vec3 von = vec3(vox, voy, voz)-vec3(vbase);

//  norm = von;

  norm = normalize(norm);

  if (abs(subAmp) > 0.001) {
  //    tn4(vp.xyz, subEps, vbase, von, noiseFreq);
  //    vbase = t4(vp.xyz, noiseFreq*subFreq);

//      tn4(vp.xyz, subEps, vbase, von, subFreq, subPos);
//      norm = norm*(subBaseNorm+clamp(1.0-subAmp*no, 0.0, 1.0))+von*subAmp*no;
//      norm = normalize(norm);
//      vp.xyz += subAmp*norm*clamp(vbase, 0.0, 1.0); // *sign(vbase)*pow(abs(vbase), 1.0);

      tn4(vp.xyz, subEps, vbase, von, subFreq, subPos);
      norm = norm*(subBaseNorm+clamp(1.0-subAmp*no, 0.0, 1.0))+von*subAmp*no;
      norm = normalize(norm);
      //vp.xyz += subAmp*norm*clamp(sin(no*0.10), 0.0, 1.0); // *sign(vbase)*pow(abs(vbase), 1.0);
      vp.xyz += subAmp*norm*clamp(sin(no*subSinFreq), 0.0, 1.0); // *sign(vbase)*pow(abs(vbase), 1.0);
      //vp.xyz += subAmp*norm*clamp(no, 0.0, 1.0); // *sign(vbase)*pow(abs(vbase), 1.0);


//     tn4ntb(vp.xyz, subEps, vbase, von, subFreq, subPos, norm.xyz, ta.xyz, tb, subAmp, noiseSub*0.0, noiseExp*0.0);
//     norm = von;
//     vp.xyz += subAmp*norm*clamp(no, 0.0, 1.0); // *sign(vbase)*pow(abs(vbase), 1.0);


      tn4(vp.xyz, subEps, vbase, von, subFreq, subPos);
      norm += von*subAmp*clamp(no, 0.0, 1.0);

  }


  norm = (modelViewMatrix*vec4(norm, 0.0)).xyz;


  norm = normalize(norm);

  tan = ta.xyz;
  bit = tb;

  if (abs(explo) > 0.001) {
      float noq = turb4((vec4(qc.xyz, 0.0)+vec4(noiseOfsX, noiseOfsY, noiseOfsZ, noisePos))*noiseFreq)*noiseAmp;
      noq = sign(noq)*pow(abs(noq), noiseExp);
      float noqo = noq;
      noq = clamp(noq, minLim, 10.0);
      vp.xyz += qc.xyz*explo*noq;
  }

  no = noo;

  if (no < (minLim*(1.0+fadeArea))) {
      no = minLim+(minLim*(1.0+fadeArea)-no)*fadeAmp;
  }

  no = clamp(no, -1000.0, 1000.0);


  vp = modelViewMatrix*vp;
  vp = projectionMatrix*vp;



  uv = tc;

  gl_Position = vp;
}
