Ok! here is my attempt with quaternion julia (fragmentarium script):

#info 4D Quaternion Julia Distance Estimator

#include "DE-Raytracer.frag"

#group 4D Quaternion Julia

// The DE implementation here is based on the implementation by Subblue:

// http://www.subblue.com/blog/2009/9/20/quaternion_julia

// which in turn is based on a CG shader by Keenan Crane.

// My implementation here is very compressed, so take a look at Subblue's above for a much clearer implementation.

//Knighty: modified the script to render orbit trap. The trap is just a cylinder

// Number of fractal iterations.

uniform int Iterations; slider[0,16,100]

// Breakout distance

uniform float Threshold; slider[0,10,100]

// Quaterion Constant

uniform vec4 C; slider[(-1,-1,-1,-1),(0.18,0.88,0.24,0.16),(1,1,1,1)]

uniform float val; slider[-2,0,2]

uniform float cylRad; slider[0,0.1,1]

uniform float cylHeight; slider[0,2,5]

float cylinder(vec3 p){

p=p.yxz;

p.x-=val;

return max(abs(p.z)-cylHeight,length(p.xy)-cylRad);

}

// The inline expanded quaterion multiplications make this DE

// look much worse than it actually is.

float DE(vec3 pos) {

float d=10000.0;

vec4 p = vec4(pos, 0.0);

vec4 dp = vec4(1.0, 0.0,0.0,0.0);

for (int i = 0; i < Iterations; i++) {

dp = 2.0* vec4(p.x*dp.x-dot(p.yzw, dp.yzw), p.x*dp.yzw+dp.x*p.yzw+cross(p.yzw, dp.yzw));

p = vec4(p.x*p.x-dot(p.yzw, p.yzw), vec3(2.0*p.x*p.yzw)) + C;

float p2 = dot(p,p);

orbitTrap = min(orbitTrap, abs(vec4(p.xyz,p2)));

float r = cylinder(p.xyz);

d=min(d,0.5 * r / length(dp));

if (p2 > Threshold) break;

}

// length(p);

return d;// 0.5 * r / length(dp);

}

#preset Default

FOV = 0.4

Eye = -3.92876,-4.64711,5.23499

Target = 1.01476,1.13689,-1.254

Up = -0.0539774,-0.724628,-0.687023

AntiAlias = 1

Detail = -3.01273

DetailAO = -0.5

FudgeFactor = 1

MaxRaySteps = 129

BoundingSphere = 3 Locked

Dither = 0.5 Locked

NormalBackStep = 1

AO = 0.529412,0.352941,0,0.7

Specular = 1.5

SpecularExp = 16

SpotLight = 1,1,1,0.38043

SpotLightDir = 0.1,0.1

CamLight = 1,1,1,1

CamLightMin = 0

Glow = 1,1,1,0.16667

GlowMax = 20

Fog = 0

HardShadow = 0

ShadowSoft = 12.9032

Reflection = 0

BaseColor = 1,1,1

OrbitStrength = 0.8

X = 0.5,0.6,0.6,0.7

Y = 1,0.6,0,0.4

Z = 0.8,0.78,1,0.5

R = 0.4,0.7,1,0.12

BackgroundColor = 1,0.666667,0

GradientBackground = 0.3

CycleColors = false

Cycles = 1.1

EnableFloor = false

FloorNormal = 0,0,0

FloorHeight = 0

FloorColor = 1,1,1

Iterations = 5

Threshold = 30

C = -0.68182,-0.63636,0,0

val = -0.144

cylRad = 0.31532

cylHeight = 5

#endpreset

#preset Round

FOV = 0.4

Eye = -1.49304,-1.07827,5.47996

Target = 0.940476,0.755427,-4.04451

Up = -0.118586,-0.968971,-0.21685

AntiAlias = 1

Detail = -3.01273

DetailAO = -0.5

FudgeFactor = 1

MaxRaySteps = 129

BoundingSphere = 3 Locked

Dither = 0.5 Locked

NormalBackStep = 1

AO = 0.529412,0.352941,0,0.7

Specular = 1.5

SpecularExp = 16

SpotLight = 1,1,1,0.38043

SpotLightDir = 0.1,0.1

CamLight = 1,1,1,1

CamLightMin = 0

Glow = 1,1,1,0.16667

GlowMax = 20

Fog = 0

HardShadow = 0

ShadowSoft = 12.9032

Reflection = 0

BaseColor = 1,1,1

OrbitStrength = 0.8

X = 0.5,0.6,0.6,0.7

Y = 1,0.6,0,0.4

Z = 0.8,0.78,1,0.5

R = 0.4,0.7,1,0.12

BackgroundColor = 1,0.666667,0

GradientBackground = 0.3

CycleColors = false

Cycles = 1.1

EnableFloor = false

FloorNormal = 0,0,0

FloorHeight = 0

FloorColor = 1,1,1

Iterations = 6

Threshold = 30

C = -0.68182,-0.63636,0,0

val = -1.84

cylRad = 0.69369

cylHeight = 0.0505

#endpreset