uniform vec2 cNearFarAB;

uniform sampler3D EarthTexture;

float
noise(vec3 cPos)
{
  int nOctaves = 5;
  vec3  cIntPos;
  vec3  cDist;
  vec3  acStep[3];

  float rSkew =   1.0 / 3.0;
  float rUnSkew = 1.0 / 6.0;

  float rNoise      = 0.0;
  float rValueScale = 1.0;
  float rTotal      = 0.0;
  float dist[4];

  for(int iOctave = 0; iOctave < nOctaves; iOctave++)
  {
    float rSum = dot(cPos, vec3(1.0));
    cIntPos = floor(cPos + vec3(rSum * rSkew));
    float rIntSum = dot(cIntPos, vec3(1.0));
    cDist = cPos - (cIntPos - vec3(rIntSum * rUnSkew));

    if(cDist.x > cDist.y)
    {
      // y cannot be the largest
      if(cDist.x > cDist.z)
      {
        acStep[0] = vec3(1.0, 0.0, 0.0);
      }
      else
      {
        acStep[0] = vec3(0.0, 0.0, 1.0);
      }
      // x cannot be the smallest
      if(cDist.z > cDist.y)
      {
        acStep[2] = vec3(0.0, 1.0, 0.0);
      }
      else
      {
        acStep[2] = vec3(0.0, 0.0, 1.0);
      }
    }
    else
    {
      // x cannot be the largest
      if(cDist.y > cDist.z)
      {
        acStep[0] = vec3(0.0, 1.0, 0.0);
      }
      else
      {
        acStep[0] = vec3(0.0, 0.0, 1.0);
      }
      // y cannot be the smallest
      if(cDist.z > cDist.x)
      {
        acStep[2] = vec3(1.0, 0.0, 0.0);
      }
      else
      {
        acStep[2] = vec3(0.0, 0.0, 1.0);
      }
    }
    acStep[1] = vec3(1.0, 1.0, 1.0) - acStep[0] - acStep[2];

    // The 3D simplex has 4 corners
    for(int iCorner = 0; iCorner < 4; iCorner++)
    {
      dist[iCorner] = length(cDist);
      float rWeight = 0.6 - length(cDist); // This corner contributes relative to the distance of the point from this corner

      if(rWeight > 0.0)
      {
        rWeight *= rWeight;
        rWeight *= rWeight;
        
        vec3 cGradient = texture3D(EarthTexture, cIntPos * (1.0 / 128.0) + vec3(0.5 / 128.0)).xyz;
        cGradient = (cGradient - vec3(0.5)) * 2.0;

        float rDot = dot(cDist, cGradient);
        rNoise += rDot * rWeight * rValueScale;
      }
      
      if(iCorner < 3)
      {
        cDist   += vec3(rUnSkew);
        cDist   -= acStep[iCorner];
        cIntPos += acStep[iCorner];
      }
    }
    rTotal += rValueScale;
    rValueScale *= 0.35;
    cPos *= 1.918281828459;
  }
  
  rNoise = 0.5 + rNoise * 64.0 / rTotal;

  return rNoise;
}

vec3
calcNormal(vec3 N, float rValue0, float rValue1, float rValue2)
{
  vec3 cDelta01 = normalize(cross(N, vec3(0.001, 1.0, 0.001)));
  vec3 cDelta02 = cross(N, cDelta01);
  
  cDelta01 -= N * (rValue1 - rValue0);
  cDelta02 -= N * (rValue2 - rValue0);
  return normalize(cross(cDelta01, cDelta02));
}

float GetRealZ(float rZ)
{
  float dist = cNearFarAB.y / (rZ - cNearFarAB.x);
  return dist;
}

float GetInvZ(float dist)
{
  float rZ = cNearFarAB.y / dist + cNearFarAB.x;
  return rZ;
}

float TransformZ(float rDeltaZ)
{
  return GetInvZ(GetRealZ(gl_FragCoord.z) - rDeltaZ);
}

 vec2 get_dist(vec3 uvw)
{
   vec3 p = floor(uvw);
   vec2 min_dq = vec2(999, 999);
   for (int i = -1; i < 2; ++i)
      for (int j = -1; j < 2; ++j)
         for (int k = -1; k < 2; ++k) {
            vec3 p2 = p + vec3(ivec3(i, j, k));
            p2.x += (texture3D(EarthTexture, p2.xyz / 128.0) + 0.1).x * 0.8;
            p2.y += (texture3D(EarthTexture, p2.xyz / 128.0 + vec3(0.01,0,0)) + 0.1).x * 0.8;
            p2.z += (texture3D(EarthTexture, p2.xyz / 128.0 + vec3(0.02,0,0)) + 0.1).x * 0.8;
            float dq = dot(p2 - uvw, p2 - uvw);
            if (dq < min_dq.x) {
               min_dq.y = min_dq.x;
               min_dq.x = dq;
            } else if (dq < min_dq[1])
               min_dq.y = dq;
         }
   return sqrt(min_dq) / 4.0;
}

float get_height(vec3 uvw)
{
   vec2 v = get_dist(uvw * 4.0);
   return clamp(v.y - v.x, 0.0, 1.0);
}


vec2 iGetIntegerAndFractional(float x)
{
  vec2
    cVec;

  cVec.x = floor(x);
  cVec.y = x - cVec.x;

  return cVec;
}

