Okay, here’s my code. I managed to use the raw format, so either an older build or a current one will render it. I tested it with sin(x) + sin(y) + sin(z), and with spherical functions. To use, just take the output file (I keep all of my scripts in C:\Python31, so my file ended up there) and load it as an 8 bit raw, and set the three resolution parameters to match ‘res’. Play with the density till it’s pretty. 
Oh, and I ran it in IDLE, but I would suppose that blender could do it, too.
To get a isosurface, just play with the color ramp. It should work, depending on the resolutions.
# 4D Function Plotter
# For use with Blender Volumetrics (2.5 prerelease version)
# Henry Schreiner
# Python 3.1
from math import *
class mat4d:
'''Make a 4D matrix with a given function'''
def __init__(self,
func: 'Has to have 3 variables'=lambda x,y,z: x+y+z,
res=9, rang=5, sph=False):
'''This uses x,y,z; or it uses r, phi, and theta if shp is True'''
(self.func,self.rang,self.sph,self.res)=(func,rang,sph,res)
self.M = [0]*res**3
ipos = lambda a: (a/(res-1)-.5)*2*rang
for n in range(res**3):
(x,y,z)=map(ipos,self.num2ind(n))
if sph: (x,y,z)=self.rec2sph(x,y,z)
try: self.M[n]=func(x,y,z)
except ZeroDivisionError: self.M[n]=0
def save(self, filename):
'''This will save a .raw.gz file for Blender RAW volumetrics'''
fob = open(filename + '.raw.gz','bw')
fob.write(bytes(map(round,self.M)))
fob.close()
def num2ind(self,a):
return (int(a/self.res**2),
int((a%self.res**2)/self.res),
a%self.res)
def ind2num(self,i,j,k):
return i*self.res**2 + j*self.res + k
def normalize(self, absolute=False, scale =255):
'''Normalizes to 0-scl, uses absolute values or total'''
minimum = min(map(abs,self.M) if absolute else self.M)
rang = max(map(abs,self.M) if absolute else self.M) - minimum
for n in range(self.res**3):
self.M[n]=((abs(self.M[n]) if absolute else self.M[n])
-minimum)/rang*scale
@staticmethod
def rec2sph(x,y,z):
'''Returns (r, phi, theta) for x,y,x'''
r= sqrt(x**2 + y**2 + z**2)
p = atan2(y,x)
t = (acos(z/r) if r else 0.0) # quick and dirty fix for r=0
return (r, p, t)
## Example usage
# Setting up wavefunction arrays
psi = [lambda r: 0, #n=0
lambda r: 1/(r*sqrt(2*pi))*sin(pi*r), #n=1
lambda r: 1/(r*sqrt(2*pi))*sin(2*pi*r)] #n=2
Y = [[lambda p,t: (1/(4*pi))**(1/2)], #lm=00
[lambda p,t: (3/(4*pi))**(1/2)*cos(t), #lm=10
lambda p,t: -(3/(8*pi))**(1/2)*sin(t)*e**(p*1j)], #lm=11
[lambda p,t: (5/(16*pi))**(1/2)*(3*cos(t)**2-1), #lm=20
lambda p,t: -(15/(8*pi))**(1/2)*sin(t)*cos(t)*e**(p*1j),#lm=21
lambda p,t: (15/(32*pi))**(1/2)*sin(t)**2*e**(p*2j)]] #lm=22
# Example using nlm=222
def wavefn(r,p,t):
(n,l,m)=(2,2,2)
if r>1: return 0.0
else: return (Y[l][m](p,t)*Y[l][m](p,t).conjugate()).real * psi[n](r)**2
if __name__ == '__main__':
import time
t = time.time()
# !!! This is how you use it!!!
q=mat4d(wavefn,res=32,rang=1,sph=True) # 1st parameter = function
q.normalize() # scales numbers to 0-255, for 8 bit use
q.save('wvfn222') # saves it as the name
t = time.time()-t
print('The time taken was:', round(t,3), 'seconds.')
Attachments
