float quartic_evaluate( float t, float c, float p, float h, float k )
{
   float a4 = 1.0-k*k;
   float a3 = -2.0*p*(1.0-k*k);
   float a2 = h*h + (1.0-k*k) * p*p - k*k*(c-h)*(c-h);
   float a1 = 2.0*k*k*(c-h)*(c-h)*p;
   float a0 = -k*k*(c-h)*(c-h)*p*p;
   return a4*t*t*t*t + a3*t*t*t + a2*t*t + a1*t + a0;
}

// Per vertex transformation to account for true refraction underwater
vec3 refraction_solve( vec3 world_co, vec3 camera_co, float water_height )
{
   if( world_co.y > water_height )
      return world_co;

   float c = camera_co.y - world_co.y;
   float p = length( world_co.xz - camera_co.xz );
   float h = water_height - world_co.y;
   float k = 1.3333/1.0005;

   float ta = ((c-h)/c) * p;
   float tb = p;

   float Fa = quartic_evaluate( ta, c,p,h,k );
   float Fb = quartic_evaluate( tb, c,p,h,k );

   float tc = 0.0;

   for( int i=0; i<10; i ++ )
   {
      tc = (0.5*Fb*ta - Fa*tb) / (0.5*Fb - Fa);
      float Fc = quartic_evaluate( tc, c,p,h,k );

      if( Fc < 0.0 )
      {
         Fa = Fc;
         ta = tc;
      }
      else
      {
         Fb = Fc;
         tb = tc;
      }
   }

   float j = c - (p / (tc / (c-h)));
   return world_co + vec3(0,j,0);
}