float3 Hit( float3 x, float3 R, samplerCUBE mp )
{
  
	float rl = readDistanceCubeMap(mp, R); // |r|
	
	float dp = rl - dot(x, R);		// parallax
	
	float3 p = x + R * dp;
	return p;
}
---------------------------------------------------------------------------------------

// This function is called several times.
float3 Hit( float3 x, float3 R, samplerCUBE mp )
{
	float rl = readDistanceCubeMap( mp, R);				// |r|
	
	float ppp = length( x ) /  readDistanceCubeMap( mp, x);			// |p|/|p’|
	float dun = 0, pun = ppp, dov = 0, pov = 0;
	float dl = rl * ( 1 - ppp );							// eq. 2
	float3 l = x + R * dl;									// ray equation

	// iteration
	for( int i = 0; i < 2; i++ )	// 2 !!!!!!!!!!!!!!!!!!!!!!!
	{
		float llp = length( l ) / readDistanceCubeMap( mp, l);		// |l|/|l’|
		if ( llp < 0.999f )									// undershooting
		{
			dun = dl; pun = llp;							// last undershooting
			dl += ( dov == 0 ) ? rl * ( 1 - llp ) :			// eq. 2
				( dl - dov ) * ( 1 - llp ) / ( llp - pov );	// eq. 3
		} else if ( llp > 1.001f )							// overshooting
		{
			dov = dl; pov = llp;							// last overshooting
			dl += ( dl -dun ) * ( 1 - llp ) / ( llp - pun );// eq. 3
		}
		l = x + R * dl;										// ray equation
	}
	return l;												// computed hit point
}
----------------------------------------------------------------------------------------

float3 Hit(float3 x, float3 R, samplerCUBE mp)
{
	R = normalize(R);
	float3 xNorm = normalize(x);
	float3 dt = (R - xNorm) / (float) LIN_ITERATIONCOUNT;
	float dx = length(cross(R,x));
	
	float3 pN = normalize(cross(R, xNorm));
	float3 N = normalize(cross(pN,R));
	float d = dot(x, N);	
	
	float3 under = xNorm * readDistanceCubeMap(mp, x);
	float3 over = R * readDistanceCubeMap(mp, R);
	
	//linear iteration	
	for(int i = 1; i < LIN_ITERATIONCOUNT; i++)
	{
		float3 dir = normalize(xNorm + dt * i);
		float dist = readDistanceCubeMap( mp, dir);
		float3 point = dir * dist;
				
		if( dot(N, point) > d)//undershooting
		{
			under = point;			
		}
		else				//overshooting
		{
			over = point;			
			i = LIN_ITERATIONCOUNT;
		}	
	}
	
	float3 dirUn = normalize(under);
	float3 dirOv = normalize(over);
	float tun = d / dot(N, dirUn);
	float pun = tun / length(under);
	float dun = length(tun * dirUn - x);
	float tov = d / dot(N, dirOv);
	float pov = tov / length(over);
	float dov = length(tov * dirOv - x);
	float tl = (d - dot(N, over)) / dot(N, under - over);
	float3 l = over + tl * (under - over);
	float dl = length(l - x);
		
	//secant iteration
	for( int i = 0; i < SECANT_ITERATIONCOUNT; i++ )
	{
		float llp = length( l ) / readDistanceCubeMap( mp, l);		// |l|/|l’|
		if ( llp < 0.999f )									// undershooting
		{
			dun = dl; pun = llp;							// last undershooting
			dl += ( dl - dov ) * ( 1 - llp ) / ( llp - pov );	// eq. 3
		} else if ( llp > 1.001f )							// overshooting
		{
			dov = dl; pov = llp;							// last overshooting
			dl += ( dl - dun ) * ( 1 - llp ) / ( llp - pun );// eq. 3
		}
		l = x + R * dl;										// ray equation
	}
	
	return l;		
	return over;	
	return (under + over)/2.0;
}
----------------------------------------------------------------------------------------------


float3 Hit(float3 x, float3 R, samplerCUBE mp)
{
	//return R + x * 0.0000000001;
	
	R = normalize(R);
	float3 xNorm = normalize(x);
	float3 dt = (R - xNorm) / (float) LIN_ITERATIONCOUNT;
	float dx = length(cross(R,x));
	
	float3 pN = normalize(cross(R, xNorm));
	float3 N = normalize(cross(pN,R));
	float D = dot(x, N);	
	
	float under = readDistanceCubeMap(mp, x);
	float over = readDistanceCubeMap(mp, R);
	float dun = 0;
	float dov = 0;
	float lun = length(x);
	float lov;
	
	float3 l;
	
	//linear iteration	
	for(int i = 1; i < LIN_ITERATIONCOUNT; i++)
	{
	
		l = normalize(xNorm + dt * i);
		float lL = readDistanceCubeMap( mp, l);
		float lambda = D / dot(N, l);
		//float lambda = (dot(x,l) - dot(R,l) * dot(x,R)) / (lL * lL - dot(R,l) * dot(R,l));
		float d = dot(l, R) * lambda - dot(x, R);
					
		if( lambda < lL)//undershooting
		{
			lun = lambda;
			dun = d;
			under = lL;			
		}
		else				//overshooting
		{
			lov = lambda;
			dov = d;
			over = lL;			
			i = LIN_ITERATIONCOUNT;
		}	
	}
	
	if(dov == 0)l = R;
	else l = x + R * dov;
	return l;
	
		
	float pun = lun / length(under);
	float pov;
	float dl;
	
	if(dov == 0)
	{
		dl = over * ( 1 - pun );
	}
	else
	{
		pov = lov / length(over);	
		dl = dun + (dun - dov) * (1 - pun)/(pun - pov);		
	}
	l = x + R * dl;	
		
	//secant iteration
	for( int i = 0; i < SECANT_ITERATIONCOUNT; i++ )
	{
		float llp = length( l ) / readDistanceCubeMap( mp, l);		// |l|/|l’|
		if ( llp < 0.999f )									// undershooting
		{
			dun = dl; pun = llp;							// last undershooting
			dl += ( dov == 0 ) ? over * ( 1 - llp ) :			// eq. 2
				( dl - dov ) * ( 1 - llp ) / ( llp - pov );	// eq. 3
		} else if ( llp > 1.001f )							// overshooting
		{
			dov = dl; pov = llp;							// last overshooting
			dl += ( dl -dun ) * ( 1 - llp ) / ( llp - pun );// eq. 3
		}
		l = x + R * dl;										// ray equation
	}
	
	return l;		
	return over;	
	return (under + over) / 2.0;
}