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"
}