//
http://www.fractalforums.com/3d-fractal-generation/kaleidoscopic-%28escape-time-ifs%29/// Return reflection matrix for plane with normal 'n'
float3x3 reflectionMatrix(float3 n)
{
const float3 n1_x = (float3)(1 - 2 * n.x * n.x, -2 * n.y * n.x, -2 * n.z * n.x);
const float3 n1_y = (float3)(-2 * n.x * n.y, 1 - 2 * n.y * n.y, -2 * n.z * n.y);
const float3 n1_z = (float3)(-2 * n.x * n.z, -2 * n.y * n.z, 1 - 2 * n.z * n.z);
}
float KaleidoIcosaDE(const float3 z0, float t)
{
const float scale = 2.393331f;
const float phi = 1.6180339887498948482045868343656f;
const float3 offset = normalize((float3)(1, phi - 1, 0));
//const float3 n1 = normalize((float3)(-phi, phi - 1, 1));
//const float3 n2 = normalize((float3)(1, -phi, phi + 1));
//const float3 offset = normalize((float3)(1, phi - 1, 0));
//const float3 offset = normalize((float3)(1.01f, phi - 1 + 0.21f + m, -0.1f + q));
const float3 n1 = normalize((float3)(-0.8652f, 0.752f, 1.25f));
const float3 n2 = normalize((float3)(0.6152f, -0.23521f, 1.325f));
//n1Mat = reflectionMatrix(n1);
const float3 n1_x = (float3)(1 - 2 * n1.x * n1.x, -2 * n1.y * n1.x, -2 * n1.z * n1.x);
const float3 n1_y = (float3)(-2 * n1.x * n1.y, 1 - 2 * n1.y * n1.y, -2 * n1.z * n1.y);
const float3 n1_z = (float3)(-2 * n1.x * n1.z, -2 * n1.y * n1.z, 1 - 2 * n1.z * n1.z);
//n2Mat = reflectionMatrix(n2);
const float3 n2_x = (float3)(1 - 2 * n2.x * n2.x, -2 * n2.y * n2.x, -2 * n2.z * n2.x);
const float3 n2_y = (float3)(-2 * n2.x * n2.y, 1 - 2 * n2.y * n2.y, -2 * n2.z * n2.y);
const float3 n2_z = (float3)(-2 * n2.x * n2.z, -2 * n2.y * n2.z, 1 - 2 * n2.z * n2.z);
const float bailout2 = 1e9f; // 1e32f
const int max_iterations = 22; // 28
int n = 0;
float3 p = z0;
while (n < max_iterations)
{
p.y = fabs(p.y);
p.z = fabs(p.z);
if (dot(p, n2) > 0)
{
//p *= n2Mat;
p = (float3)(dot(p, n2_x), dot(p, n2_y), dot(p, n2_z));
}
p.x = fabs(p.x);
p.z = fabs(p.z);
if (dot(p, n1) > 0)
{
//p *= n1Mat;
p = (float3)(dot(p, n1_x), dot(p, n1_y), dot(p, n1_z));
}
p = p * scale - offset * (scale - 1);
if (dot(p, p) > bailout2) break;
n++;
}
return length(p) * pow(scale, (float)-n);
}