#version 410 core
#define WELL_GL_PATCH_SIZE 16

layout(quads, equal_spacing, cw) in;

uniform vec4 cRadiusScale;
uniform vec3 cSceneScale;
uniform mat4 modelViewProjectionMatrix;
uniform mat4 modelViewMatrixInverse;
uniform mat4 modelViewMatrix;

in vec4 outpos_ES_in[];
in vec2 outAziIncl_ES_in[];
in float outProp_ES_in[];

out float rLogValueTransferFunction;
out float rTexU;

out vec3 vDelta;
out vec3 vCenter;
out float vLength;
out float rDiscard;


void main()
{
  float
    rIndex = gl_TessCoord.y * float(WELL_GL_PATCH_SIZE - 1);

  int
    iIndex = int(floor(rIndex + 0.01f));

  vec3
    cCurrVertex = outpos_ES_in[iIndex].xyz;

  vec2
    cAziIncl = outAziIncl_ES_in[iIndex];

  rLogValueTransferFunction = outProp_ES_in[iIndex];

  float
    rRadiusLog = outpos_ES_in[iIndex].w;

  float
    rInnerOuterRadiusFac = 1.0;
  
  rDiscard = 1.0;

  // are we at novalue?
  if (rRadiusLog > 2.0)
  {
    rInnerOuterRadiusFac = 0.0;
    rDiscard = 0.0;

    int
      iReadIndex = iIndex;

    // get previous
    if (iIndex > 0)
    {
      rRadiusLog = outpos_ES_in[iIndex - 1].w;
      iReadIndex = iIndex - 1;
    }

    // If not valid, try to get next.
    if (rRadiusLog > 2.0)
    {
      if (iIndex < (WELL_GL_PATCH_SIZE - 1))
      {
        rRadiusLog = outpos_ES_in[iIndex + 1].w;
        iReadIndex = iIndex + 1;
      }

      if (rRadiusLog > 2.0)
      {
        rRadiusLog = 0.0;
      }
    }

    cCurrVertex = outpos_ES_in[iReadIndex].xyz;
    cAziIncl = outAziIncl_ES_in[iReadIndex];
    rLogValueTransferFunction = outProp_ES_in[iReadIndex];
  }

  vec2
    cAzimuthVector = vec2(sin(cAziIncl.x), cos(cAziIncl.x));

  vec3
    cView = modelViewMatrixInverse[3].xyz - cCurrVertex.xyz,
    cViewOrthographic = modelViewMatrixInverse[2].xyz;

  //cView;// = mix(cView, cViewOrthographic, min(0.999999, projectionMatrix[3].w));

  mat3 
    cMat;

  cMat[2] = normalize(vec3(cAzimuthVector.xy * sin(cAziIncl.y), -cos(cAziIncl.y))  / cSceneScale);
  //cMat[0] = normalize(cross(cMat[2], vec3(cView / cSceneScale) + vec3(1e-16)));
  cMat[0] = vec3(cAzimuthVector.y, -cAzimuthVector.x, 0.0);
  cMat[1] = cross(cMat[2], cMat[0]);



  float
    rTessXFac;

  if (gl_TessCoord.y == 0.0f)
  {
    rTessXFac = gl_TessLevelOuter[1];
  }
  else
  {
    rTessXFac = gl_TessLevelOuter[3];
  }
  
  float
    rTessUnit = 1.0 / rTessXFac;

  rTessXFac =  rTessXFac / ((rTessXFac - 2.0) * 0.5);

  float
    rAngleFac = gl_TessCoord.x;

  // outer circle
  if (rAngleFac == 0.5)
  {
    gl_Position = vec4(modelViewProjectionMatrix * vec4(cCurrVertex.xyz, 1.0));
  }
  else
  {
    if (rAngleFac < 0.5)
    {
      rAngleFac *= rTessXFac;
    }
    else //if (rAngleFac >= 0.5)
    {
      rInnerOuterRadiusFac = 0.0;
      rAngleFac -= 0.5 + rTessUnit;
      rAngleFac *= rTessXFac;
      rAngleFac = 1.0 - rAngleFac;
    }

    float
      rInnerRadius = mix(cRadiusScale.x,cRadiusScale.y,rRadiusLog),
      rOuterRadius = mix(cRadiusScale.z,cRadiusScale.w,rRadiusLog);

    float
      rFac = 2.0 * rAngleFac * 3.14159265;

    rFac += cAziIncl.x;//REPLACE_CONSISTENT_ASMIUTH
    rFac -= 0.5 * 3.14159265;

    float
      rRadius = mix(rInnerRadius, rOuterRadius, rInnerOuterRadiusFac);
  
    vec2 cDir = vec2(sin(rFac), cos(rFac));
    vec3 cPos = cMat * vec3(cDir * rRadius, 0.0);

    cPos *= cSceneScale;
    cPos += cCurrVertex.xyz;

    vCenter = vec3(modelViewMatrix* vec4(cCurrVertex.xyz , 1.0));
    vDelta = vec3(modelViewMatrix * vec4(cPos.xyz, 1.0)) - vCenter;
    vLength = length(vDelta);

    gl_Position = modelViewProjectionMatrix * vec4(cPos, 1.0);
  }
}