Building on from my previous question, how can I plot something more complex than a sphere (and not as easy to split into a set of parametric equations)?
Say, for example
$$2x^2+y^2-5z^2+z-7x=16$$
Graphed on Desmos (use the z-slider to change the z value) for those who cannot as easily visualize the 3D surface.
As far as I know, there isn't an add-on or method that can graph more complex mathematical functions in the way you want. As an alternative, you can put your equation into nodes and then render it with a Volume Scatter shader.
Begin with the default cube. Add a new material, and add a Texture Coordinate node (Shift + A > Input > Texture Coordinate) and Separate XYZ node (Shift + A > Converter > Separate XYZ):
Then, create your equation using Math nodes (Shift + A > Converter > Math). Here, I created your example equation:
Now, add a Color Mix node (Shift + A > Color > MixRGB) and set it to Multiply. Plug the texture coordinates in to one slot, and a Value node in to the other. This will scale our window:
Finally, multiply the output by 10 to increase the value, then plug it into the Density of a Volume Scatter and a Volume Absorption shader added together:
This will nicely display the function when rendered, though the shading isn't ideal:
click for larger image
To optimize rendering, change the following settings:
Although it doesn't give you a workable mesh, this can definitely give you a renderable graph.
We can also improve the look of it by coloring it:
click for larger image
A possible approch that avoids solving the equation would consist in interating through all possible points in the 3D space, check whether or not in each one of them the equation is verified, and if it is create a vertex there.
As in the 3D space there are infinite points, we must discretize it and limits the boundaries of the investigation domain. For the sake of semplicity, I choose to consider, for example, only integer numbers from -100 to +100 for each axis. Here's the code for such iteration:
for x in range(-100,100):
for y in range(-100,100):
for z in range(-100,100):
Then, for each vector we'll calculate the left member of the equation and check if it's equal to the right member of the equation previusly defined. As we discretized the domain, it's best to give to the result a bit of "tollerance" in verifying the equation (that's why I'm using the "+/-toll" trick upon the right member verification).
function=2*x**2+y**2-5*z**2+z-7*x #function here
if((function<(goal+toll))and(function>(goal-toll))):
bm.verts.new([x,y,z])
Here's a possible code with all the lines needed to create an object in the scene:
import bpy
import bmesh
mesh = bpy.data.meshes.new("mesh")
obj = bpy.data.objects.new("Function", mesh)
scene = bpy.context.scene
scene.objects.link(obj)
scene.objects.active = obj
obj.select = True
mesh = bpy.context.object.data
bm = bmesh.new()
goal=16 #right member of the equation
toll=5 #threshold
for x in range(-100,100):
for y in range(-100,100):
for z in range(-100,100):
function=2*x**2+y**2-5*z**2+z-7*x #function here
if((function<(goal+toll))and(function>(goal-toll))):
bm.verts.new([x,y,z])
bm.to_mesh(mesh)
bm.free()
Examples of usage:
function=2*x**2+y**2-5*z**2+z-7*x
Note: as you can see, it takes a while to iterate through all the 8 000 000 positions...
This script only provides some of the vertices due to the approximation, but, with a proper resolution, you can obtain a good representation of the surface.
The meshing of the vertices into a surface needs it's own approach. This "brute force plotting" I'm proposing does not go well with the "on the fly" definitions of faces. The set of point is not structured in any way. You'll need something to mesh the pointcloud.
Here's another version of the code where you can better control the domain size and resolution using floats to inspect the region closer to the origin in high definition:
import bpy
import bmesh
max=7
min=-7
prec=10 #precision of the calc (1=1 units, 10 = 0.1 units, 100 = 0.01 units...)
tol=0.5 #tollerance
def f(x, y, z): return 2*x**2+y**2-5*z**2+z-7*x -16 #the function = 0
mesh = bpy.data.meshes.new("mesh")
obj = bpy.data.objects.new("XYZ Function", mesh)
scene = bpy.context.scene
scene.objects.link(obj)
scene.objects.active = obj
obj.select = True
mesh = bpy.context.object.data
bm = bmesh.new()
for x in [float(j)/prec for j in range(min*prec, max*prec+1,1)]:
for y in [float(j)/prec for j in range(min*prec, max*prec+1,1)]:
for z in [float(j)/prec for j in range(min*prec, max*prec+1,1)]:
if(abs(f(x,y,z))<tol):
bm.verts.new([x,y,z])
bm.to_mesh(mesh)
bm.free()
Here's a possible result:
You can for example export points to Meshlab and take advantage of the built in tools. Here's an example using the previusly obtained set of points.
The equation in the form in the question is a quadric surface. By completing the square we can translate the axis and reduce to form
(x-7/4)² - 49/16 + 1/2y² - 5/2(z - 1/10)² + 5/200 - 8 = 0
or
X² + 0.5Y² - 2.5Z² - - 11.0375 = 0
A sample of the surface created by iterating through z (an xy plane) then finding the 0, 1 or 2 roots for each y using the quadratic formula to solve for x. Notice the missing faces where x=0. Can run thru solving for y and have missing faces around y=0. Overlaying the two gives the shape, but not a pretty mesh.
I've put in both the original equation f(x, y, z) and the translated form f2(x, y, z).
import bpy
import bmesh
from math import sqrt
def f(x, y, z):
return 2 * x * x + y * y - 5 * z * z + z -7 * x -16
# f transformed to have origin (0, 0, 0)
def f2(x, y, z):
x * x + 0.5 * y * y - 2.5 * z * z - 11.0375
def quad_root(a, b, c):
axis = -b / (2 * a)
check = b * b - 4 * a * c
if abs(check) < 0.0001:
return True, axis, axis
if check < 0:
return False, None, None
dist = sqrt(check) / (2 * a)
return True, axis - dist, axis + dist
bm = bmesh.new()
dom = [m for m in range(-10, 11)]
use_f2 = False # use f2()
planes = []
# coefficents of x
if use_f2:
a, b = 1, 0 # f2()
else:
a, b = 2, -7 # f()
for z in dom:
# for each xy plane
neg_line = []
pos_line = []
for y in dom:
# for each line in plane
if use_f2:
# constant f2()
c = 0.5 * y * y - 2.5 * z * z -11.0375
else:
# constant f()
c = y * y - 5 * z * z + z -16
# find roots
has_solution, neg_root, pos_root = quad_root(a, b, c)
if has_solution:
neg_root = bm.verts.new((neg_root, y, z))
pos_root = bm.verts.new((pos_root, y, z))
pos_line.append(pos_root)
neg_line.append(neg_root)
planes.append([neg_line, pos_line])
# skin it
if len(planes) > 1:
for i in range(len(planes) - 1):
n0, p0 = planes[i]
n1, p1 = planes[i + 1]
# pos
for j in range(len(n0) - 1):
verts = [n0[j], n1[j], n1[j+1], n0[j + 1]]
if None not in verts:
f = bm.faces.new(verts)
verts = [p0[j], p1[j], p1[j+1], p0[j + 1]]
if None not in verts:
f = bm.faces.new(verts)
# make an object
scene = bpy.context.scene
mesh = bpy.data.meshes.new("Quadric")
bm.to_mesh(mesh)
ob = bpy.data.objects.new("Quadric", mesh)
scene.objects.link(ob)
scene.objects.active = ob
ob.scale *= 0.01
To solve for y for f(x, y, z) edit
a, b = 1, 0
c = 2 * x * x - 5 * z * z + z -7 * x - 16
has_solution, neg_root, pos_root = quad_root(a, b, c)
if has_solution:
#print(neg_root, pos_root)
neg_root = bm.verts.new((x, neg_root, z))
pos_root = bm.verts.new((x, pos_root, z))
Result after running in both x and y and joining
