pgd_UnrollCardioid { ; Unrolls the cardioid of the M-set into a straight line. transform: #pixel = (1 - sqr((0,1)*(4/(#pixel-(0,2))-(0,1))-1))/4 default: title = "Unroll Cardioid" } pgd_ComponentCenter { ; Centers the fractal on a root of (((c^2+c)^2+c)^2+c ... = 0. ; This will be the middle of a bud or minibrot if the fractal ; is the default Mandelbrot fractal. ; Newton's method is used to locate the root. If it fails to ; converge, the fractal will end up centered on the origin. ; The component centered will have a period that divides the ; "Order of component" parameter. global: complex s = @seed complex os bool found = false float p = 0.1 int i = 0 WHILE (i < @numi) p = p * 0.1 i = i + 1 ENDWHILE float q = 1/p i = 0 WHILE ((i < @maxi) && (!found)) i = i + 1 os = s complex sd = 1 int j = 1 bool far = false WHILE ((j < @numi) && !far) j = j + 1 sd = 2*sd*s + 1 s = sqr(s) + os IF (@trap) IF (|sd| > q) far = true ENDIF ENDIF ENDWHILE ; s is now m(s) and sd is m'(s), with m(x) = (x^2+x)^2+x... s = os - (s/sd) ; Apply Newton's method once found = (|s - os| < p) ENDWHILE IF (!found) s = 0 ENDIF transform: #pixel = #pixel + s default: title = "M-Set Component Centers" param numi caption = "Order of component" default = 6 min = 1 endparam param seed caption = "Seed" hint = "Determines which component is selected." default = (4,7) endparam param maxi caption = "Maximum iterations" default = 1000 min = 1 endparam param trap caption = "Trap points" default = false endparam } pgd_RiemannSphere { ; Riemann sphere visualizations and mapping projections ; Treats the Riemann sphere as an actual globe in 3-space, ; and visualizes the globe with map projections like ; Mercator, etc. ; ; Spherical metric: ; ds = 2|dw|/(1 + |w|^2) ; s = 2atan |w| ; The distance from 0 to w with the spherical metric is ; 2atan |w|. ; Given x, y, z = sqrt(1 - x^2 - y^2) on the unit ; sphere in R^3, the arc length from (0,0,1) to (x,y,z) ; on the sphere is the fraction s of a complete circle ; traversed (with a complete circle being s = 2pi). In ; this case, put r = sqrt(x^2 + y^2) so ; z = sqrt(1 - r^2) and we have z = cos(s), r = sin(s). ; Thus |w| = tan(s/2) = sin(s/2)/cos(s/2) and ; r = sin(s) = 2sin(s/2)cos(s/2). ; To send (x, y, z) to w is thus to send (x , y) to ; (x , y)|w|/r. |w|/r = sin(s/2)/2sin(s/2)cos^2(s/2) ; = 1/2cos^2(s/2) ; z = cos(s) = 2cos^2(s/2) - 1 and ; 2cos^2(s/2) = 1 + z, so ; |w|/r = 1/(1 + z) = 1/(1 + sqrt(1 - r^2). ; ; Thus the upper half of a sphere in R^3 is mapped to ; the unit disk by taking the (x, y) coordinates of the ; point projected from the sphere onto the unit disk ; and then scale by 1/(1 + sqrt(1 - r^2)). The lower ; half is mapped similarly, with the scale factor now ; 1/(1 - sqrt(1 - r^2)), which notably sends r = 0 (the ; sphere's south pole) to the point at infinity. ; ; To turn a projection of that sphere into w is our ; task here. We'll use a conformal projection, since ; conformal mappings best preserve the form of ; Mandelbrot and Julia fractals. Mercator is a good ; choice. The x-axis is simply equal to longitude, so ; using (u, v) as the map coordinates, u = atan2(x + iy) ; = atan2(w); the other axis (v) is modified latitide. ; The latitude l = pi/2 - s. In the southern hemisphere ; v = -ln(tan(pi/4 - l/2)) = -ln(tan(s/2)). To invert ; these, we use (x + iy) = re^(iu) so we simply need to ; determine r or even |w| in terms of v. Since ; |w| = tan(s/2) this is easy: v = -ln|w| and therefore ; w = e^(iu)/e^v. If v = 0 this lies on the unit ; circle. Increasing v gives increasing |w|. In the ; northern hemisphere, v = ln(tan(pi/4 + l/2)) which is ; ln(tan(pi/2 - s)), but our tan(x - y) trig identity ; is useless for evaluating this because of all the ; tan(pi/2)s that come up. But the mapping w -> 1/w is ; an isometry for the spherical metric, so replacing w ; with 1/w replaces |w| with 1/|w| and exchanges ; hemispheres as though the sphere were rotated about ; the line through (1,0,0) and (-1,0,0). We can then ; use the southern hemisphere formula -- s is now the ; distance from the north pole for 1/w, which is (as ; w -> 1/w is an isometry) the distance from the south ; pole for w. Swapping hemispheres negates our ; v-coordinate, so we get |1/w| = 1/e^-v, or |w| = 1/e^v ; meaning that our equation w = e^(iu)/e^v is in fact ; valid for both hemispheres. ; Collapsing the exponentials gives w = e^(-v + iu), ; which means our Mercator projection ends up simply ; being (u + iv) = iln(w), the complex logarithm rotated ; 90 degrees. ; ; The Lambert conformal conic projection of a sphere ; seems interesting at first, as well, but when the sphere's ; north pole is used, it turns out to simply turn the ; sphere back into the complex plane -- u + iv = w. ; Using it at the south pole gives u + iv = 1/w. ; ; Maps that show the whole sphere in a finite region ; are also desirable. These cannot be conformal, since ; conformal mappings are complex analytic so a conformal ; map from w to u + iv that omits more than 2 points ; of the (u,v) plane is constant, i.e. one point u + iv ; corresponds to all w values! To get a feel for the ; Riemann sphere's geometry we will probably want ; orthographic projection and at least one equal-area ; projection. Mollweide, which maps a sphere to an ; ellipse, seems a good bet. ; ; The orthographic projection about the north pole is ; easy: u = x, v = y, and thus ; w = (u + iv)/(1 + sqrt(1 - r^2)) ; = (u + iv)/(1 + sqrt(1 - u^2 - v^2)) ; To get the south pole, change w above to 1/w. ; ; Mollweide: ; If we set t = sin^-1(v/sqrt(2)), we get ; l = sin^-1[(2t + sin(2t))/pi] and ; atan2(w) = pi*u/(2sqrt(2)cos(t)) ; We have |w| = tan(s/2) ; = tan(pi/4 - l/2) ; = (tan(pi/4) - tan(l/2))/(1 + tan(pi/4)tan(l/2)) ; = (1 - tan(l/2))/(1 + tan(l/2)) ; = (1 - sin(l/2)/cos(l/2))/(1 + sin(l/2)/cos(l/2)) ; = (cos(l/2) - sin(l/2))/(cos(l/2) + sin(l/2)) ; = (2cos^2(l/2) - 2sin(l/2)cos(l/2))/(2cos^2(l/2) + 2sin(l/2)cos(l/2)) ; = (cos^2(l/2) - 2cos(l/2)sin(l/2) + sin^2(l/2))/(cos^2(l/2) - sin^2(l/2)) ; = (1 - sin(l))/cos(l) ; = (1 - sin(l))/sqrt(1 - sin^2(l)) ; |w|^2 = (1 - sin(l))^2/(1 - sin^2(l)) ; = (1 - sin(l))(1 - sin(l))/(1 - sin(l))(1 + sin(l)) ; = (1 - sin(l))/(1 + sin(l)) ; = (pi - 2t - sin(2t))/(pi + 2t + sin(2t)) ; sin(2t) = 2sin(t)cos(t) = sqrt(2)cos(t)v ; cos^2(t) = 1 - sin^2(t) = 1 - v^2/2 so ; sin(2t) = sqrt(2)sqrt(1 - v^2/2)v ; = sqrt(2 - v^2)v ; |w|^2 = (pi - 2sin^-1(v/sqrt(2)) - sqrt(2 - v^2)v)/(pi + 2sin^-1(v/sqrt(2)) + sqrt(2 - v^2)v) ; atan2(w) = (pi*u)/2sqrt(2 - v^2) ; w = sqrt((pi - 2sin^-1(v/sqrt(2)) - sqrt(2 - v^2)v)/(pi + 2sin^-1(v/sqrt(2)) + sqrt(2 - v^2)v))e^(pi*ui/2sqrt(2 - v^2)) ; 2 - v^2 must not become negative so v^2 <= 2 and ; -sqrt(2) <= v <= sqrt(2). pi*ui/2sqrt(2 - v^2) is ; periodic, with -1 <= u/2sqrt(2 - v^2) <= 1 giving ; one repetition. The ellipse we are interested in ; thus satisfies u^2 <= 4(2 - v^2), or u^2 + 4v^2 <= 8. ; u^2 + 4v^2 > 8 can thus be made #solid. ; The extent of the ellipse ranges from -sqrt(2) to ; sqrt(2) in the v direction; the maximum extent in the ; u direction occurs for v = 0, and is then from ; -2sqrt(2) 2sqrt(2). It seems sensible to reduce all ; coordinates by sqrt(2), which gives a new equation: ; w = sqrt((pi - 2sin^-1(v) - 2sqrt(1 - v^2)v)/(pi + 2sin^-1(v) + 2sqrt(1 - v^2)v))e^(pi*ui/2sqrt(1 - v^2)) ; which gives an ellipse with -2 <= u <= 2, -1 <= v <= 1. ; u^2 + 4v^2 > 4 is the solid region outside the ; ellipse now. ; ; NOTE: The code below all reverses this so that the ; point at infinity is the north pole, not the south! transform: complex p = 0 IF (@type == 0) ; Mercator is just a rotated complex exp. IF (abs(real(#pixel)) > #pi/2) #solid = true ELSE p = exp((0,-2)*#pixel) ENDIF ELSEIF (@type == 1) IF (|#pixel - 1| <= 1) ; Northern hemisphere p = -(1 + sqrt(1 - |#pixel - 1|))/(#pixel - 1) ELSEIF (|#pixel + 1| <= 1) ; Southern hemisphere p = (#pixel + 1)/(1 + sqrt(1 - |#pixel + 1|)) ELSE #solid = true ENDIF ELSEIF (@type == 2) float u = -real(#pixel) float v = -imag(#pixel) IF (sqr(u) + 4*sqr(v) > 4) #solid = true ELSE p = sqrt((#pi - 2*asin(v) - 2*v*sqrt(1 - sqr(v)))/(#pi + 2*asin(v) + 2*v*sqrt(1 - sqr(v))))*exp(#pi*(0,1)*u/(2*sqrt(1 - sqr(v)))) ENDIF ENDIF p = p*@scale*exp((0,1)*pi*@zmer/180) IF (@turn) p = 1/(p + 1) - 1/2 ENDIF #pixel = p default: title = "Riemann Sphere Map Projections" param type caption = "Projection type" enum = "Mercator" "Orthographic" "Mollweide" default = 0 endparam param zmer caption = "Zero Meridian" default = 0.0 hint = "Change this to change the center meridian in \ Mercator and Mollweide projections and the \ horizontal meridian in orthographic. Angle is \ in degrees." endparam param turn caption = "Turn sphere" default = false hint = "Turning this on rotates the sphere to put 1 \ and -1 at the poles instead of 0 and infinity." endparam param scale caption = "Scale fractal" default = 1.0 hint = "Turning this on scales the fractal being \ mapped to the sphere." endparam } pgd_MandelSpeedup { ; Obliterates the cardioid and largest disk of the ; standard M-set, making seahorse valley and other ; areas compute faster. ; Don't use in layers that are inside-coloring these ; areas. ; Example: parameter set Elephant in pdf.upr uses this ; transformation. Switch to guessing and disable the ; MandelSpeedup transformation and it takes a 1.5GHz ; Athlon 30 seconds to generate. Enable MandelSpeedup and ; switch back to 1-pass and the fractal comes up in 10 ; seconds on the same hardware AND IS MORE ACCURATE. ; Three-fold speedups are possible for Mandelbrot zooms on ; the boundary of the cardioid or largest disk. transform: IF (|1 + sqrt(1 - 4*#pixel)| <= 1) #solid = true ELSEIF (|1 - sqrt(1 - 4*#pixel)| <= 1) #solid = true ELSEIF (|#pixel + 1| <= 0.0625) #solid = true ELSE #pixel = #pixel ENDIF default: title = "Mandelbrot Masking Speedup" }