Thanks, here's the code:

Vector class:

class Vector(){

public:

func Init(float xx, float yy, float zz)

x = xx

y = yy

z = zz

endfunc

float func Dot(Vector other)

return x*other.x + y*other.y + z*other.z

endfunc

func Divide(float f)

x = x / f

y = y / f

z = z / f

endfunc

func Multiply(Vector a, float f)

x = a.x * f

y = a.y * f

z = a.z * f

endfunc

func MultiplyEquals(float f)

x = x * f

y = y * f

z = z * f

endfunc

func Normalise()

float mag = sqrt(x*x + y * y + z*z)

if mag>0

x = x / mag

y = y / mag

z = z / mag

endif

endfunc

float func Magnitude()

return sqrt(x*x + y*y + z*z)

endfunc

func Add(Vector a, Vector b)

x = a.x + b.x

y = a.y + b.y

z = a.z + b.z

endfunc

func AddEquals(Vector a)

x = x + a.x

y = y + a.y

z = z + a.z

endfunc

func Subtract(Vector a, Vector b)

x = a.x - b.x

y = a.y - b.y

z = a.z - b.z

endfunc

func Cross(Vector a, Vector b)

float _x = (a.y * b.z) - (a.z * b.y)

float _y = (a.z * b.x) - (a.x * b.z)

float _z = (a.x * b.y) - (a.y * b.x)

x = _x

y = _y

z = _z

endfunc

float x

float y

float z

}

Initialisation:

Vector corners[4] ; unit length tetrahedron corners

corners[0] = new Vector()

corners[1] = new Vector()

corners[2] = new Vector()

corners[3] = new Vector()

corners[0].Init(0, 1, 0)

corners[1].Init(0, -1/3, sqrt(8/9))

corners[2].Init(-sqrt(2/3), -1/3, -(1/3) * sqrt(2))

corners[3].Init(sqrt(2/3), -1/3, -(1/3) * sqrt(2))

Vector normals[4]

Vector crosses[4,3] ; crosses[face index][corner index]. cross product of corners for inside-outside measures

Vector C = new Vector()

int triangles[4,3] ; triangles[face index, corner of face]

triangles[0, 0] = 0 ; wish I knew how to do an initialiser list

triangles[0, 1] = 1

triangles[0, 2] = 2

triangles[1, 0] = 0

triangles[1, 1] = 2

triangles[1, 2] = 3

triangles[2, 0] = 0

triangles[2, 1] = 3

triangles[2, 2] = 1

triangles[3, 0] = 2

triangles[3, 1] = 1

triangles[3, 2] = 3

; for use later

Vector point = new Vector()

Vector corn0 = new Vector()

Vector corn1 = new Vector()

Vector corn2 = new Vector()

Vector bary = new Vector()

int it = 0

repeat

normals[it] = new Vector()

crosses[it, 0] = new Vector()

crosses[it, 1] = new Vector()

crosses[it, 2] = new Vector()

Vector cr1 = new Vector()

cr1.Subtract(corners[triangles[it, 2]], corners[triangles[it, 0]])

Vector cr2 = new Vector()

cr2.Subtract(corners[triangles[it, 1]], corners[triangles[it, 0]])

normals[it].Cross(cr1, cr2)

normals[it].Normalise()

; leave this divide out below for an alternative where face of tetrahedron has magnitude 1, rather than corners. Produces slightly different shape

normals[it].Divide(normals[it].Dot(corners[triangles[it, 0]])) ; normal is not unit length so that a corner point gives dot of 1

crosses[it, 0].Cross(corners[triangles[it, 2]], corners[triangles[it, 1]])

crosses[it, 1].Cross(corners[triangles[it, 0]], corners[triangles[it, 2]])

crosses[it, 2].Cross(corners[triangles[it, 1]], corners[triangles[it, 0]])

until (it=it+1)>=4

int squareMap[4,4,3] ; squareMap[quadrant, face index, corner of quadrant]

squareMap[0, 0, 0] = 0

squareMap[0, 0, 1] = 1

squareMap[0, 0, 2] = 2

