bugman


« Reply #28 on: November 18, 2009, 06:11:17 PM » 

OK, well here are my first attempts to create 3D Glynn fractals using the triplex formulas. The first image is a slice that was ray traced using the escape time method. But it doesn't show any tree branch details so it's pretty boring. The other images are point clouds that were calculated using the MIIM method. Finding all the unique valid roots was a challenge. In general, I found that there are at least 2 roots when x < 0 and at least one root when x > 0. But there is also a strange conic section region containing 3 roots when z² > x² + y². It should also be noted that convergence is very sensitive to the initial seed value. I found the best results when the seed value is inside one of the cones of the conic section. But this only finds one half of the overall fractal. You can find the other half by placing the seed in the other cone, but I preferred to leave this part out, because it gives a nice crosssection so you can see inside the fractal. In these images c={0.2,0,0}.
Here is some Mathematica code: (* runtime: 15 seconds *) power = 1.5; norm[x_] := x.x; TriplexPow[{x_, y_, z_}, n_] := Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi},phi = n ArcSin[z/r]; r^n{Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}]; NumRootsGlynn[{x_, y_, z_}] := If[z^2 > x^2 + y^2, 3, If[x < 0, 2, 1]]; TriplexRootGlynn[{x_, y_, z_}, k_] := Module[{r = Sqrt[x^2 + y^2 + z^2], ktheta, kphi, theta, phi}, {ktheta, kphi} = If[k == 0, {0, 0}, If[z^2 > x^2 + y^2, If[k == 1, If[x <0 && y < 0, {2, 0}, {0.5, If[z > 0, 0.5, 2.5]}], If[x < 0 && y > 0, {1, 0}, {2.5, If[z > 0, 0.5, 2.5]}]], {If[y < 0, 2, 1], 0}]]; theta = (ArcTan[x, y] + ktheta Pi)/power; phi = (ArcSin[z/r] +kphi Pi)/power; r^(1/power){Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}]; c = {0.2, 0, 0}; imax = 100; dpmax = 50.0; p = Table[{0.61, 0.0, 0.42}, {imax}]; dp = Table[1, {imax}]; roots = Table[0, {imax}]; plist = {}; i = 2; While[i > 1 && Length[plist] < 10000, p[[ i]] = TriplexRootGlynn[p[[ i  1]]  c, roots[[ i]]]; dp[[ i]] = Abs[power]Sqrt[norm[TriplexPow[p[[ i]], power  1]]]dp[[ i  1]];plist = Append[plist, p[[ i]]]; prune = (i == imax  dp[[ i]] > dpmax); If[prune, While[i > 1 && roots[[ i]] == NumRootsGlynn[p[[ i  1]]  c]  1, roots[[ i]] = 0; i]; roots[[ i]]++, i++; roots[[ i]] = 0]]; Show[Graphics3D[{PointSize[0.005], Point /@ plist}]]
For those of you who don't use Mathematica, here is some C++ code: //Glynn roots: {x,y,z}^(1/1.5) = {x,y,z}^(2/3) int NumRootsGlynn(triplex p) { // finds number of branches for a given triplex p bool cone=(sqr(p.z)>sqr(p.x)+sqr(p.y)); return (cone?3:(p.x<0.0?2:1)); } triplex RootGlynn(triplex p, int k) {double n=1.5; bool cone=(sqr(p.z)>sqr(p.x)+sqr(p.y)); double ktheta, kphi; if(k==0) {ktheta=kphi=0.0;} else { if(!cone) {ktheta=(p.y<0.0?2.0:1.0); kphi=0.0;} else { if(k==1) { if(p.x<0.0 && p.y<0.0) {ktheta=2.0; kphi=0.0;} else {ktheta=0.5; kphi=(p.z>0.0?0.5:2.5);} } else { if(p.x<0.0 && p.y>0.0) {ktheta=1.0; kphi=0.0;} else {ktheta=2.5; kphi=(p.z>0.0?0.5:2.5);} } } } double r=abs(p), theta=(atan2(p.y,p.x)+ktheta*pi)/n, phi=(asin(p.z/r)+kphi*pi)/n; double cosphi=cos(phi); return pow(r,1.0/n)*triplex(cos(theta)*cosphi,sin(theta)*cosphi,sin(phi)); }
