Long Beach: adjacent the Queen Mary (former Spruce Goose garage)
Another take on linear algebra involves importing from what we call Martian Math. A goal is to connect both to science fiction and to the notion of something alien, as in unfamiliar. The coordinate system and vector apparatus shared below is anything but mainstream.
We call it Quadrays. The picture below is a screen shot, so clicking on the links only takes you to the archive for the picture.
import qrays
from qrays import Qvector as Martian
from qrays import Vector
from itertools import permutations as perm
dir_4 = {spoke for spoke in perm((0,0,0,1))} # dupeless
dir_4
{(0, 0, 0, 1), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0)}
Four rays splay from the origin to carve space into four quadrants. Linear combinations of vectors in these four directions span all of space.
These four elementary rays need not be unit length. Rather, the bounding tetrahedron has unit edges and unit volume.
In a lowest terms representation of any point, at least one coordinate is always zero, the ray not directly bounding the quadrant wherein the point resides. The remaining three vectors add together positively to reach the point i.e. quadray coordinates in lowest terms are always non-negative.
a,b,c,d = [Martian(coords) for coords in dir_4]
a + b + c + d
ivm_vector(a=0, b=0, c=0, d=0)
a.length()
0.6123724356957945
(a - b).length()
1.0
a.angle(b)
import sympy as sp
root2 = sp.sqrt(2)
import qrays
import importlib
importlib.reload(qrays)
from qrays import Vector, Qvector as Martian
def to_xyz(a, b, c, d):
k = root2/4
x = k * (a - b - c + d)
y = k * (a - b + c - d)
z = k * (a + b - c - d)
return x,y,z
x,y,z = to_xyz(1,0,0,0) # convert from Martian
(x,y,z)
(sqrt(2)/4, sqrt(2)/4, sqrt(2)/4)
def to_qray(x, y, z):
"""return (a, b, c, d) quadray based on current (x, y, z)"""
k = root2
x_ge_0 = 1 if x >=0 else 0
y_ge_0 = 1 if y >=0 else 0
z_ge_0 = 1 if z >=0 else 0
x_lt_0 = 1 if x < 0 else 0
y_lt_0 = 1 if y < 0 else 0
z_lt_0 = 1 if z < 0 else 0
a = k * (x_ge_0 * x + y_ge_0 * y + z_ge_0 * z)
b = k * (x_lt_0 * -x + y_lt_0 * -y + z_ge_0 * z)
c = k * (x_lt_0 * -x + y_ge_0 * y + z_lt_0 * -z)
d = k * (x_ge_0 * x + y_lt_0 * -y + z_lt_0 * -z)
# put in canonical form
m = min((a,b,c,d))
a,b,c,d = a-m, b-m, c-m, d-m
return (a, b, c, d)
to_qray(x, y, z) # convert to Martian
(1, 0, 0, 0)
v = Vector((x,y,z))
v
xyz_vector(x=sqrt(2)/4, y=sqrt(2)/4, z=sqrt(2)/4)
q = v.quadray() # convert to Martian
q
ivm_vector(a=1, b=0, c=0, d=0)
q.xyz
xyz_vector(x=0.25*sqrt(2), y=0.25*sqrt(2), z=0.25*sqrt(2))
Two of one quadray, added to one of two others, with a fourth set to zero, equals one of the twelve directions from a ball center to a neighboring ball center within the IVM (isotropic vector matrix). The ball packing in question is known as the CCP (cubic close packing).
dir_12 = {spoke for spoke in perm((0,1,1,2))} # dupeless
dir_12
{(0, 1, 1, 2), (0, 1, 2, 1), (0, 2, 1, 1), (1, 0, 1, 2), (1, 0, 2, 1), (1, 1, 0, 2), (1, 1, 2, 0), (1, 2, 0, 1), (1, 2, 1, 0), (2, 0, 1, 1), (2, 1, 0, 1), (2, 1, 1, 0)}
from pov import POV_Vector, pov_header
beacon = [Martian((t)) for t in dir_12]
xyz_beacon = [POV_Vector(v.xyz.x, v.xyz.y, v.xyz.z) for v in beacon]
# POV-Ray
edge_color = "rgb <1, 0.4, 0>"
edge_radius= 0.03
vert_color = "rgb <0, 0, 1>"
vert_radius= 0.05
with open("spokes.pov", 'w') as output:
print(pov_header, file=output) # print to the output file
for v in xyz_beacon:
v.draw_edge(edge_color, edge_radius, output)
v.draw_vert(vert_color, vert_radius, output)
print("spokes.pov ready for rendering")
spokes.pov ready for rendering
o,p,q,r,s,t,u,v,w,x,y,z = [Martian((t)) for t in dir_12]
In the case of the XYZ coordinate system, three vectors from the origin, in x, y and z directions, define a 90-90-90 degree corner. If the lengths of the three vectors are a, b, c then the volume of the resulting parallelopiped is their product.
In the case of the IVM, consider any three vectors defining a 60-60-60 degree corner with lengths a, b and c. Their volume is the resulting tetrahedron, and is likewise their product.
origin = Martian((0,0,0,0))
The absolute value of the determinant of a matrix, times 1/4, gives the volume, in tetravolumes, of the tetrahedron determined by four quadrays (a, b, c, d), each with four coordinates (e.g. a0, a1, a2, a3).
Lets find three vectors from our 12 above, that form a 60-60-60 corner. We may then scale each to a different length and compute the resulting volume. The quadray to the origin, (0,0,0,0), will be one of the four points.
o.angle(s)
s.angle(v)
o.angle(v)
from sympy import Matrix
# try varying scale factors
e0, e1, e2 = o*3, s*3, v*3 # per picture: 2, 2, 5
The four corners of a tetrahedron with a 60-60-60 degree corner at the origin, e0, e1, e2.
uvt = Matrix([[*origin.coords,1], # origin
[*e0.coords, 1], # e0
[*e1.coords, 1], # e1
[*e2.coords, 1], # e2
[1,1,1,1,0]])
uvt
abs(uvt.det())/4
s.area(v)
Now that we have the notion of "tetravolumes" firmly nailed down, we're able to create a new volumes table:
# try varying scale factors
e0, e1, e2 = o, s, v # 1, 1, 1
def tetravol(v0, v1, v2):
uvt = Matrix([[*origin.coords,1], # origin
[*v0.coords, 1], # e0
[*v1.coords, 1], # e1
[*v2.coords, 1], # e2
[1,1,1,1,0]])
return abs(uvt.det())/4
def frustrum(start, end):
s0 = o*start
s1 = s*start
s2 = v*start
e0 = o*end
e1 = s*end
e2 = v*end
return tetravol(e0,e1,e2) - tetravol(s0, s1, s2)
tip = tetravol(e0, e1, e2)
tip
frust = frustrum(2, 4)
frust
frustrum(0, 1) + frustrum(1, 2) + frustrum(2, 3) + frustrum(3, 4) + frustrum(4, 5)
import pandas as pd
import numpy as np
notches = np.linspace(0, 5, 101)[:-1]
notches
array([0. , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 , 0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1. , 1.05, 1.1 , 1.15, 1.2 , 1.25, 1.3 , 1.35, 1.4 , 1.45, 1.5 , 1.55, 1.6 , 1.65, 1.7 , 1.75, 1.8 , 1.85, 1.9 , 1.95, 2. , 2.05, 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 , 2.55, 2.6 , 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05, 3.1 , 3.15, 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 , 3.65, 3.7 , 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15, 4.2 , 4.25, 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55, 4.6 , 4.65, 4.7 , 4.75, 4.8 , 4.85, 4.9 , 4.95])
delta = notches[1]-notches[0]
delta
0.05
df = pd.DataFrame({"notch":notches,
"next_notch":notches+delta})
df.head()
notch | next_notch | |
---|---|---|
0 | 0.00 | 0.05 |
1 | 0.05 | 0.10 |
2 | 0.10 | 0.15 |
3 | 0.15 | 0.20 |
4 | 0.20 | 0.25 |
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 100 entries, 0 to 99 Data columns (total 2 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 notch 100 non-null float64 1 next_notch 100 non-null float64 dtypes: float64(2) memory usage: 1.7 KB
df["Frustrum"] = [float(frustrum(df.notch[i], df.next_notch[i])) for i in df.index]
df
notch | next_notch | Frustrum | |
---|---|---|---|
0 | 0.00 | 0.05 | 0.000125 |
1 | 0.05 | 0.10 | 0.000875 |
2 | 0.10 | 0.15 | 0.002375 |
3 | 0.15 | 0.20 | 0.004625 |
4 | 0.20 | 0.25 | 0.007625 |
... | ... | ... | ... |
95 | 4.75 | 4.80 | 3.420125 |
96 | 4.80 | 4.85 | 3.492125 |
97 | 4.85 | 4.90 | 3.564875 |
98 | 4.90 | 4.95 | 3.638375 |
99 | 4.95 | 5.00 | 3.712625 |
100 rows × 3 columns
df.Frustrum.cumsum()
0 0.000125 1 0.001000 2 0.003375 3 0.008000 4 0.015625 ... 95 110.592000 96 114.084125 97 117.649000 98 121.287375 99 125.000000 Name: Frustrum, Length: 100, dtype: float64
df["Cumulative"] = df.Frustrum.cumsum()
df
notch | next_notch | Frustrum | Cumulative | |
---|---|---|---|---|
0 | 0.00 | 0.05 | 0.000125 | 0.000125 |
1 | 0.05 | 0.10 | 0.000875 | 0.001000 |
2 | 0.10 | 0.15 | 0.002375 | 0.003375 |
3 | 0.15 | 0.20 | 0.004625 | 0.008000 |
4 | 0.20 | 0.25 | 0.007625 | 0.015625 |
... | ... | ... | ... | ... |
95 | 4.75 | 4.80 | 3.420125 | 110.592000 |
96 | 4.80 | 4.85 | 3.492125 | 114.084125 |
97 | 4.85 | 4.90 | 3.564875 | 117.649000 |
98 | 4.90 | 4.95 | 3.638375 | 121.287375 |
99 | 4.95 | 5.00 | 3.712625 | 125.000000 |
100 rows × 4 columns
df["Difference"] = df.Frustrum.diff()
df
notch | next_notch | Frustrum | Cumulative | Difference | |
---|---|---|---|---|---|
0 | 0.00 | 0.05 | 0.000125 | 0.000125 | NaN |
1 | 0.05 | 0.10 | 0.000875 | 0.001000 | 0.00075 |
2 | 0.10 | 0.15 | 0.002375 | 0.003375 | 0.00150 |
3 | 0.15 | 0.20 | 0.004625 | 0.008000 | 0.00225 |
4 | 0.20 | 0.25 | 0.007625 | 0.015625 | 0.00300 |
... | ... | ... | ... | ... | ... |
95 | 4.75 | 4.80 | 3.420125 | 110.592000 | 0.07125 |
96 | 4.80 | 4.85 | 3.492125 | 114.084125 | 0.07200 |
97 | 4.85 | 4.90 | 3.564875 | 117.649000 | 0.07275 |
98 | 4.90 | 4.95 | 3.638375 | 121.287375 | 0.07350 |
99 | 4.95 | 5.00 | 3.712625 | 125.000000 | 0.07425 |
100 rows × 5 columns
0.002375 - 0.000875
0.0015
df["DiffDiff"] = df.Difference.diff()
df
notch | next_notch | Frustrum | Cumulative | Difference | DiffDiff | |
---|---|---|---|---|---|---|
0 | 0.00 | 0.05 | 0.000125 | 0.000125 | NaN | NaN |
1 | 0.05 | 0.10 | 0.000875 | 0.001000 | 0.00075 | NaN |
2 | 0.10 | 0.15 | 0.002375 | 0.003375 | 0.00150 | 0.00075 |
3 | 0.15 | 0.20 | 0.004625 | 0.008000 | 0.00225 | 0.00075 |
4 | 0.20 | 0.25 | 0.007625 | 0.015625 | 0.00300 | 0.00075 |
... | ... | ... | ... | ... | ... | ... |
95 | 4.75 | 4.80 | 3.420125 | 110.592000 | 0.07125 | 0.00075 |
96 | 4.80 | 4.85 | 3.492125 | 114.084125 | 0.07200 | 0.00075 |
97 | 4.85 | 4.90 | 3.564875 | 117.649000 | 0.07275 | 0.00075 |
98 | 4.90 | 4.95 | 3.638375 | 121.287375 | 0.07350 | 0.00075 |
99 | 4.95 | 5.00 | 3.712625 | 125.000000 | 0.07425 | 0.00075 |
100 rows × 6 columns
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 100 entries, 0 to 99 Data columns (total 6 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 notch 100 non-null float64 1 next_notch 100 non-null float64 2 Frustrum 100 non-null float64 3 Cumulative 100 non-null float64 4 Difference 99 non-null float64 5 DiffDiff 98 non-null float64 dtypes: float64(6) memory usage: 4.8 KB
df.plot(x="next_notch", y=["Frustrum", "Difference", "DiffDiff"]);
df.plot(x="next_notch", y=["Frustrum", "Difference"]);
frustrum(0, 1) + frustrum(1, 5)
frustrum(1, 5)
In the explorations below, we're continuing to slice our tetrahedron parallel to any face.
from sympy.utilities.lambdify import lambdify
pd.set_option("display.precision", 10)
X = sp.Symbol('X')
notches = np.linspace(0, 5, 101)
notches
array([0. , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 , 0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1. , 1.05, 1.1 , 1.15, 1.2 , 1.25, 1.3 , 1.35, 1.4 , 1.45, 1.5 , 1.55, 1.6 , 1.65, 1.7 , 1.75, 1.8 , 1.85, 1.9 , 1.95, 2. , 2.05, 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 , 2.55, 2.6 , 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05, 3.1 , 3.15, 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 , 3.65, 3.7 , 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15, 4.2 , 4.25, 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55, 4.6 , 4.65, 4.7 , 4.75, 4.8 , 4.85, 4.9 , 4.95, 5. ])
linear = 0.015 * X - 0.00075
quadratic = 0.15*X**2 - 0.0075*X + 0.000125
third_power = X**3
df2 = pd.DataFrame({"notch":notches})
df2
notch | |
---|---|
0 | 0.00 |
1 | 0.05 |
2 | 0.10 |
3 | 0.15 |
4 | 0.20 |
... | ... |
96 | 4.80 |
97 | 4.85 |
98 | 4.90 |
99 | 4.95 |
100 | 5.00 |
101 rows × 1 columns
df2["Frustrum"] = [float(frustrum(df2.notch[i]-0.05, df2.notch[i])) for i in df2.index]
df2.iloc[0,1]=0
df2
notch | Frustrum | |
---|---|---|
0 | 0.00 | 0.000000 |
1 | 0.05 | 0.000125 |
2 | 0.10 | 0.000875 |
3 | 0.15 | 0.002375 |
4 | 0.20 | 0.004625 |
... | ... | ... |
96 | 4.80 | 3.420125 |
97 | 4.85 | 3.492125 |
98 | 4.90 | 3.564875 |
99 | 4.95 | 3.638375 |
100 | 5.00 | 3.712625 |
101 rows × 2 columns
linear_f = lambdify(X, linear, 'numpy')
df2['linear'] = linear_f(df2.notch)
df2
notch | Frustrum | linear | |
---|---|---|---|
0 | 0.00 | 0.000000 | -0.00075 |
1 | 0.05 | 0.000125 | 0.00000 |
2 | 0.10 | 0.000875 | 0.00075 |
3 | 0.15 | 0.002375 | 0.00150 |
4 | 0.20 | 0.004625 | 0.00225 |
... | ... | ... | ... |
96 | 4.80 | 3.420125 | 0.07125 |
97 | 4.85 | 3.492125 | 0.07200 |
98 | 4.90 | 3.564875 | 0.07275 |
99 | 4.95 | 3.638375 | 0.07350 |
100 | 5.00 | 3.712625 | 0.07425 |
101 rows × 3 columns
df2['DiffFrustrum'] = df2.Frustrum.diff()
df2
notch | Frustrum | linear | DiffFrustrum | |
---|---|---|---|---|
0 | 0.00 | 0.000000 | -0.00075 | NaN |
1 | 0.05 | 0.000125 | 0.00000 | 0.000125 |
2 | 0.10 | 0.000875 | 0.00075 | 0.000750 |
3 | 0.15 | 0.002375 | 0.00150 | 0.001500 |
4 | 0.20 | 0.004625 | 0.00225 | 0.002250 |
... | ... | ... | ... | ... |
96 | 4.80 | 3.420125 | 0.07125 | 0.071250 |
97 | 4.85 | 3.492125 | 0.07200 | 0.072000 |
98 | 4.90 | 3.564875 | 0.07275 | 0.072750 |
99 | 4.95 | 3.638375 | 0.07350 | 0.073500 |
100 | 5.00 | 3.712625 | 0.07425 | 0.074250 |
101 rows × 4 columns
df2['Total_Volume'] = df2.Frustrum.cumsum()
df2
notch | Frustrum | linear | DiffFrustrum | Total_Volume | |
---|---|---|---|---|---|
0 | 0.00 | 0.000000 | -0.00075 | NaN | 0.000000 |
1 | 0.05 | 0.000125 | 0.00000 | 0.000125 | 0.000125 |
2 | 0.10 | 0.000875 | 0.00075 | 0.000750 | 0.001000 |
3 | 0.15 | 0.002375 | 0.00150 | 0.001500 | 0.003375 |
4 | 0.20 | 0.004625 | 0.00225 | 0.002250 | 0.008000 |
... | ... | ... | ... | ... | ... |
96 | 4.80 | 3.420125 | 0.07125 | 0.071250 | 110.592000 |
97 | 4.85 | 3.492125 | 0.07200 | 0.072000 | 114.084125 |
98 | 4.90 | 3.564875 | 0.07275 | 0.072750 | 117.649000 |
99 | 4.95 | 3.638375 | 0.07350 | 0.073500 | 121.287375 |
100 | 5.00 | 3.712625 | 0.07425 | 0.074250 | 125.000000 |
101 rows × 5 columns
Lets construct the unit XYZ cube with normal vectors, convert them to Quadrays, and feed them to our tetravolumes formula.
To express the result in cubic volumes, we will need to convert out of tetravolumes using the Synergetics Constant S3 i.e. $\sqrt{9/8}$.
The relationship between XYZ and IVM that the Martians + Earthings have constructed (in the science fiction story behind Martian Math), assumes a CCP ball size in common, i.e. an IVM ball of radius R, diameter D.
The unit XYZ cube has edges R, whereas the IVM tetrahedron has edges D. Nevertheless, the $R^{3}$ cube has volume greater than $D^{3}$ by a scale factor of S3.
In the qrays module, the edges between any two elementary quadray tips is unity (the unity-2 of 2R). With respect to 2R, the XYZ cube has edges half that length, or R. So the X, Y and Z vectors of length R get entered with length 1/2 with respect to the IVM prime vector of length 1.
ex = qrays.Vector((1/2,0,0))
ey = qrays.Vector((0,1/2,0))
ez = qrays.Vector((0,0,1/2))
ex.angle(ey).evalf()
ex.angle(ez).evalf()
ey.angle(ez).evalf()
qx = ex.quadray()
qy = ey.quadray()
qz = ez.quadray()
qx.angle(qy).evalf()
S3 = sp.sqrt(9/8)
S3.evalf(30)
def corner_vol(v0, v1, v2):
uvt = Matrix([[*origin.coords,1], # origin
[*v0.coords, 1], # v0
[*v1.coords, 1], # v1
[*v2.coords, 1], # v2
[1,1,1,1,0]])
return abs(uvt.det())/4
def tetra_vol(v0, v1, v2, v3):
uvt = Matrix([[*v0.coords,1], # v0
[*v1.coords, 1], # v1
[*v2.coords, 1], # v2
[*v3.coords, 1], # v3
[1,1,1,1,0]])
return abs(uvt.det())/4
Lets test our volume function with the original basic Quadray tetrahedron of edges D and tetravolume 1.
qa = Martian((1,0,0,0))
qb = Martian((0,1,0,0))
qc = Martian((0,0,1,0))
qd = Martian((0,0,0,1))
tetra_vol(qa, qb, qc, qd)
Confirm or qx, qy, qz are the expected normal vectors of length R.
qx.xyz
xyz_vector(x=0.500000000000000, y=0, z=0)
qy.xyz
xyz_vector(x=0, y=0.500000000000000, z=0)
qz.xyz
xyz_vector(x=0, y=0, z=0.500000000000000)
The 90-90-90 tip tetrahedron is one 1/4 of 2/3 of the total cube, i.e. 4 such tips apply to an internal tetrahedron of 1/3 the total cube's volume. I.E. 1/6th the total cube.
So multiply the result we get, for the cube's tetrahedron tip (closing the lid on qx, qy, qz), by 6 to get the total cube volume.
(corner_vol(qx, qy, qz) * 6).evalf()
But we're still in tetravolumes.
It takes fewer XYZ unit cubes (edges R) than IVM unit tetrahedrons (edges D) to fill the same volume, i.e. the XYZ unit cube is bigger by a scale factor of S3.
So use 1/S3 when going from tetravolumes to cubic volumes i.e. it takes fewer of the latter so your constant is < 1.
(corner_vol(qx, qy, qz) * 6 * 1/S3).evalf()