squareMap[0, 1, 0] = 0

squareMap[0, 1, 1] = 2

squareMap[0, 1, 2] = 3

squareMap[0, 2, 0] = 0

squareMap[0, 2, 1] = 3

squareMap[0, 2, 2] = 1

squareMap[0, 3, 0] = 0

squareMap[0, 3, 1] = 3

squareMap[0, 3, 2] = 1

squareMap[1, 0, 0] = 1

squareMap[1, 0, 1] = 0

squareMap[1, 0, 2] = 3

squareMap[1, 1, 0] = 2

squareMap[1, 1, 1] = 0

squareMap[1, 1, 2] = 1

squareMap[1, 2, 0] = 3

squareMap[1, 2, 1] = 0

squareMap[1, 2, 2] = 2

squareMap[1, 3, 0] = 3

squareMap[1, 3, 1] = 0

squareMap[1, 3, 2] = 2

squareMap[2, 0, 0] = 2

squareMap[2, 0, 1] = 3

squareMap[2, 0, 2] = 0

squareMap[2, 1, 0] = 3

squareMap[2, 1, 1] = 1

squareMap[2, 1, 2] = 0

squareMap[2, 2, 0] = 1

squareMap[2, 2, 1] = 2

squareMap[2, 2, 2] = 0

squareMap[2, 3, 0] = 1

squareMap[2, 3, 1] = 2

squareMap[2, 3, 2] = 0

squareMap[3, 0, 0] = 1

squareMap[3, 0, 1] = 3

squareMap[3, 0, 2] = 2

squareMap[3, 1, 0] = 2

squareMap[3, 1, 1] = 1

squareMap[3, 1, 2] = 3

squareMap[3, 2, 0] = 3

squareMap[3, 2, 1] = 2

squareMap[3, 2, 2] = 1

squareMap[3, 3, 0] = 3

squareMap[3, 3, 1] = 2

squareMap[3, 3, 2] = 1

Inner loop:

; convert c and z to vectors, for convenience

C.x = real(cri)

C.y = imag(cri)

C.z = cj

point.x = real(zri)

point.y = imag(zri)

point.z = zj

float magnitude = 1

int ij = 0 ; ij is just a unique name for the iterator over each tetrahedron face

repeat

magnitude = point.Dot(normals[ij]) ; normals isn't unit length so this gives radial magnitude in shape of tetrahedron

if (magnitude > 0) ; on right side of triangle

; TDL super simple way to get barycentric coordinates of triangle

bary.x = crosses[ij, 0].Dot(point)

bary.y = crosses[ij, 1].Dot(point)

bary.z = crosses[ij, 2].Dot(point)

if (bary.x >= 0 && bary.y >= 0 && bary.z >= 0) ; correct face is found

bary.Divide(bary.x+bary.y+bary.z) ; barycentrics sum to 1

int n = 0

if (bary.x >= 0.5) ; quadrant 0

n = 0

bary.x = bary.x - 0.5

elseif (bary.y >= 0.5) ; quadrant 1

n = 1

bary.y = bary.y - 0.5

elseif (bary.x + bary.y <= 0.5) ; quadrant 2

n = 2

else ; quadrant 3 (middle)

n = 3

float baryz = 1 - (bary.x + bary.y)

bary.y = 0.5 - bary.x

bary.x = 0.5 - baryz

endif

; Z^2 + c

bary.x = bary.x * 2 ; double the size here

bary.y = bary.y * 2

bary.z = 1 - (bary.x + bary.y)

magnitude = magnitude*magnitude;

corn0.Multiply(corners[squareMap[n, ij, 0]], bary.x)

corn1.Multiply(corners[squareMap[n, ij, 1]], bary.y)

corn2.Multiply(corners[squareMap[n, ij, 2]], bary.z)

point.Add(corn0, corn1)

point.AddEquals(corn2)

point.MultiplyEquals(magnitude)

point.AddEquals(C)

ij=10 ; break out of loop since the face was found

endif

endif

until (ij = ij+1)>=4

zri = point.x + flip(point.y) ; convert back to z

zj = point.z