//

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);

}