text:
function overlap_cuboid_binary(	xAB::Float64,
							yAB::Float64,
							zAB::Float64,
							a11::Float64,
							a12::Float64,
							a13::Float64,
							a21::Float64,
							a22::Float64,
							a23::Float64,
							a31::Float64,
							a32::Float64,
							a33::Float64,
							b11::Float64,
							b12::Float64,
							b13::Float64,
							b21::Float64,
							b22::Float64,
							b23::Float64,
							b31::Float64,
							b32::Float64,
							b33::Float64,
							r1a::Float64,
							r2a::Float64,
							r3a::Float64,
							r1b::Float64,
							r2b::Float64,
							r3b::Float64)

	overlapfun::Float64 = 1.0

	cA::Array{Float64, 1} = zeros(8)
	cB::Array{Float64, 1} = zeros(8)

	cA = [ - a11*(a11*r1a + a12*r2a + a13*r3a) - a21*(a21*r1a + a22*r2a + a23*r3a) - a31*(a31*r1a + a32*r2a + a33*r3a), - a11*(a11*r1a + a12*r2a - a13*r3a) - a21*(a21*r1a + a22*r2a - a23*r3a) - a31*(a31*r1a + a32*r2a - a33*r3a), - a11*(a11*r1a - a12*r2a + a13*r3a) - a21*(a21*r1a - a22*r2a + a23*r3a) - a31*(a31*r1a - a32*r2a + a33*r3a), a11*(a12*r2a - a11*r1a + a13*r3a) + a21*(a22*r2a - a21*r1a + a23*r3a) + a31*(a32*r2a - a31*r1a + a33*r3a), - a11*(a12*r2a - a11*r1a + a13*r3a) - a21*(a22*r2a - a21*r1a + a23*r3a) - a31*(a32*r2a - a31*r1a + a33*r3a), a11*(a11*r1a - a12*r2a + a13*r3a) + a21*(a21*r1a - a22*r2a + a23*r3a) + a31*(a31*r1a - a32*r2a + a33*r3a), a11*(a11*r1a + a12*r2a - a13*r3a) + a21*(a21*r1a + a22*r2a - a23*r3a) + a31*(a31*r1a + a32*r2a - a33*r3a), a11*(a11*r1a + a12*r2a + a13*r3a) + a21*(a21*r1a + a22*r2a + a23*r3a) + a31*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - a11*(b11*r1b - xAB + b12*r2b + b13*r3b) - a21*(b21*r1b - yAB + b22*r2b + b23*r3b) - a31*(b31*r1b - zAB + b32*r2b + b33*r3b), a11*(xAB - b11*r1b - b12*r2b + b13*r3b) + a21*(yAB - b21*r1b - b22*r2b + b23*r3b) + a31*(zAB - b31*r1b - b32*r2b + b33*r3b), a11*(xAB - b11*r1b + b12*r2b - b13*r3b) + a21*(yAB - b21*r1b + b22*r2b - b23*r3b) + a31*(zAB - b31*r1b + b32*r2b - b33*r3b), a11*(xAB - b11*r1b + b12*r2b + b13*r3b) + a21*(yAB - b21*r1b + b22*r2b + b23*r3b) + a31*(zAB - b31*r1b + b32*r2b + b33*r3b), a11*(xAB + b11*r1b - b12*r2b - b13*r3b) + a21*(yAB + b21*r1b - b22*r2b - b23*r3b) + a31*(zAB + b31*r1b - b32*r2b - b33*r3b), a11*(xAB + b11*r1b - b12*r2b + b13*r3b) + a21*(yAB + b21*r1b - b22*r2b + b23*r3b) + a31*(zAB + b31*r1b - b32*r2b + b33*r3b), a11*(xAB + b11*r1b + b12*r2b - b13*r3b) + a21*(yAB + b21*r1b + b22*r2b - b23*r3b) + a31*(zAB + b31*r1b + b32*r2b - b33*r3b), a11*(xAB + b11*r1b + b12*r2b + b13*r3b) + a21*(yAB + b21*r1b + b22*r2b + b23*r3b) + a31*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ - a12*(a11*r1a + a12*r2a + a13*r3a) - a22*(a21*r1a + a22*r2a + a23*r3a) - a32*(a31*r1a + a32*r2a + a33*r3a), - a12*(a11*r1a + a12*r2a - a13*r3a) - a22*(a21*r1a + a22*r2a - a23*r3a) - a32*(a31*r1a + a32*r2a - a33*r3a), - a12*(a11*r1a - a12*r2a + a13*r3a) - a22*(a21*r1a - a22*r2a + a23*r3a) - a32*(a31*r1a - a32*r2a + a33*r3a), a12*(a12*r2a - a11*r1a + a13*r3a) + a22*(a22*r2a - a21*r1a + a23*r3a) + a32*(a32*r2a - a31*r1a + a33*r3a), - a12*(a12*r2a - a11*r1a + a13*r3a) - a22*(a22*r2a - a21*r1a + a23*r3a) - a32*(a32*r2a - a31*r1a + a33*r3a), a12*(a11*r1a - a12*r2a + a13*r3a) + a22*(a21*r1a - a22*r2a + a23*r3a) + a32*(a31*r1a - a32*r2a + a33*r3a), a12*(a11*r1a + a12*r2a - a13*r3a) + a22*(a21*r1a + a22*r2a - a23*r3a) + a32*(a31*r1a + a32*r2a - a33*r3a), a12*(a11*r1a + a12*r2a + a13*r3a) + a22*(a21*r1a + a22*r2a + a23*r3a) + a32*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - a12*(b11*r1b - xAB + b12*r2b + b13*r3b) - a22*(b21*r1b - yAB + b22*r2b + b23*r3b) - a32*(b31*r1b - zAB + b32*r2b + b33*r3b), a12*(xAB - b11*r1b - b12*r2b + b13*r3b) + a22*(yAB - b21*r1b - b22*r2b + b23*r3b) + a32*(zAB - b31*r1b - b32*r2b + b33*r3b), a12*(xAB - b11*r1b + b12*r2b - b13*r3b) + a22*(yAB - b21*r1b + b22*r2b - b23*r3b) + a32*(zAB - b31*r1b + b32*r2b - b33*r3b), a12*(xAB - b11*r1b + b12*r2b + b13*r3b) + a22*(yAB - b21*r1b + b22*r2b + b23*r3b) + a32*(zAB - b31*r1b + b32*r2b + b33*r3b), a12*(xAB + b11*r1b - b12*r2b - b13*r3b) + a22*(yAB + b21*r1b - b22*r2b - b23*r3b) + a32*(zAB + b31*r1b - b32*r2b - b33*r3b), a12*(xAB + b11*r1b - b12*r2b + b13*r3b) + a22*(yAB + b21*r1b - b22*r2b + b23*r3b) + a32*(zAB + b31*r1b - b32*r2b + b33*r3b), a12*(xAB + b11*r1b + b12*r2b - b13*r3b) + a22*(yAB + b21*r1b + b22*r2b - b23*r3b) + a32*(zAB + b31*r1b + b32*r2b - b33*r3b), a12*(xAB + b11*r1b + b12*r2b + b13*r3b) + a22*(yAB + b21*r1b + b22*r2b + b23*r3b) + a32*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ - a13*(a11*r1a + a12*r2a + a13*r3a) - a23*(a21*r1a + a22*r2a + a23*r3a) - a33*(a31*r1a + a32*r2a + a33*r3a), - a13*(a11*r1a + a12*r2a - a13*r3a) - a23*(a21*r1a + a22*r2a - a23*r3a) - a33*(a31*r1a + a32*r2a - a33*r3a), - a13*(a11*r1a - a12*r2a + a13*r3a) - a23*(a21*r1a - a22*r2a + a23*r3a) - a33*(a31*r1a - a32*r2a + a33*r3a), a13*(a12*r2a - a11*r1a + a13*r3a) + a23*(a22*r2a - a21*r1a + a23*r3a) + a33*(a32*r2a - a31*r1a + a33*r3a), - a13*(a12*r2a - a11*r1a + a13*r3a) - a23*(a22*r2a - a21*r1a + a23*r3a) - a33*(a32*r2a - a31*r1a + a33*r3a), a13*(a11*r1a - a12*r2a + a13*r3a) + a23*(a21*r1a - a22*r2a + a23*r3a) + a33*(a31*r1a - a32*r2a + a33*r3a), a13*(a11*r1a + a12*r2a - a13*r3a) + a23*(a21*r1a + a22*r2a - a23*r3a) + a33*(a31*r1a + a32*r2a - a33*r3a), a13*(a11*r1a + a12*r2a + a13*r3a) + a23*(a21*r1a + a22*r2a + a23*r3a) + a33*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - a13*(b11*r1b - xAB + b12*r2b + b13*r3b) - a23*(b21*r1b - yAB + b22*r2b + b23*r3b) - a33*(b31*r1b - zAB + b32*r2b + b33*r3b), a13*(xAB - b11*r1b - b12*r2b + b13*r3b) + a23*(yAB - b21*r1b - b22*r2b + b23*r3b) + a33*(zAB - b31*r1b - b32*r2b + b33*r3b), a13*(xAB - b11*r1b + b12*r2b - b13*r3b) + a23*(yAB - b21*r1b + b22*r2b - b23*r3b) + a33*(zAB - b31*r1b + b32*r2b - b33*r3b), a13*(xAB - b11*r1b + b12*r2b + b13*r3b) + a23*(yAB - b21*r1b + b22*r2b + b23*r3b) + a33*(zAB - b31*r1b + b32*r2b + b33*r3b), a13*(xAB + b11*r1b - b12*r2b - b13*r3b) + a23*(yAB + b21*r1b - b22*r2b - b23*r3b) + a33*(zAB + b31*r1b - b32*r2b - b33*r3b), a13*(xAB + b11*r1b - b12*r2b + b13*r3b) + a23*(yAB + b21*r1b - b22*r2b + b23*r3b) + a33*(zAB + b31*r1b - b32*r2b + b33*r3b), a13*(xAB + b11*r1b + b12*r2b - b13*r3b) + a23*(yAB + b21*r1b + b22*r2b - b23*r3b) + a33*(zAB + b31*r1b + b32*r2b - b33*r3b), a13*(xAB + b11*r1b + b12*r2b + b13*r3b) + a23*(yAB + b21*r1b + b22*r2b + b23*r3b) + a33*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ - b11*(a11*r1a + a12*r2a + a13*r3a) - b21*(a21*r1a + a22*r2a + a23*r3a) - b31*(a31*r1a + a32*r2a + a33*r3a), - b11*(a11*r1a + a12*r2a - a13*r3a) - b21*(a21*r1a + a22*r2a - a23*r3a) - b31*(a31*r1a + a32*r2a - a33*r3a), - b11*(a11*r1a - a12*r2a + a13*r3a) - b21*(a21*r1a - a22*r2a + a23*r3a) - b31*(a31*r1a - a32*r2a + a33*r3a), b11*(a12*r2a - a11*r1a + a13*r3a) + b21*(a22*r2a - a21*r1a + a23*r3a) + b31*(a32*r2a - a31*r1a + a33*r3a), - b11*(a12*r2a - a11*r1a + a13*r3a) - b21*(a22*r2a - a21*r1a + a23*r3a) - b31*(a32*r2a - a31*r1a + a33*r3a), b11*(a11*r1a - a12*r2a + a13*r3a) + b21*(a21*r1a - a22*r2a + a23*r3a) + b31*(a31*r1a - a32*r2a + a33*r3a), b11*(a11*r1a + a12*r2a - a13*r3a) + b21*(a21*r1a + a22*r2a - a23*r3a) + b31*(a31*r1a + a32*r2a - a33*r3a), b11*(a11*r1a + a12*r2a + a13*r3a) + b21*(a21*r1a + a22*r2a + a23*r3a) + b31*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - b11*(b11*r1b - xAB + b12*r2b + b13*r3b) - b21*(b21*r1b - yAB + b22*r2b + b23*r3b) - b31*(b31*r1b - zAB + b32*r2b + b33*r3b), b11*(xAB - b11*r1b - b12*r2b + b13*r3b) + b21*(yAB - b21*r1b - b22*r2b + b23*r3b) + b31*(zAB - b31*r1b - b32*r2b + b33*r3b), b11*(xAB - b11*r1b + b12*r2b - b13*r3b) + b21*(yAB - b21*r1b + b22*r2b - b23*r3b) + b31*(zAB - b31*r1b + b32*r2b - b33*r3b), b11*(xAB - b11*r1b + b12*r2b + b13*r3b) + b21*(yAB - b21*r1b + b22*r2b + b23*r3b) + b31*(zAB - b31*r1b + b32*r2b + b33*r3b), b11*(xAB + b11*r1b - b12*r2b - b13*r3b) + b21*(yAB + b21*r1b - b22*r2b - b23*r3b) + b31*(zAB + b31*r1b - b32*r2b - b33*r3b), b11*(xAB + b11*r1b - b12*r2b + b13*r3b) + b21*(yAB + b21*r1b - b22*r2b + b23*r3b) + b31*(zAB + b31*r1b - b32*r2b + b33*r3b), b11*(xAB + b11*r1b + b12*r2b - b13*r3b) + b21*(yAB + b21*r1b + b22*r2b - b23*r3b) + b31*(zAB + b31*r1b + b32*r2b - b33*r3b), b11*(xAB + b11*r1b + b12*r2b + b13*r3b) + b21*(yAB + b21*r1b + b22*r2b + b23*r3b) + b31*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ - b12*(a11*r1a + a12*r2a + a13*r3a) - b22*(a21*r1a + a22*r2a + a23*r3a) - b32*(a31*r1a + a32*r2a + a33*r3a), - b12*(a11*r1a + a12*r2a - a13*r3a) - b22*(a21*r1a + a22*r2a - a23*r3a) - b32*(a31*r1a + a32*r2a - a33*r3a), - b12*(a11*r1a - a12*r2a + a13*r3a) - b22*(a21*r1a - a22*r2a + a23*r3a) - b32*(a31*r1a - a32*r2a + a33*r3a), b12*(a12*r2a - a11*r1a + a13*r3a) + b22*(a22*r2a - a21*r1a + a23*r3a) + b32*(a32*r2a - a31*r1a + a33*r3a), - b12*(a12*r2a - a11*r1a + a13*r3a) - b22*(a22*r2a - a21*r1a + a23*r3a) - b32*(a32*r2a - a31*r1a + a33*r3a), b12*(a11*r1a - a12*r2a + a13*r3a) + b22*(a21*r1a - a22*r2a + a23*r3a) + b32*(a31*r1a - a32*r2a + a33*r3a), b12*(a11*r1a + a12*r2a - a13*r3a) + b22*(a21*r1a + a22*r2a - a23*r3a) + b32*(a31*r1a + a32*r2a - a33*r3a), b12*(a11*r1a + a12*r2a + a13*r3a) + b22*(a21*r1a + a22*r2a + a23*r3a) + b32*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - b12*(b11*r1b - xAB + b12*r2b + b13*r3b) - b22*(b21*r1b - yAB + b22*r2b + b23*r3b) - b32*(b31*r1b - zAB + b32*r2b + b33*r3b), b12*(xAB - b11*r1b - b12*r2b + b13*r3b) + b22*(yAB - b21*r1b - b22*r2b + b23*r3b) + b32*(zAB - b31*r1b - b32*r2b + b33*r3b), b12*(xAB - b11*r1b + b12*r2b - b13*r3b) + b22*(yAB - b21*r1b + b22*r2b - b23*r3b) + b32*(zAB - b31*r1b + b32*r2b - b33*r3b), b12*(xAB - b11*r1b + b12*r2b + b13*r3b) + b22*(yAB - b21*r1b + b22*r2b + b23*r3b) + b32*(zAB - b31*r1b + b32*r2b + b33*r3b), b12*(xAB + b11*r1b - b12*r2b - b13*r3b) + b22*(yAB + b21*r1b - b22*r2b - b23*r3b) + b32*(zAB + b31*r1b - b32*r2b - b33*r3b), b12*(xAB + b11*r1b - b12*r2b + b13*r3b) + b22*(yAB + b21*r1b - b22*r2b + b23*r3b) + b32*(zAB + b31*r1b - b32*r2b + b33*r3b), b12*(xAB + b11*r1b + b12*r2b - b13*r3b) + b22*(yAB + b21*r1b + b22*r2b - b23*r3b) + b32*(zAB + b31*r1b + b32*r2b - b33*r3b), b12*(xAB + b11*r1b + b12*r2b + b13*r3b) + b22*(yAB + b21*r1b + b22*r2b + b23*r3b) + b32*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ - b13*(a11*r1a + a12*r2a + a13*r3a) - b23*(a21*r1a + a22*r2a + a23*r3a) - b33*(a31*r1a + a32*r2a + a33*r3a), - b13*(a11*r1a + a12*r2a - a13*r3a) - b23*(a21*r1a + a22*r2a - a23*r3a) - b33*(a31*r1a + a32*r2a - a33*r3a), - b13*(a11*r1a - a12*r2a + a13*r3a) - b23*(a21*r1a - a22*r2a + a23*r3a) - b33*(a31*r1a - a32*r2a + a33*r3a), b13*(a12*r2a - a11*r1a + a13*r3a) + b23*(a22*r2a - a21*r1a + a23*r3a) + b33*(a32*r2a - a31*r1a + a33*r3a), - b13*(a12*r2a - a11*r1a + a13*r3a) - b23*(a22*r2a - a21*r1a + a23*r3a) - b33*(a32*r2a - a31*r1a + a33*r3a), b13*(a11*r1a - a12*r2a + a13*r3a) + b23*(a21*r1a - a22*r2a + a23*r3a) + b33*(a31*r1a - a32*r2a + a33*r3a), b13*(a11*r1a + a12*r2a - a13*r3a) + b23*(a21*r1a + a22*r2a - a23*r3a) + b33*(a31*r1a + a32*r2a - a33*r3a), b13*(a11*r1a + a12*r2a + a13*r3a) + b23*(a21*r1a + a22*r2a + a23*r3a) + b33*(a31*r1a + a32*r2a + a33*r3a)]
	cB = [ - b13*(b11*r1b - xAB + b12*r2b + b13*r3b) - b23*(b21*r1b - yAB + b22*r2b + b23*r3b) - b33*(b31*r1b - zAB + b32*r2b + b33*r3b), b13*(xAB - b11*r1b - b12*r2b + b13*r3b) + b23*(yAB - b21*r1b - b22*r2b + b23*r3b) + b33*(zAB - b31*r1b - b32*r2b + b33*r3b), b13*(xAB - b11*r1b + b12*r2b - b13*r3b) + b23*(yAB - b21*r1b + b22*r2b - b23*r3b) + b33*(zAB - b31*r1b + b32*r2b - b33*r3b), b13*(xAB - b11*r1b + b12*r2b + b13*r3b) + b23*(yAB - b21*r1b + b22*r2b + b23*r3b) + b33*(zAB - b31*r1b + b32*r2b + b33*r3b), b13*(xAB + b11*r1b - b12*r2b - b13*r3b) + b23*(yAB + b21*r1b - b22*r2b - b23*r3b) + b33*(zAB + b31*r1b - b32*r2b - b33*r3b), b13*(xAB + b11*r1b - b12*r2b + b13*r3b) + b23*(yAB + b21*r1b - b22*r2b + b23*r3b) + b33*(zAB + b31*r1b - b32*r2b + b33*r3b), b13*(xAB + b11*r1b + b12*r2b - b13*r3b) + b23*(yAB + b21*r1b + b22*r2b - b23*r3b) + b33*(zAB + b31*r1b + b32*r2b - b33*r3b), b13*(xAB + b11*r1b + b12*r2b + b13*r3b) + b23*(yAB + b21*r1b + b22*r2b + b23*r3b) + b33*(zAB + b31*r1b + b32*r2b + b33*r3b)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a11*b31 - a31*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a21*b31 - a31*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b21 - a21*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a11*b31 - a31*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a21*b31 - a31*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b21 - a21*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a11*b31 - a31*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a21*b31 - a31*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b21 - a21*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a11*b31 - a31*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a21*b31 - a31*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b21 - a21*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2)]
	cB = [ ((a11*b31 - a31*b11)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a21*b31 - a31*b21)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b21 - a21*b11)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2), ((a21*b31 - a31*b21)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) - ((a11*b31 - a31*b11)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2) + ((a11*b21 - a21*b11)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a11*b21 - a21*b11)^2 + (a11*b31 - a31*b11)^2 + (a21*b31 - a31*b21)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a11*b32 - a31*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a21*b32 - a31*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b22 - a21*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a11*b32 - a31*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a21*b32 - a31*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b22 - a21*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a11*b32 - a31*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a21*b32 - a31*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b22 - a21*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a11*b32 - a31*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a21*b32 - a31*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b22 - a21*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2)]
	cB = [ ((a11*b32 - a31*b12)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a21*b32 - a31*b22)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b22 - a21*b12)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2), ((a21*b32 - a31*b22)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) - ((a11*b32 - a31*b12)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2) + ((a11*b22 - a21*b12)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a11*b22 - a21*b12)^2 + (a11*b32 - a31*b12)^2 + (a21*b32 - a31*b22)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a11*b33 - a31*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a21*b33 - a31*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b23 - a21*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a11*b33 - a31*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a21*b33 - a31*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b23 - a21*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a11*b33 - a31*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a21*b33 - a31*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b23 - a21*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a11*b33 - a31*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a21*b33 - a31*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b23 - a21*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2)]
	cB = [ ((a11*b33 - a31*b13)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a21*b33 - a31*b23)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b23 - a21*b13)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2), ((a21*b33 - a31*b23)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) - ((a11*b33 - a31*b13)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2) + ((a11*b23 - a21*b13)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a11*b23 - a21*b13)^2 + (a11*b33 - a31*b13)^2 + (a21*b33 - a31*b23)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a12*b31 - a32*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a22*b31 - a32*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b21 - a22*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a12*b31 - a32*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a22*b31 - a32*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b21 - a22*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a12*b31 - a32*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a22*b31 - a32*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b21 - a22*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a12*b31 - a32*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a22*b31 - a32*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b21 - a22*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2)]
	cB = [ ((a12*b31 - a32*b11)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a22*b31 - a32*b21)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b21 - a22*b11)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2), ((a22*b31 - a32*b21)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) - ((a12*b31 - a32*b11)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2) + ((a12*b21 - a22*b11)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a12*b21 - a22*b11)^2 + (a12*b31 - a32*b11)^2 + (a22*b31 - a32*b21)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a12*b32 - a32*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a22*b32 - a32*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b22 - a22*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a12*b32 - a32*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a22*b32 - a32*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b22 - a22*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a12*b32 - a32*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a22*b32 - a32*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b22 - a22*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a12*b32 - a32*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a22*b32 - a32*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b22 - a22*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2)]
	cB = [ ((a12*b32 - a32*b12)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a22*b32 - a32*b22)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b22 - a22*b12)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2), ((a22*b32 - a32*b22)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) - ((a12*b32 - a32*b12)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2) + ((a12*b22 - a22*b12)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a12*b22 - a22*b12)^2 + (a12*b32 - a32*b12)^2 + (a22*b32 - a32*b22)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a12*b33 - a32*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a22*b33 - a32*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b23 - a22*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a12*b33 - a32*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a22*b33 - a32*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b23 - a22*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a12*b33 - a32*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a22*b33 - a32*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b23 - a22*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a12*b33 - a32*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a22*b33 - a32*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b23 - a22*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2)]
	cB = [ ((a12*b33 - a32*b13)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a22*b33 - a32*b23)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b23 - a22*b13)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2), ((a22*b33 - a32*b23)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) - ((a12*b33 - a32*b13)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2) + ((a12*b23 - a22*b13)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a12*b23 - a22*b13)^2 + (a12*b33 - a32*b13)^2 + (a22*b33 - a32*b23)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a13*b31 - a33*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a23*b31 - a33*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b21 - a23*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a13*b31 - a33*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a23*b31 - a33*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b21 - a23*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a13*b31 - a33*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a23*b31 - a33*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b21 - a23*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a13*b31 - a33*b11)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a23*b31 - a33*b21)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b21 - a23*b11)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2)]
	cB = [ ((a13*b31 - a33*b11)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a23*b31 - a33*b21)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b21 - a23*b11)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2), ((a23*b31 - a33*b21)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) - ((a13*b31 - a33*b11)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2) + ((a13*b21 - a23*b11)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a13*b21 - a23*b11)^2 + (a13*b31 - a33*b11)^2 + (a23*b31 - a33*b21)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a13*b32 - a33*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a23*b32 - a33*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b22 - a23*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a13*b32 - a33*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a23*b32 - a33*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b22 - a23*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a13*b32 - a33*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a23*b32 - a33*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b22 - a23*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a13*b32 - a33*b12)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a23*b32 - a33*b22)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b22 - a23*b12)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2)]
	cB = [ ((a13*b32 - a33*b12)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a23*b32 - a33*b22)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b22 - a23*b12)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2), ((a23*b32 - a33*b22)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) - ((a13*b32 - a33*b12)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2) + ((a13*b22 - a23*b12)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a13*b22 - a23*b12)^2 + (a13*b32 - a33*b12)^2 + (a23*b32 - a33*b22)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	cA = [ ((a13*b33 - a33*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a23*b33 - a33*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b23 - a23*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a13*b33 - a33*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a23*b33 - a33*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b23 - a23*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a13*b33 - a33*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a23*b33 - a33*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b23 - a23*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a13*b33 - a33*b13)*(a22*r2a - a21*r1a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a23*b33 - a33*b23)*(a12*r2a - a11*r1a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b23 - a23*b13)*(a32*r2a - a31*r1a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(a11*r1a - a12*r2a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(a21*r1a - a22*r2a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(a31*r1a - a32*r2a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(a11*r1a + a12*r2a - a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(a21*r1a + a22*r2a - a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(a31*r1a + a32*r2a - a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(a11*r1a + a12*r2a + a13*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(a21*r1a + a22*r2a + a23*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(a31*r1a + a32*r2a + a33*r3a))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2)]
	cB = [ ((a13*b33 - a33*b13)*(b21*r1b - yAB + b22*r2b + b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a23*b33 - a33*b23)*(b11*r1b - xAB + b12*r2b + b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b23 - a23*b13)*(b31*r1b - zAB + b32*r2b + b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB - b11*r1b - b12*r2b + b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB - b21*r1b - b22*r2b + b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB - b31*r1b - b32*r2b + b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB - b11*r1b + b12*r2b - b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB - b21*r1b + b22*r2b - b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB - b31*r1b + b32*r2b - b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB - b11*r1b + b12*r2b + b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB - b21*r1b + b22*r2b + b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB - b31*r1b + b32*r2b + b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB + b11*r1b - b12*r2b - b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB + b21*r1b - b22*r2b - b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB + b31*r1b - b32*r2b - b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB + b11*r1b - b12*r2b + b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB + b21*r1b - b22*r2b + b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB + b31*r1b - b32*r2b + b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB + b11*r1b + b12*r2b - b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB + b21*r1b + b22*r2b - b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB + b31*r1b + b32*r2b - b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2), ((a23*b33 - a33*b23)*(xAB + b11*r1b + b12*r2b + b13*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) - ((a13*b33 - a33*b13)*(yAB + b21*r1b + b22*r2b + b23*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2) + ((a13*b23 - a23*b13)*(zAB + b31*r1b + b32*r2b + b33*r3b))/((a13*b23 - a23*b13)^2 + (a13*b33 - a33*b13)^2 + (a23*b33 - a33*b23)^2)^(1/2)]
	if minimum(cA) >= maximum(cB) || maximum(cA) <= minimum(cB)
		overlapfun = 0.0
		return overlapfun
	end

	return overlapfun
end


text:
using Thebes, Luxor, Colors

@png begin
    setlinejoin("bevel")
    helloworld()
    perspective(300)
    eyepoint(800, 700, 600)
    uppoint(0, 0, 10)
    setopacity(0.7)
    k = 50
    xmax = 6π
    ymax = 6π
    st = 0.15
    f(x, y) = 2(sin(x) * cos(y)) + (cos(x) * sin(y))
    setline(.5)
    fontsize(5)
    for x in -xmax:st:xmax
        for y in -ymax:st:ymax
            sethue(HSB(360rescale(f(x, y), -2, 2), .7, .7))
            p1 = Point3D(k * x,        k * y,       k * f(x, y))
            p2 = Point3D(k * (x ),     k * (y + st), k * f(x, y +st))
            p3 = Point3D(k * (x + st),  k * (y + st), k * f(x + st, y + st))
            p4 = Point3D(k * (x + st),  k * y,       k * f(x + st, y))
            pin([p1, p2, p3, p4], gfunction = (p3s, p2s) -> begin
                poly(p2s, close=true, :fill)
                sethue("red")
                poly(p2s, close=true, :stroke)
            end)
        end
    end
    axes3D()
end 500 500 "/tmp/t.png"

ext: jl


text:
""" Generate and parse Julia programs based on diagrams.
"""
module Programs

using Reexport

include("GenerateJuliaPrograms.jl")
include("ParseJuliaPrograms.jl")
include("RelationalPrograms.jl")
include("DiagrammaticPrograms.jl")

@reexport using .GenerateJuliaPrograms
@reexport using .ParseJuliaPrograms
@reexport using .RelationalPrograms
@reexport using .DiagrammaticPrograms

end


text:
#src this layout is heavily inspired by how examples in https://github.com/Ferrite-FEM/Ferrite.jl are generated

# # Mutagenesis Example
# Following example demonstrates learning to [predict the mutagenicity on Salmonella typhimurium](https://relational.fit.cvut.cz/dataset/Mutagenesis) (dataset is stored in json format [in MLDatasets.jl](https://juliaml.github.io/MLDatasets.jl/stable/datasets/Mutagenesis/) for your convenience).

#md # !!! tip
#md #     This example is also available as a Jupyter notebook, feel free to run it yourself:
#md #     [`mutagenesis.ipynb`](@__NBVIEWER_ROOT_URL__/examples/mutagenesis.ipynb)

#nb # We start by installing JsonGrinder and few other packages we need for the example.
#nb # Julia Ecosystem follows philosophy of many small single-purpose composable packages
#nb # which may be different from e.g. python where we usually use fewer larger packages.
#nb using Pkg
#nb pkg"add JsonGrinder#master MLDatasets Flux Mill#master MLDataPattern Statistics ChainRulesCore"

# Here we include libraries all necessary libraries
using JsonGrinder, MLDatasets, Flux, Mill, MLDataPattern, Statistics, ChainRulesCore

# Here we load all samples.
train_x, train_y = MLDatasets.Mutagenesis.traindata();
test_x, test_y = MLDatasets.Mutagenesis.testdata();

# We define some basic parameters for the construction and training of the neural network.
# Minibatch size is self-explanatory, iterations is number of iterations of gradient descent
# Neurons is number of neurons in hidden layers for each version of part of the neural network.
minibatchsize = 100
iterations = 5_000
neurons = 20

# We create the schema of the training data, which is the first important step in using the JsonGrinder.
# This computes both the structure (also known as JSON schema) and histogram of occurrences of individual values in the training data.
sch = JsonGrinder.schema(train_x)

# Then we use it to create the extractor converting jsons to Mill structures.
# The `suggestextractor` is executed below with default setting, but it allows you heavy customization.
extractor = suggestextractor(sch)

# We convert jsons to mill data samples and prepare list of classes. This classification problem is two-class, but we want to infer it from labels.
# The extractor is callable, so we can pass it vector of samples to obtain vector of structures with extracted features.
train_data = extractor.(train_x)
test_data = extractor.(test_x)
labelnames = unique(train_y)

# # Create the model
# We create the model reflecting structure of the data
model = reflectinmodel(sch, extractor,
	layer -> Dense(layer, neurons, relu),
	bag -> SegmentedMeanMax(bag),
	fsm = Dict("" => layer -> Dense(layer, length(labelnames))),
)
# this allows us to create model flexibly, without the need to hardcode individual layers.
# Individual arguments of `reflectinmodel` are explained in [Mill.jl documentation](https://CTUAvastLab.github.io/Mill.jl/dev/manual/reflectin/#Model-Reflection). But briefly: for every numeric array in the sample, model will create a dense layer with `neurons` neurons (20 in this example). For every vector of observations (called bag in Multiple Instance Learning terminology), it will create aggregation function which will take mean, maximum of feature vectors and concatenate them. The `fsm` keyword argument basically says that on the end of the NN, as a last layer, we want 2 neurons `length(labelnames)` in the output layer, not 20 as in the intermediate layers.

# # Train the model
# Then, we define few handy functions and a loss function, which is categorical crossentropy in our case.
loss(x,y) = Flux.logitcrossentropy(inference(x), Flux.onehotbatch(y, labelnames))
inference(x::AbstractMillNode) = model(x).data
inference(x::AbstractVector{<:AbstractMillNode}) = inference(reduce(catobs, x))
accuracy(x,y) = mean(labelnames[Flux.onecold(inference(x))] .== y)
loss(xy::Tuple) = loss(xy...)
@non_differentiable Base.reduce(catobs, x::AbstractVector{<:AbstractMillNode})

# And we can add a callback which will be printing train and test accuracy during the training
# and then we can start trining
cb = () -> begin
	train_acc = accuracy(train_data, train_y)
	test_acc = accuracy(test_data, test_y)
	println("accuracy: train = $train_acc, test = $test_acc")
end

# Lastly we turn our training data to minibatches, and we can start training
minibatches = RandomBatches((train_data, train_y), size = minibatchsize, count = iterations)
Flux.Optimise.train!(loss, Flux.params(model), minibatches, ADAM(), cb = Flux.throttle(cb, 2))
# We can see the accuracy rising and obtaining over 98% on training set quite quickly, and on test set we get over 70%.

# # Classify test set
# The Last part is inference on test data.
probs = softmax(inference(test_data))
o = Flux.onecold(probs)
pred_classes = labelnames[o]
mean(pred_classes .== test_y)

# `pred_classes` contains the predictions for our test set.
# we see the accuracy is around 75% on test set
# predicted classes for test set
pred_classes
# Ground truth classes for test set
test_y
# probabilities for test set
probs

# We can look at individual samples. For instance, some sample from test set is
test_data[2]

# and the corresponding classification is
pred_classes[2]

# if you want to see the probability distribution, it can be obtained by applying `softmax` to the output of the network.
softmax(model(test_data[2]).data)

# so we can see that the probability that given sample is `mutagenetic` is almost 1.


text:
T = 2
netload_min = 0.0 * ones(T)
netload_max = 10.0 * ones(T)
p_mkt = [2.0, 3.0]
p_tou = [1.0, 3.0]

# build houses
b1 = DR.DER.Battery(
    index=1,
    T=T,
    dt=1.0,
    soc_min=zeros(T), soc_max=10.0*ones(T), soc_init=0.0, selfdischargerate=1.0,
    charge_pwr_min=0.0, charge_pwr_max=6.0, discharge_pwr_min=0.0, discharge_pwr_max=6.0,
    charge_eff=1.0, discharge_eff=1.0
)
b2 = DR.DER.Battery(
    index=2,
    T=T,
    dt=1.0,
    soc_min=zeros(T), soc_max=10.0*ones(T), soc_init=0.0, selfdischargerate=1.0,
    charge_pwr_min=0.0, charge_pwr_max=6.0, discharge_pwr_min=0.0, discharge_pwr_max=6.0,
    charge_eff=1.0, discharge_eff=1.0
)
l1 = DR.DER.FixedLoad(index=1, T=T, load=ones(T))
l2 = DR.DER.FixedLoad(index=2, T=T, load=2.0*ones(T))
h1 = DR.DER.House(index=1, T=T, netload_min=zeros(T),netload_max=5.0*ones(T), price=0.0*p_tou, appliances=[b1, l1])
h2 = DR.DER.House(index=2, T=T, netload_min=zeros(T),netload_max=5.0*ones(T), price=0.0*p_tou, appliances=[b2, l2])

# Build aggregator
agg = DR.Aggregator(T, netload_min, netload_max, p_mkt, [h1, h2])

# Build centralized model
m, var2idx, row2idx = DR.buildmodel(agg, GLPKSolverMIP())
MPB.optimize!(m)

@printf("Obj: %f\n", MPB.getobjval(m))

# Solve DW relaxation


text:
#!/usr/bin/env julia

# This file tests generation of multiblock VTK files.
# It also tests the non-default values of the optional parameters of vtk_grid,
# i.e., compress=false and append=false.

using WriteVTK
using StaticArrays: SVector
import OffsetArrays: OffsetArray

const FloatType = Float32
const vtm_filename_noext = "multiblock"

function first_block_data()
    Ni, Nj, Nk = 10, 15, 20

    x = zeros(FloatType, Ni, Nj, Nk)
    y = copy(x)
    z = copy(x)

    # Just to test subarrays:
    qall = zeros(FloatType, 3Ni, 2Nj, 2Nk)
    q = view(qall, 1:Ni, 1:Nj, 1:Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r::FloatType = 1 + (i - 1)/(Ni - 1)
        th::FloatType = 3*pi/4 * (j-1) / (Nj-1)
        x[i, j, k] = r * cos(th)
        y[i, j, k] = r * sin(th)
        z[i, j, k] = (k-1) / (Nk-1)
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function second_block_data()
    Ni, Nj, Nk = 20, 16, 12

    x = zeros(FloatType, Ni)
    y = zeros(FloatType, Nj)
    z = zeros(FloatType, Nk)

    [x[i] =  1 +   (i-1) / (Ni-1) for i = 1:Ni]
    [y[j] = -2 + 2*(j-1) / (Nj-1) for j = 1:Nj]
    [z[k] =        (k-1) / (Nk-1) for k = 1:Nk]

    q = zeros(FloatType, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function third_block_data()
    # This is basically the same as in the unstructured example.
    Ni, Nj, Nk = 40, 50, 20
    indices = LinearIndices((1:Ni, 1:Nj, 1:Nk))

    # Create points and point data.
    pts = Array{FloatType}(undef, 3, Ni, Nj, Nk)
    pdata = Array{FloatType}(undef, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r = 1 + (i - 1)/(Ni - 1)
        th = 3*pi/4 * (j - 1)/(Nj - 1)
        pts[1, i, j, k] = r * cos(th)
        pts[2, i, j, k] = r * sin(th)
        pts[3, i, j, k] = -2 + (k - 1)/(Nk - 1)
        pdata[i, j, k] = i*i + k*sqrt(j)
    end

    # Create cells (all hexahedrons in this case) and cell data.
    celltype = VTKCellTypes.VTK_HEXAHEDRON
    cells = MeshCell[]
    cdata = FloatType[]

    for k = 2:Nk, j = 2:Nj, i = 2:Ni
        # Define connectivity of cell.
        inds = SVector{8, Int32}(
            indices[i-1, j-1, k-1],
            indices[i  , j-1, k-1],
            indices[i  , j  , k-1],
            indices[i-1, j  , k-1],
            indices[i-1, j-1, k  ],
            indices[i  , j-1, k  ],
            indices[i  , j  , k  ],
            indices[i-1, j  , k  ])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i*j*k)
    end

    return pts, cells, pdata, cdata
end

function fourth_block_data()
    # Create unstructured example with boundaries.
    imax, jmax, kmax = 40, 50, 20
    indices = LinearIndices((1:imax, 1:jmax, 1:kmax))

    # Create points and point data.
    pts_vol = Array{FloatType}(undef, 3, imax, jmax, kmax)
    pdata_vol = Array{FloatType}(undef, imax, jmax, kmax)

    x = pts_vol[1, :, :, :] .= reshape(range(0.0, stop=1.0, length=imax), imax, 1, 1)
    y = pts_vol[2, :, :, :] .= reshape(range(-1.0, stop=0.0, length=jmax), 1, jmax, 1)
    z = pts_vol[3, :, :, :] .= reshape(range(3.0, stop=4.0, length=kmax), 1, 1, kmax)

    @. pdata_vol = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    # Create cells (all tetrahedra in this case) and cell data.
    celltype = VTKCellTypes.VTK_TETRA
    cells_vol = MeshCell[]
    cdata_vol = FloatType[]

    for k = 1:kmax-1, j = 1:jmax-1, i = 1:imax-1
        # 1, 3, 4, 5
        inds = SVector{4, Int32}(
            indices[i, j, k],
            indices[i+1, j, k+1],
            indices[i, j, k+1],
            indices[i, j+1, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

        # 8, 5, 4, 3
        inds = SVector{4, Int32}(
            indices[i, j+1, k+1],
            indices[i, j+1, k],
            indices[i, j, k+1],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

        # 8, 7, 5, 3
        inds = SVector{4, Int32}(
            indices[i, j+1, k+1],
            indices[i+1, j+1, k+1],
            indices[i, j+1, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

        # 7, 6, 5, 3
        inds = SVector{4, Int32}(
            indices[i+1, j+1, k+1],
            indices[i+1, j+1, k],
            indices[i, j+1, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

        # 6, 2, 5, 3
        inds = SVector{4, Int32}(
            indices[i+1, j+1, k],
            indices[i+1, j, k],
            indices[i, j+1, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

        # 5, 2, 1, 3
        inds = SVector{4, Int32}(
            indices[i, j+1, k],
            indices[i+1, j, k],
            indices[i, j, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_vol, c)
        push!(cdata_vol, i*j*k)

    end

    # Create the block for the imin patch.

    indices = LinearIndices((1:1, 1:jmax, 1:kmax))

    # Create points and point data for the imin patch.
    pts_imin = Array{FloatType}(undef, 3, 1, jmax, kmax)
    pdata_imin = Array{FloatType}(undef, 1, jmax, kmax)

    x = pts_imin[1, :, :, :] .= pts_vol[1, 1:1, :, :]
    y = pts_imin[2, :, :, :] .= pts_vol[2, 1:1, :, :]
    z = pts_imin[3, :, :, :] .= pts_vol[3, 1:1, :, :]

    @. pdata_imin = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_imin = MeshCell[]
    cdata_imin = FloatType[]

    for k = 1:kmax-1, j = 1:jmax-1, i = 1:1
        # 1, 3, 4, 5
        inds = SVector{3, Int32}(
            indices[i, j, k],
            indices[i, j, k+1],
            indices[i, j+1, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_imin, c)
        push!(cdata_imin, i*j*k)

        # 8, 5, 4, 3
        inds = SVector{3, Int32}(
            indices[i, j+1, k+1],
            indices[i, j+1, k],
            indices[i, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_imin, c)
        push!(cdata_imin, i*j*k)

    end

    # Create the block for the imax patch.

    indices = LinearIndices((imax:imax, 1:jmax, 1:kmax))
    # Use an OffsetArray to make the indexing easier.
    indices = OffsetArray(indices, (imax:imax, 1:jmax, 1:kmax))

    # Create points and point data for the imin patch.
    pts_imax = Array{FloatType}(undef, 3, 1, jmax, kmax)
    pdata_imax = Array{FloatType}(undef, 1, jmax, kmax)

    x = pts_imax[1, :, :, :] .= pts_vol[1, imax:imax, :, :]
    y = pts_imax[2, :, :, :] .= pts_vol[2, imax:imax, :, :]
    z = pts_imax[3, :, :, :] .= pts_vol[3, imax:imax, :, :]

    @. pdata_imax = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_imax = MeshCell[]
    cdata_imax = FloatType[]

    for k = 1:kmax-1, j = 1:jmax-1, i = imax-1:imax-1
        # 7, 6, 5, 3
        inds = SVector{3, Int32}(
            indices[i+1, j+1, k+1],
            indices[i+1, j+1, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_imax, c)
        push!(cdata_imax, i*j*k)

        # 6, 2, 5, 3
        inds = SVector{3, Int32}(
            indices[i+1, j+1, k],
            indices[i+1, j, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_imax, c)
        push!(cdata_imax, i*j*k)

    end

    # Create the block for the jmin patch.

    indices = LinearIndices((1:imax, 1:1, 1:kmax))

    # Create points and point data for the jmin patch.
    pts_jmin = Array{FloatType}(undef, 3, imax, 1, kmax)
    pdata_jmin = Array{FloatType}(undef, imax, 1, kmax)

    x = pts_jmin[1, :, :, :] .= pts_vol[1, 1:imax, 1:1, :]
    y = pts_jmin[2, :, :, :] .= pts_vol[2, 1:imax, 1:1, :]
    z = pts_jmin[3, :, :, :] .= pts_vol[3, 1:imax, 1:1, :]

    @. pdata_jmin = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_jmin = MeshCell[]
    cdata_jmin = FloatType[]

    for k = 1:kmax-1, j = 1:1, i = 1:imax-1
        # 1, 3, 4, 5
        inds = SVector{3, Int32}(
            indices[i, j, k],
            indices[i+1, j, k+1],
            indices[i, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_jmin, c)
        push!(cdata_jmin, i*j*k)

        # 5, 2, 1, 3
        inds = SVector{3, Int32}(
            indices[i+1, j, k],
            indices[i, j, k],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_jmin, c)
        push!(cdata_jmin, i*j*k)

    end

    # Create the block for the jmax patch.

    indices = LinearIndices((1:imax, jmax:jmax, 1:kmax))
    # Use an OffsetArray to make the indexing easier.
    indices = OffsetArray(indices, (1:imax, jmax:jmax, 1:kmax))

    # Create points and point data for the imin patch.
    pts_jmax = Array{FloatType}(undef, 3, imax, 1, kmax)
    pdata_jmax = Array{FloatType}(undef, imax, 1, kmax)

    x = pts_jmax[1, :, :, :] .= pts_vol[1, :, jmax:jmax, :]
    y = pts_jmax[2, :, :, :] .= pts_vol[2, :, jmax:jmax, :]
    z = pts_jmax[3, :, :, :] .= pts_vol[3, :, jmax:jmax, :]

    @. pdata_jmax = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_jmax = MeshCell[]
    cdata_jmax = FloatType[]

    for k = 1:kmax-1, j = jmax-1:jmax-1, i = 1:imax-1
        # 8, 7, 5, 3
        inds = SVector{3, Int32}(
            indices[i, j+1, k+1],
            indices[i+1, j+1, k+1],
            indices[i, j+1, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_jmax, c)
        push!(cdata_jmax, i*j*k)

        # 7, 6, 5, 3
        inds = SVector{3, Int32}(
            indices[i+1, j+1, k+1],
            indices[i+1, j+1, k],
            indices[i, j+1, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_jmax, c)
        push!(cdata_jmax, i*j*k)

    end

    # Create the block for the kmin patch.

    indices = LinearIndices((1:imax, 1:jmax, 1:1))

    # Create points and point data for the kmin patch.
    pts_kmin = Array{FloatType}(undef, 3, imax, jmax, 1)
    pdata_kmin = Array{FloatType}(undef, imax, jmax, 1)

    x = pts_kmin[1, :, :, :] .= pts_vol[1, :, :, 1:1]
    y = pts_kmin[2, :, :, :] .= pts_vol[2, :, :, 1:1]
    z = pts_kmin[3, :, :, :] .= pts_vol[3, :, :, 1:1]

    @. pdata_kmin = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_kmin = MeshCell[]
    cdata_kmin = FloatType[]

    for k = 1:1, j = 1:jmax-1, i = 1:imax-1
        # 6, 2, 5, 3
        inds = SVector{3, Int32}(
            indices[i+1, j+1, k],
            indices[i+1, j, k],
            indices[i, j+1, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_kmin, c)
        push!(cdata_kmin, i*j*k)

        # 5, 2, 1, 3
        inds = SVector{3, Int32}(
            indices[i, j+1, k],
            indices[i+1, j, k],
            indices[i, j, k])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_kmin, c)
        push!(cdata_kmin, i*j*k)

    end

    # Create the block for the kmax patch.

    indices = LinearIndices((1:imax, 1:jmax, kmax:kmax))
    # Use an OffsetArray to make the indexing easier.
    indices = OffsetArray(indices, (1:imax, 1:jmax, kmax:kmax))

    # Create points and point data for the imin patch.
    pts_kmax = Array{FloatType}(undef, 3, imax, jmax, 1)
    pdata_kmax = Array{FloatType}(undef, imax, jmax, 1)

    x = pts_kmax[1, :, :, :] .= pts_vol[1, :, :, kmax:kmax]
    y = pts_kmax[2, :, :, :] .= pts_vol[2, :, :, kmax:kmax]
    z = pts_kmax[3, :, :, :] .= pts_vol[3, :, :, kmax:kmax]

    @. pdata_kmax = sin(2*pi*x)*sin(4*pi*y)*sin(6*pi*z)

    celltype = VTKCellTypes.VTK_TRIANGLE
    cells_kmax = MeshCell[]
    cdata_kmax = FloatType[]

    for k = kmax-1:kmax-1, j = 1:jmax-1, i = 1:imax-1
        # 8, 5, 4, 3
        inds = SVector{3, Int32}(
            indices[i, j+1, k+1],
            indices[i, j, k+1],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_kmax, c)
        push!(cdata_kmax, i*j*k)

        # 8, 7, 5, 3
        inds = SVector{3, Int32}(
            indices[i, j+1, k+1],
            indices[i+1, j+1, k+1],
            indices[i+1, j, k+1])

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells_kmax, c)
        push!(cdata_kmax, i*j*k)

    end

    return ((pts_vol, cells_vol, pdata_vol, cdata_vol),
            (pts_imin, cells_imin, pdata_imin, cdata_imin),
            (pts_imax, cells_imax, pdata_imax, cdata_imax),
            (pts_jmin, cells_jmin, pdata_jmin, cdata_jmin),
            (pts_jmax, cells_jmax, pdata_jmax, cdata_jmax),
            (pts_kmin, cells_kmin, pdata_kmin, cdata_kmin),
            (pts_kmax, cells_kmax, pdata_kmax, cdata_kmax))

end

function main()
    # Generate data.
    x1, y1, z1, q1 = first_block_data()
    x2, y2, z2, q2 = second_block_data()
    points3, cells3, q3, c3 = third_block_data()
    ((points4, cells4, q4, c4),
     (points41, cells41, q41, c41),
     (points42, cells42, q42, c42),
     (points43, cells43, q43, c43),
     (points44, cells44, q44, c44),
     (points45, cells45, q45, c45),
     (points46, cells46, q46, c46)) = fourth_block_data()

    # Create multiblock file.
    @time outfiles = vtk_multiblock(vtm_filename_noext) do vtm
        # First block.
        vtk = vtk_grid(vtm, x1, y1, z1; compress=false, append=false)
        vtk["q_values"] = q1

        # Second block.
        vtk = vtk_grid(vtm, x2, y2, z2; compress=false)
        vtk["q_values"] = q2

        # Third block.
        vtk = vtk_grid(vtm, points3, cells3; append=false)
        vtk["q_values"] = q3
        vtk["c_values"] = c3

        # Fourth block is a collection of multiblock files itself.
        block = multiblock_add_block(vtm, "multiblock_4")
        vtk = vtk_grid(block, "multiblock_4_volume", points4, cells4)
        vtk["q_values"] = q4
        vtk["c_values"] = c4
        vtk = vtk_grid(block, "multiblock_4_imin", points41, cells41)
        vtk["q_values"] = q41
        vtk["c_values"] = c41
        vtk = vtk_grid(block, "multiblock_4_imax", points42, cells42)
        vtk["q_values"] = q42
        vtk["c_values"] = c42
        vtk = vtk_grid(block, "multiblock_4_jmin", points43, cells43)
        vtk["q_values"] = q43
        vtk["c_values"] = c43
        vtk = vtk_grid(block, "multiblock_4_jmax", points44, cells44)
        vtk["q_values"] = q44
        vtk["c_values"] = c44
        vtk = vtk_grid(block, "multiblock_4_kmin", points45, cells45)
        vtk["q_values"] = q45
        vtk["c_values"] = c45
        vtk = vtk_grid(block, "multiblock_4_kmax", points46, cells46)
        vtk["q_values"] = q46
        vtk["c_values"] = c46
    end
    println("Saved:  ", join(outfiles, "  "))

    return outfiles::Vector{String}
end

main()