A little while ago Jim Muth commented on his disappointment in the results
he had in the "quat" formula built in to Fractint, suggesting that not
enough degrees of freedom were usable to produce all the 2-D unequivalent
slicing views of the set in an R4 space, and pointing out that the analogous
Julia formula had more degrees of freedom in choice of parameters. He went
on to relate that he used his own formula for exploring the hypercomplex
equivalent, which he later posted. I outlined a possible algorithm to
exhaustively and uniquely specify all the quat slices of this nature, and
have now attempted to implement it for quaternion Julia and Mandelbrot sets.
I would appreciate it if those reading would check the procedure for these
formulae and and see if I have made errors in writing them; there are
settings which do what Jim talked about before, ie. looking either identical
to the well known complex M-set or a slice through a rotation of it, there
are other views very unlike such an object, and there is a strange behavior
for a certain parameter arrangement. So please take a look at this and
comment, and point out any obvious errors I've made in the implementation.
The algorithm for producing the slices is this: choose an arbitrary vector
in R4 specified by p1 and p2 (as four dimensional spherical coordinates);
establish a basis in the 3-space orthogonal to the vector with the vector
head as origin; read a 3-D spherical coordinate direction from p3; choose
the plane normal to this direction passing through the 3-space origin as the
slicing plane and figure 2 orthogonal unit vectors on that plane; then the
formula generates on each pixel pass a quaternion defined as the origin plus
a linear combination of the unit vectors with coefficients as in the parts
of the "pixel" variable. That quaternion is either the value of C, for
M-sets, or the starting value of z, for Julia sets.
The strange behavior is that when a zero value is entered for the length of
the vector rho in real(p1), inside points are plotted properly but when the
length of rho is nonzero, even very small, an outside color is plotted in
the very same region. It doesn't seem reasonable to me that a tiny change
in the origin of the slicing plane would always cause the iterations to
diverge, and only where the inside points used to be. So this is strange.
Also I would appreciate it if someone would check that I used the correct
coefficients in going from the 4-space to the 3-space. Thanks.
Here are the formulae, and 3 pars for the QM-set which if valid seem to show
that it is not simply a rotation of the complex M-set-- the stretching seems
too complicated to me at small scale to be due to rotation alone, but that
may simply be an artifact of the slicing. Comments are welcome.
---Hiram
frm: MandelQuatAnySlice{ ; the quaternion M-set sliced by arbitrary plane in
R4 space. Fractint z-value at each iteration is
; z calculated to be the projection of the quaternion value onto that plane
(to allow orbital images to be viewed through the
; passes=o option). Plane specified by R4 vector spherical coordinates and
; a normal direction in the 3-space orthogonal to the vector:
; p1=(rho(real),omega(-90..90degrees,4-space orientation
angle,"ternitude")),
; p2=(phi4(-90..90deg,"longitude"),theta4(-180..180deg),"latitude"),
; p3=(phi3(-90..90deg,normal long),theta3(-180..180deg,normal lat)),
; p4=first two components of the quaternion perturbation
; p5=last two components of the quaternion perturbation
; z-value each iteration is the projection of the quat onto the plane
IF(initialized==0)
initialized=1
bailout=64
degtorad=pi/180
pert_xy=p4,pert_zw=p5
rho=real(p1),omega4=imag(p1)*degtorad,phi4=real(p2)*degtorad,theta4=imag(p2)*degtorad
phi3=real(p3)*degtorad,theta3=imag(p3)*degtorad
Pnxy=cos(omega4)*cos(phi4)*(cos(theta4)+flip(sin(theta4))) ; unit vector in
direction of point
Pnzw=cos(omega4)*sin(phi4)+flip(sin(omega4))
x4=real(Pnxy),y4=imag(Pnxy),z4=real(Pnzw),w4=imag(Pnzw)
IF( w4==0 && z4==0) ; then 4-comp expression of basis in the 3-space must
use x4 and y4 instead
Bx4uxy= flip(conj(Pnxy)), Bx4uzw= 0 ; unnormalized basis x-unit vector in
the orthogonal 3-space
By4uxy= w4*Pnxy, By4uzw= flip( -|Pnxy| )
Bz4uxy= |Pnxy|*(z4*Pnxy), Bz4uzw= |Pnxy|*(-sqr(w4)-|Pnxy|+flip(z4*w4))
ELSE ; use w4 and z4 for the basis
Bx4uxy= 0, Bx4uzw= flip(conj(Pnzw)),
By4uxy= flip( -|Pnzw| ), By4uzw= y4*Pnzw
Bz4uxy= |Pnzw|*(-sqr(y4)-|Pnzw|+flip(x4*y4)), Bz4uzw= |Pnzw|*(x4*Pnzw)
ENDIF
x3=cos(phi3)*cos(theta3),y3=cos(phi3)*sin(theta3),z3=sin(phi3) ; direction
cosines in the 3-space
IF( z3==0 && y3==0) ; set plane's basis in the 3-space
Bx3uxy=z3, Bx3uz=-x3
By3uxy=x3*y3-flip(sqr(x3)+sqr(z3)),By3uz=y3*z3
ELSE
Bx3uxy=flip(-z3),Bx3uz=y3
By3uxy= -sqr(y3) -sqr(z3) + flip(x3*y3), By3uz=x3*z3
ENDIF
len_x4=sqrt(|Bx4uxy|+|Bx4uzw|),len_y4=sqrt(|By4uxy|+|By4uzw|),len_z4=sqrt(|Bz4uxy|+|Bz4uzw|)
len_x3=sqrt(|Bx3uxy|+|Bx3uz|),len_y3=sqrt(|By3uxy|+|By3uz|)
Bx4xy=Bx4uxy/len_x4,Bx4zw=Bx4uzw/len_x4 ; normalization
By4xy=By4uxy/len_y4,By4zw=By4uzw/len_y4
Bz4xy=Bz4uxy/len_z4,Bz4zw=Bz4uzw/len_z4
Bx3xy=Bx3uxy/len_x3,Bx3z=Bx3uz/len_x3
By3xy=By3uxy/len_y3,By3z=By3uz/len_y3
; sum the contributions of each 3-space component to each 4-space component
of plane's unit axes:
Bxx=real(Bx3xy)*real(Bx4xy)+imag(Bx3xy)*real(By4xy)+Bx3z*real(Bz4xy) ; quat
components of Bx unit vector of
Bxy=real(Bx3xy)*imag(Bx4xy)+imag(Bx3xy)*imag(By4xy)+Bx3z*imag(Bz4xy) ;
planar variation
Bxz=real(Bx3xy)*real(Bx4zw)+imag(Bx3xy)*real(By4zw)+Bx3z*real(Bz4zw)
Bxw=real(Bx3xy)*imag(Bx4zw)+imag(Bx3xy)*imag(By4zw)+Bx3z*imag(Bz4zw)
Byx=real(By3xy)*real(Bx4xy)+imag(By3xy)*real(By4xy)+By3z*real(Bz4xy) ; quat
components of By unit vector of
Byy=real(By3xy)*imag(Bx4xy)+imag(By3xy)*imag(By4xy)+By3z*imag(Bz4xy) ;
planar variation
Byz=real(By3xy)*real(Bx4zw)+imag(By3xy)*real(By4zw)+By3z*real(Bz4zw)
Byw=real(By3xy)*imag(Bx4zw)+imag(By3xy)*imag(By4zw)+By3z*imag(Bz4zw)
Pxy=rho*Pnxy, Pzw=rho*Pnzw ; quat components of the planar origin
ENDIF
Q_xy=pert_xy,Q_zw=pert_zw ; corresponds to z in usual Mandelbrot image
C_xy=Pxy+real(pixel)*Bxx+imag(pixel)*Byx+flip(real(pixel)*Bxy+imag(pixel)*Byy)
; corresponds to usual C
C_zw=Pzw+real(pixel)*Bxz+imag(pixel)*Byz+flip(real(pixel)*Bxw+imag(pixel)*Byw)
:
Q_xy=sqr(Q_xy)-|Q_zw|+C_xy, Q_zw=2*real(Q_xy)*Q_zw+C_zw ; iteration Q->Q^2+C
z=real(Q_xy)*Bxx+imag(Q_xy)*Bxy+real(Q_zw)*Bxz+imag(Q_zw)*Bxw ; dot product
with planar x-direction
z=z+flip(real(Q_xy)*Byx+imag(X_xy)*Byy+real(Q_zw)*Byz+imag(Q_zw)*Byw) ; dot
product with planar y-direction
|Q_xy|+|Q_zw| <= bailout
}
frm: JuliaQuatAnySlice{ ; the quaternion J-set sliced by arbitrary plane in
R4 space. Fractint z-value at each iteration is
; z calculated to be the projection of the quaternion value onto that plane
(to allow orbital images to be viewed through the
; passes=o option). Plane specified by R4 vector spherical coordinates and
; a normal direction in the 3-space orthogonal to the vector:
; p1=(rho(real),omega(-90..90degrees,4-space orientation
angle,"ternitude")),
; p2=(phi4(-90..90deg,"longitude"),theta4(-180..180deg),"latitude"),
; p3=(phi3(-90..90deg,normal long),theta3(-180..180deg,normal lat)),
; p4=first two components of the "C" quaternion for the set
; p5=last two components of the "C" quaternion for the set
; z-value each iteration is the projection of the quat onto the plane
IF(initialized==0)
initialized=1
IF(real(p4)==0)
bailout=64
ELSE
bailout=real(p4)
ENDIF
degtorad=pi/180
rho=real(p1),omega4=imag(p1)*degtorad,phi4=real(p2)*degtorad,theta4=imag(p2)*degtorad
phi3=real(p3)*degtorad,theta3=imag(p3)*degtorad
Pnxy=cos(omega4)*cos(phi4)*(cos(theta4)+flip(sin(theta4))) ; unit vector in
direction of point
Pnzw=cos(omega4)*sin(phi4)+flip(sin(omega4))
x4=real(Pnxy),y4=imag(Pnxy),z4=real(Pnzw),w4=imag(Pnzw)
IF( w4==0 && z4==0) ; then 4-comp expression of basis in the 3-space must
use x4 and y4 instead
Bx4uxy= flip(conj(Pnxy)), Bx4uzw= 0 ; unnormalized basis x-unit vector in
the orthogonal 3-space
By4uxy= w4*Pnxy, By4uzw= flip( -|Pnxy| )
Bz4uxy= |Pnxy|*(z4*Pnxy), Bz4uzw= |Pnxy|*(-sqr(w4)-|Pnxy|+flip(z4*w4))
ELSE ; use w4 and z4 for the basis
Bx4uxy= 0, Bx4uzw= flip(conj(Pnzw)),
By4uxy= flip( -|Pnzw| ), By4uzw= y4*Pnzw
Bz4uxy= |Pnzw|*(-sqr(y4)-|Pnzw|+flip(x4*y4)), Bz4uzw= |Pnzw|*(x4*Pnzw)
ENDIF
x3=cos(phi3)*cos(theta3),y3=cos(phi3)*sin(theta3),z3=sin(phi3) ; direction
cosines in the 3-space
IF( z3==0 && y3==0) ; set plane's basis in the 3-space
Bx3uxy=z3, Bx3uz=-x3
By3uxy=x3*y3-flip(sqr(x3)+sqr(z3)),By3uz=y3*z3
ELSE
Bx3uxy=flip(-z3),Bx3uz=y3
By3uxy= -sqr(y3) -sqr(z3) + flip(x3*y3), By3uz=x3*z3
ENDIF
len_x4=sqrt(|Bx4uxy|+|Bx4uzw|),len_y4=sqrt(|By4uxy|+|By4uzw|),len_z4=sqrt(|Bz4uxy|+|Bz4uzw|)
len_x3=sqrt(|Bx3uxy|+|Bx3uz|),len_y3=sqrt(|By3uxy|+|By3uz|)
Bx4xy=Bx4uxy/len_x4,Bx4zw=Bx4uzw/len_x4 ; normalization
By4xy=By4uxy/len_y4,By4zw=By4uzw/len_y4
Bz4xy=Bz4uxy/len_z4,Bz4zw=Bz4uzw/len_z4
Bx3xy=Bx3uxy/len_x3,Bx3z=Bx3uz/len_x3
By3xy=By3uxy/len_y3,By3z=By3uz/len_y3
; sum the contributions of each 3-space component to each 4-space component
of plane's unit axes:
Bxx=real(Bx3xy)*real(Bx4xy)+imag(Bx3xy)*real(By4xy)+Bx3z*real(Bz4xy) ; quat
components of Bx unit vector of
Bxy=real(Bx3xy)*imag(Bx4xy)+imag(Bx3xy)*imag(By4xy)+Bx3z*imag(Bz4xy) ;
planar variation
Bxz=real(Bx3xy)*real(Bx4zw)+imag(Bx3xy)*real(By4zw)+Bx3z*real(Bz4zw)
Bxw=real(Bx3xy)*imag(Bx4zw)+imag(Bx3xy)*imag(By4zw)+Bx3z*imag(Bz4zw)
Byx=real(By3xy)*real(Bx4xy)+imag(By3xy)*real(By4xy)+By3z*real(Bz4xy) ; quat
components of By unit vector of
Byy=real(By3xy)*imag(Bx4xy)+imag(By3xy)*imag(By4xy)+By3z*imag(Bz4xy) ;
planar variation
Byz=real(By3xy)*real(Bx4zw)+imag(By3xy)*real(By4zw)+By3z*real(Bz4zw)
Byw=real(By3xy)*imag(Bx4zw)+imag(By3xy)*imag(By4zw)+By3z*imag(Bz4zw)
Pxy=rho*Pnxy, Pzw=rho*Pnzw ; quat components of the planar origin
C_xy=p4,C_zw=p5
ENDIF
Q_xy=Pxy+real(pixel)*Bxx+imag(pixel)*Byx+flip(real(pixel)*Bxy+imag(pixel)*Byy)
; corresponds to usual C
Q_zw=Pzw+real(pixel)*Bxz+imag(pixel)*Byz+flip(real(pixel)*Bxw+imag(pixel)*Byw)
:
Q_xy=sqr(Q_xy)-|Q_zw|+C_xy, Q_zw=2*real(Q_xy)*Q_zw+C_zw ; iteration Q->Q^2+C
z=real(Q_xy)*Bxx+imag(Q_xy)*Bxy+real(Q_zw)*Bxz+imag(Q_zw)*Bxw ; dot product
with planar x-direction
z=z+flip(real(Q_xy)*Byx+imag(X_xy)*Byy+real(Q_zw)*Byz+imag(Q_zw)*Byw) ; dot
product with planar y-direction
|Q_xy|+|Q_zw| <= bailout
}
par: NotAMinibrot {
reset=2004 type=formula formulafile=quat.frm
formulaname=mandelquatanyslice passes=1
center-mag=-0.85477129473239550/+0.40709301723382000/266.6654
params=0/0/0/90/15/-60/0/0/0/0 float=y maxiter=15000 inside=0
outside=summ periodicity=0 rseed=28455 cyclerange=2/253
colors=000Gfrww0<3>wY0wR0wK0<3>UCU<3>Y6jZ4n_2ra0w<14>P8tO8tN9t<3>JBtHBtG\
CsECsCCs<2>CIsDJsDJs<44>girgirhjrhjrikr<9>qqtqrtrrt<3>uuuvvvvwv<2>yyvzzw\
zyq<3>xuRwtKwtK<3>wsKwsKwsKvrK<3>uiHugGudFtbE<6>sOCsMCsKC<2>sECsCCqCC<2>\
kCCiCCiCC<39>fClfCmfCnfCofCp<4>cFrbGraHr`Ir_Jr_Js<17>JasIbsHcsGdsFes<3>D\
ksCmsCosCqsCss<10>CsY
}
par: AParticularView { ; of the QM-set, using outside=summ to show the
probabilit
; y that the distortion of the M-set in the slice is far
; more complicated than likely achievable with a simple
; rotation.
reset=2004 type=formula formulafile=quat.frm
formulaname=MandelQuatAnySlice passes=1 center-mag=0/0/0.6666667
params=0/0/75/135/-135/-24/0/0 float=y maxiter=600 inside=0
outside=summ logmap=yes periodicity=0 rseed=28455 cyclerange=2/253
colors=zzz000ww0<3>wY0wR0wK0<3>UCU<3>Y6jZ4n_2ra0w<16>N9tM9tLAtKAtJBt<2>E\
CsCCsCEsCGsCIs<41>cfrdfregregrfhrgirgir<3>jkrjlsklp<3>ppaqpZsqVtrSusOwtK\
<17>wtKwsKwsK<2>wsKvrKvpK<3>ugGudFtbEt`D<5>sOCsMCsKC<2>sECsCCqCC<2>kCCiC\
CiCC<40>fCmfCnfCofCpfCq<10>YLsXMsWNs<14>IbsHcsGdsFesEfs<3>CmsCosCqsCss<1\
0>CsY
}
par: ExtremeStretching { ; happens in parts of this slice of the Quaternion
M-set
reset=2004 type=formula formulafile=quat.frm
formulaname=mandelquatanyslice passes=1
center-mag=-0.484191/-0.633079/0.3625693
params=0/0/0/90/15/-60/0/0/0/0 float=y maxiter=1200 inside=0
outside=summ periodicity=0 rseed=28455 cyclerange=2/253
colors=000GfrVXe<2>B5O<12>9_b9ac9de<2>9kh8nj8nj<42>gGIhGHiFH<2>lCElCElDE\
<54>bsNbtNbuNbvNawO<50>rffsegseg<2>tdhtdhsdg<31>kjZkjZkjZ<3>ikXd_P<3>LjR\
GmSGmS<3>HmTHmUHmUImVlte<5>nwwhnqaekCs_CsY
}