float myRandomMagic(int i, int j, int s)
{
  vec3
    cIntPos;
  
  cIntPos = vec3(i, j, s);

  return texture3D(EarthTexture, cIntPos * (1.0 / 128.0) + vec3(0.5 / 128.0)).x;
}

vec4 dnoise3f(const float x, const float y, const float z )
{
    int   i, j, k;
    float u, v, w;

    vec2
      cIU = iGetIntegerAndFractional(x),
      cJV = iGetIntegerAndFractional(y),
      cKW = iGetIntegerAndFractional(z);

    i = int(cIU.x);
    j = int(cJV.x);
    k = int(cKW.x);

    u = cIU.y;
    v = cJV.y;
    w = cKW.y;
    
    float du = 30.0*u*u*(u*(u-2.0)+1.0);
    float dv = 30.0*v*v*(v*(v-2.0)+1.0);
    float dw = 30.0*w*w*(w*(w-2.0)+1.0);

    u = u*u*u*(u*(u*6.0-15.0)+10.0);
    v = v*v*v*(v*(v*6.0-15.0)+10.0);
    w = w*w*w*(w*(w*6.0-15.0)+10.0);

    float a = myRandomMagic( i+0, j+0, k+0 );
    float b = myRandomMagic( i+1, j+0, k+0 );
    float c = myRandomMagic( i+0, j+1, k+0 );
    float d = myRandomMagic( i+1, j+1, k+0 );
    float e = myRandomMagic( i+0, j+0, k+1 );
    float f = myRandomMagic( i+1, j+0, k+1 );
    float g = myRandomMagic( i+0, j+1, k+1 );
    float h = myRandomMagic( i+1, j+1, k+1 );

    float k0 =   a;
    float k1 =   b - a;
    float k2 =   c - a;
    float k3 =   e - a;
    float k4 =   a - b - c + d;
    float k5 =   a - c - e + g;
    float k6 =   a - b - e + f;
    float k7 = - a + b + c - d + e - f - g + h;

    vec4
      cReturn;

    cReturn = vec4(k0 + k1*u + k2*v + k3*w + k4*u*v + k5*v*w + k6*w*u + k7*u*v*w,
                   du * (k1 + k4*v + k6*w + k7*v*w),
                   dv * (k2 + k5*w + k4*u + k7*w*u),
                   dw * (k3 + k6*u + k5*v + k7*u*v));

    return cReturn;
}

float
fbm(vec3 cPos)
{
  float
    x = cPos.x;
  float
    y = cPos.y;
  float
    z = cPos.z;

    float f = 0.0;
    float w = 0.5;
    float dx = 0.0;
    float dy = 0.0;
    float dz = 0.0;
    for( int i=0; i < 8 ; i++ )
    {
      vec4
        cResult = dnoise3f(x, y, z );
        dx += cResult.y;
        dy += cResult.z;
        dz += cResult.w;
        f += w * cResult.x / (1.0 + dx*dx + dy*dy + dz*dz); // replace with "w * n[0]" for a classic fbm()
	w *= 0.5;
	x *= 2.0;
	y *= 2.0;
	z *= 2.0;
    }
    
  return f;
}

float
ffunc(vec3 cPos, float rFac, bool isMulColor)
{
  float
    rU = rFac;

  rU = (rU - cScalarVertexValueRangeAndNoValue.x) * cScalarVertexValueRangeAndNoValue.y;

  vec3
    cScaledPos = cPos * vec3(100.0, 100.0, 100.0);
    //cScaledPos = cPos * vec3(0.02,0.02,4.0);

  float f = fbm(cScaledPos * vec3(0.7, 0.7, 0.8));
 

  float
    rMaterial = f / 0.7;

  rMaterial *= 0.04;
  rU *= 0.96;

  rMaterial += rU;

  if (isMulColor)
  {
    gl_FragColor.rgb *= texture1D(textureScalarVertex, rMaterial).rgb;
  }

  f-=rU * 0.6;
  if (f<0)
  {
    f *= -0.1;
  }
  else
  {
    f *= -0.2;
  }

  rU = smoothstep(0,1,rMaterial);
  rU *= rU;

  float
    rSand = noise(cScaledPos * vec3(7,7,4.0) + vec3(0,0,f));//texture3D(EarthTexture, fmod(cScaledPos * 0.5,1.0));

  rSand *= rSand;
  rSand -= 0.21 + rMaterial * 0.1;
  rSand *= 2.0;
  if (rSand > 0.0) rSand = 0.0;

  rU -= 0.25;
  rU *= 0.4;
  rU = min(max(rU,0.0),1.0);

  f = (rSand - f) * rU + f;

  float
    rCrack = 1.0 - texture1D(textureScalarVertex, rMaterial).a;
  
  f -= pow(rCrack, 1.5) * 0.35;

  return f * 10.0;

  /*hape
    rSand = texture3D(EarthTexture, fmod(cScaledPos * 2.0,1.0));

  rSand = pow(rSand, 4.0) * -0.2;

  f = (rSand - f) * rU + f;

  float
    rCrack = 1.0 - texture1D(textureScalarVertex, rU).a;
  
  f -= pow(rCrack, 8.0);

  return f * 6.0;
  */
}


