comment { dmj-pub.ucl 2.1 Coloring methods for Ultra Fractal 2 by Damien M. Jones April 2, 2000 For more information about this formula collection, please visit its home page: http://www.fractalus.com/ultrafractal/dmj-pub-uf.htm Many of these coloring algorithms are general-purpose tools which work with a variety of fractals. A few are more specialized. Huge portions of this compilation build upon the techniques explored by others. } dmj-Bifurcation { ; ; Intended to color bifurcation plots. ; It assumes the fractal formula will have placed the count of ; orbit points landing within the pixel area into real(z). ; It then performs some simple scaling of the result. ; final: #index = real(#z)/(#maxit-@filter) default: title = "Bifurcation" helpfile = "dmj-pub\dmj-pub-uf-bifurcationc.htm" param filter caption = "Filter Iterations" default = 200 hint = "Specifies the number of iterations to skip before \ plotting orbit points. (This should match what is on \ the Formula tab.)" endparam } dmj-CarlsonAtan { ; ; This is the Paul Carlson ATAN method. ; Basically, it is atan(z[n] - z[n-1]). ; On most fractals it will give slightly different ; results than the FractInt ATAN, but on some ; (especially Nova and other convergent fractals) ; it yields substantially different images. ; init: complex zold = (0,0) float d = 0.0 loop: zold = #z ; save current orbit value final: zold = #z - zold ; take final z and subtract previous z d = atan2(zold) ; retrieve the angle of that IF (d < 0) ; it's negative d = d + #pi * 2 ; make it positive ENDIF #index = d / (#pi * 2) ; scale so full angle range is 0..1 default: title = "Atan (Paul Carlson)" helpfile = "dmj-pub\dmj-pub-uf-atan.htm" } dmj-CarlsonAtan2 { ; ; This is a variation on the Paul Carlson ATAN ; method. Basically, it is atan(z[n] - z[n-2]). ; init: complex zold = (0,0) complex zold2 = (0,0) float d = 0.0 loop: zold2 = zold ; save previous iteration zold = #z ; save current iteration final: zold = #z - zold2 ; take final z and subtract 2nd previous z d = atan2(zold) ; determine angle of zold IF (d < 0) ; it's negative d = d + #pi * 2 ; make it positive ENDIF #index = d / (#pi * 2) ; scale so full angle range is 0..1 default: title = "Atan2 (Paul Carlson variant)" helpfile = "dmj-pub\dmj-pub-uf-atan.htm" } dmj-Cilia(OUTSIDE) { ; ; This coloring method is an implementation of the ; "Cilia" coloring type from the KPT5 FraxPlorer. ; The technique was developed by Ben Weiss of ; MetaCreations; I've added several new modes and ; parameters to it. Many thanks to Ben for confirming ; my guesses as to how this algorithm works. ; ; Because this method relies on smoothed iteration ; values, it will work best with Mandelbrot and ; derived types. It also DOES NOT WORK YET with ; fractals with exponents other than 2! Feel free ; to experiment regardless, but the results will ; not be seamless. ; init: complex il = 1/log(2.0) ; Inverse log (power). float lp = log(log(@bailout)) ; log(log bailout). complex z1 = (0,0) ; old iterations complex z2 = (0,0) complex z3 = (0,0) complex z4 = (0,0) loop: IF (@olditer > 0) ; need old iterations z4 = z3 ; save them z3 = z2 z2 = z1 z1 = #z ENDIF final: IF (@olditer == 0) ; select the iteration we want z4 = #z ELSEIF (@olditer == 1) z4 = z1 ELSEIF (@olditer == 2) z4 = z2 ELSEIF (@olditer == 3) z4 = z3 ENDIF ; determine potential and decomposition ; this gives us a "coordinate" in d and f, ranging from 0 to 1 float f = real(il*lp - il*log(log(cabs(#z)))) ; fractional iteration float d = atan2(z4) ; angle of z IF (d < 0) ; it's negative d = d + #pi * 2 ; make it positive ENDIF d = d / (#pi * 2) ; scale so full angle range is 0..1 ; apply skew d = d + @skew*f d = d - floor(d) ; f = f - floor(f) ; convert decomp and pot to tileable value d ; this is where the patterns are made IF (@ciliamode == 0) ; sawtooth IF (d < 0.5) d = d * (f+1) ELSE d = (d-0.5) * (f+1) - f/2 + 0.5 ENDIF d = abs(1-d*2) ELSEIF (@ciliamode == 1) ; interpolated float d1 = abs(1-d*2) float d2 = abs(1-(d*real(2.0)-floor(d*real(2.0)))*2) d = d1 + f*(d2-d1) ELSEIF (@ciliamode == 2) ; sine waves float d1 = cos(d*#pi*2)*0.5 + 0.5 float d2 = cos(d*#pi*2*real(2.0))*0.5 + 0.5 d = d1 + f*(d2-d1) ELSEIF (@ciliamode == 3) ; cells d = ((0.25-sqr(d-0.5)) * (0.25-sqr(f-0.5))) * 16 ELSEIF (@ciliamode == 4) ; bricks float d1 = 1 - abs(d-0.5)*2 float d2 = 1 - abs(f-0.5)*2 IF (d1 < d2) d = d1 ELSE d = d2 ENDIF ELSEIF (@ciliamode == 5) ; diamonds d = 0.5 - abs(d - 0.5) IF (f < 2/3) d = d + f*0.75 d = abs(1-d*2) ELSE IF (d < f*0.75-0.5) d = f*3-2-d*4 ELSE d = 1-(0.5-d)/(1-f*0.75) ENDIF ENDIF ENDIF #index = d * 0.99 ; reduce it just a bit to prevent ; the colors from wrapping default: title = "Cilia" helpfile = "dmj-pub\dmj-pub-uf-cilia.htm" ; param power ; caption = "Exponent" ; default = (2,0) ; hint = "This should be set to match the exponent of the \ ; formula you are using. For Mandelbrot, this is 2." ; endparam param bailout caption = "Bailout value" default = 65536.0 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab." endparam param ciliamode caption = "Cilia Mode" default = 0 enum = "sawtooth" "interpolated" "sine waves" "cells" "bricks" "diamonds" hint = "This is the periodic function used to generate the \ cilia pieces. Each has a different 'flavor'." endparam param skew caption = "Cilia Skew" default = 0.0 hint = "This skews the pattern. Use multiples of 0.5 for \ complete continuity." endparam param olditer caption = "Use Older Iterations" default = 0 enum = "no" "1 older" "2 older" "3 older" "4 older" hint = "Uses older iterations for decomposition angle; this \ has the effect of spreading out the pattern. Best \ used to compensate for higher bailouts." endparam } dmj-Curvature { ; ; This formula averages curvature over the course ; of all iterations by averaging the difference in ; absolute value of angles between orbit steps. The ; technique was suggested by Prof. Javier Barrallo. ; ; Note that the results are often very similar to ; Kerry Mitchell's Triangle Inequality Average method, ; but are slightly more angular. ; init: complex zold = (0,0) complex zold2 = (0,0) float a = 0.0 float a2 = 0.0 int i = 0 loop: a2 = a IF (i >= 2 && @aflavor != 2) ; zold and zold2 are valid a = a + abs(atan2((#z-zold)/(zold-zold2))) ; update average ENDIF i = i + 1 ; count the iteration zold2 = zold ; save the orbit values zold = #z final: IF (@aflavor == 0) float il = 1/log(@power) float lp = log(log(@bailout)/2.0) float f = il*lp - il*log(log(cabs(#z))) a = a / i a2 = a2 / (i-1) #index = (a2 + (a-a2) * (f+1)) / #pi ELSEIF (@aflavor == 1) #index = (a/i) / #pi + 1 ELSE a = atan2((#z-zold)/(zold-zold2)) ; update average #index = a / #pi + 1 ENDIF default: title = "Curvature Average" helpfile = "dmj-pub\dmj-pub-uf-ca.htm" param aflavor caption = "Average Flavor" default = 0 enum = "smooth" "average" "last" endparam param power caption = "Exponent" default = 2.0 hint = "This should be set to match the exponent of the \ formula you are using. For Mandelbrot, this is 2." endparam param bailout caption = "Bailout" default = 1e20 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab. Use a very high bailout!" endparam } dmj-Decomp { ; ; This coloring method simulates the decomp ; option in FractInt. Note that this will ; give smooth coloring; this may make matching ; some uses of decomp difficult. ; ; Note also that unlike the Compatibility ; formula that comes with Ultra Fractal, this ; method can be used for inside coloring, too. ; final: float d = atan2(#z) ; get angle of z IF (d < 0) ; it's negative d = d + #pi * 2 ; make it positive ENDIF #index = d / (#pi * 2) * @portion / 256 ; scale so full angle range is 0..1 default: title = "Final Decomposition" helpfile = "dmj-pub\dmj-pub-uf-decomp.htm" param portion caption = "Palette Portion" default = 256.0 hint = "Indicates the portion of the palette to use. \ 256 means use the whole palette, 128 means use \ half of it." endparam } dmj-Distance(OUTSIDE) { ; ; This is the distance-estimator coloring method ; for Mandelbrot and other z^n formula types ; (Phoenix, Julia). Results on other types may ; be unpredictable, but might be interesting. ; ; To match the coloring of FractInt's distance ; estimator method, use a Sqr transfer function ; and a coloring style of "distance". This is not, ; however, a perfect match. ; ; The "real" and "imaginary" coloring styles are ; not simply the real and imaginary components ; of the distance/angle calculations. These were ; exceptionally dull. I've mangled these modes to ; try to produce something more interesting, but ; generally these were an experiment that didn't ; work at all like I expected them to. ; init: complex dz = (0,0) loop: dz = @power * #z^(@power-1) * dz + 1 final: IF (@style == 0) ; distance #index = (@power*log(cabs(#z)) * cabs(#z) / cabs(dz))^(1/@power) ELSEIF (@style == 1) ; angle #index = atan2(#z/dz)/#pi + 1 ELSEIF (@style == 2) ; real #index = (abs(real(#z/dz)))^(1/exp(@power)) * log(#numiter) ELSEIF (@style == 3) ; imaginary #index = (abs(imag(#z/dz)))^(1/exp(@power)) * log(#numiter) ENDIF default: title = "Distance Estimator (Mandelbrot)" helpfile = "dmj-pub\dmj-pub-uf-dem.htm" param style caption = "Coloring Style" default = 0 enum = "distance" "angle" "real" "imaginary" hint = "Selects the value to use from the distance estimator \ function." endparam param power caption = "Exponent" default = 2.0 hint = "This should be set to match the exponent of the \ formula you are using. For Mandelbrot, this is 2." endparam } dmj-fBm { ; ; Fractional Brownian Motion ; This is an adaptation of Intel's real-time MMX fBm ; algorithm, which is derived from Ken Musgrave's fBm ; algorithm, which is built on Perlin's noise basis ; function. It provides bumpy coloring that can be ; rendered at high resolution. The algorithm is fractal. ; ; You will find that this coloring has sharp breaks at ; iteration boundaries. That's because it was created to ; color "flat" areas, not those twisted by fractal ; formulas. Interesting effects can be obtained if you can ; arrange for all pixels to use the same number of fractal ; formula iterations. ; ; I'm not going to try to provide a complete description ; of how this algorithm works, because it has several ; layers. I was able to find a fairly decent description ; on Intel's site (which is where the outline of this ; code comes from) but it is written for assembly language ; programmers and not mathematicians (let alone laypeople) ; and I can't recommend it. I'll write better docs when ; I do the full help file package. ; final: complex r = (0,1) ^ (@angle / 90.0) complex r2 = (0,1) ^ (@anglestep / 90.0) complex p = #z * @scale * r + @offset float sum = 0.0 float freq = 1.0 int i = @octaves WHILE (i > 0) ; determine integer coordinate for corners of square ; surrounding p float bx0 = floor(real(p)) % 256 float by0 = floor(imag(p)) % 256 IF (bx0 < 0) bx0 = bx0 + 256 ENDIF IF (by0 < 0) by0 = by0 + 256 ENDIF float bx1 = (bx0 + 1) % 256 float by1 = (by0 + 1) % 256 float rx0 = real(p) - floor(real(p)) float ry0 = imag(p) - floor(imag(p)) float rx1 = rx0 - 1 float ry1 = ry0 - 1 ; create a "random" index for each corner ; (this is where Intel's version differs from Perlin's; ; I used Intel's version because it doesn't require a ; pre-computed random table, which is difficult to manage ; in UF.) float b00 = (bx0^@power % 65536 + by0)^@power % 65536 float b10 = (bx1^@power % 65536 + by0)^@power % 65536 float b01 = (bx0^@power % 65536 + by1)^@power % 65536 float b11 = (bx1^@power % 65536 + by1)^@power % 65536 ; produce a "random" vector for each corner float g_b00_0 = (b00)^@power*0.25 % 512 - 256 float g_b10_0 = (b10)^@power*0.25 % 512 - 256 float g_b01_0 = (b01)^@power*0.25 % 512 - 256 float g_b11_0 = (b11)^@power*0.25 % 512 - 256 float g_b00_1 = (b00+1)^@power*0.25 % 512 - 256 float g_b10_1 = (b10+1)^@power*0.25 % 512 - 256 float g_b01_1 = (b01+1)^@power*0.25 % 512 - 256 float g_b11_1 = (b11+1)^@power*0.25 % 512 - 256 ; normalize each vector float d = 0.0; d = 1 / sqrt(sqr(g_b00_0) + sqr(g_b00_1)) g_b00_0 = g_b00_0 * d g_b00_1 = g_b00_1 * d d = 1 / sqrt(sqr(g_b10_0) + sqr(g_b10_1)) g_b10_0 = g_b10_0 * d g_b10_1 = g_b10_1 * d d = 1 / sqrt(sqr(g_b01_0) + sqr(g_b01_1)) g_b01_0 = g_b01_0 * d g_b01_1 = g_b01_1 * d d = 1 / sqrt(sqr(g_b11_0) + sqr(g_b11_1)) g_b11_0 = g_b11_0 * d g_b11_1 = g_b11_1 * d ; produce colors for each corner float u1 = rx0 * g_b00_0 + ry0 * g_b00_1 float v1 = rx1 * g_b10_0 + ry0 * g_b10_1 float u2 = rx0 * g_b01_0 + ry1 * g_b01_1 float v2 = rx1 * g_b11_0 + ry1 * g_b11_1 ; interpolate between corners using ; bilinear filtering float sx = sqr(rx0) * (3 - rx0*2) float sy = sqr(ry0) * (3 - ry0*2) float a = u1 + sx*(v1-u1) float b = u2 + sx*(v2-u2) sum = sum + (a + sy*(b-a))*freq freq = freq * @step p = p * r2 / @step i = i - 1 ENDWHILE #index = (sum + 1) * 0.5 default: title = "Fractional Brownian Motion" helpfile = "dmj-pub\dmj-pub-uf-fbm.htm" param offset caption = "Offset" default = (0,0) hint = "This is the offset of the pattern. You can use this to shift \ the pattern around on the complex plane." endparam param scale caption = "Scale" default = 1.0 hint = "This is the overall scale of the noise." endparam param angle caption = "Rotation" default = 0.0 hint = "This is the angle, in degrees, of the noise." endparam param step caption = "Scale Step" default = 0.5 hint = "This is the step in scale between noise iterations." endparam param anglestep caption = "Rotation Step" default = 37.0 hint = "This is the angle, in degrees, to rotate between noise \ iterations." endparam param octaves caption = "Octaves" default = 7 min = 1 hint = "This is the number of iterations of the noise formula." endparam param power caption = "Exponent" default = 2.0 hint = "This is the exponent used to scramble numbers." endparam } dmj-fBm2 { ; ; Fractional Brownian Motion II ; This is an adaptation of Intel's real-time MMX fBm ; algorithm, which is derived from Ken Musgrave's fBm ; algorithm, which is built on Perlin's noise basis ; function. It provides bumpy coloring that can be ; rendered at high resolution. ; ; This variant uses the fBm texture as a cumulative ; effect at each iteration. As such, you can usually ; get away with using only a few octaves at each ; iteration. ; init: complex r = (0,1) ^ (@angle / 90.0) complex r2 = (0,1) ^ (@anglestep / 90.0) float sum = 0.0 loop: float fade2 = 1-cabs(#z)/@fade IF (fade2 < 0) fade2 = 0 ELSE fade2 = sqr(fade2) ENDIF complex p = #z * @scale * r + @offset float freq = 1.0 int i = @octaves WHILE (i > 0) ; determine integer coordinate for corners of square ; surrounding p float bx0 = floor(real(p)) % 256 float by0 = floor(imag(p)) % 256 IF (bx0 < 0) bx0 = bx0 + 256 ENDIF IF (by0 < 0) by0 = by0 + 256 ENDIF float bx1 = (bx0 + 1) % 256 float by1 = (by0 + 1) % 256 float rx0 = real(p) - floor(real(p)) float ry0 = imag(p) - floor(imag(p)) float rx1 = rx0 - 1 float ry1 = ry0 - 1 ; create a "random" index for each corner ; (this is where Intel's version differs from Perlin's; ; I used Intel's version because it doesn't require a ; pre-computed random table, which is difficult to manage ; in UF.) float b00 = (bx0^@power % 65536 + by0)^@power % 65536 float b10 = (bx1^@power % 65536 + by0)^@power % 65536 float b01 = (bx0^@power % 65536 + by1)^@power % 65536 float b11 = (bx1^@power % 65536 + by1)^@power % 65536 ; produce a "random" vector for each corner float g_b00_0 = (b00)^@power*0.25 % 512 - 256 float g_b10_0 = (b10)^@power*0.25 % 512 - 256 float g_b01_0 = (b01)^@power*0.25 % 512 - 256 float g_b11_0 = (b11)^@power*0.25 % 512 - 256 float g_b00_1 = (b00+1)^@power*0.25 % 512 - 256 float g_b10_1 = (b10+1)^@power*0.25 % 512 - 256 float g_b01_1 = (b01+1)^@power*0.25 % 512 - 256 float g_b11_1 = (b11+1)^@power*0.25 % 512 - 256 ; normalize each vector float d = 0.0; d = 1 / sqrt(sqr(g_b00_0) + sqr(g_b00_1)) g_b00_0 = g_b00_0 * d g_b00_1 = g_b00_1 * d d = 1 / sqrt(sqr(g_b10_0) + sqr(g_b10_1)) g_b10_0 = g_b10_0 * d g_b10_1 = g_b10_1 * d d = 1 / sqrt(sqr(g_b01_0) + sqr(g_b01_1)) g_b01_0 = g_b01_0 * d g_b01_1 = g_b01_1 * d d = 1 / sqrt(sqr(g_b11_0) + sqr(g_b11_1)) g_b11_0 = g_b11_0 * d g_b11_1 = g_b11_1 * d ; produce colors for each corner float u1 = rx0 * g_b00_0 + ry0 * g_b00_1 float v1 = rx1 * g_b10_0 + ry0 * g_b10_1 float u2 = rx0 * g_b01_0 + ry1 * g_b01_1 float v2 = rx1 * g_b11_0 + ry1 * g_b11_1 ; interpolate between corners using ; bilinear filtering float sx = sqr(rx0) * (3 - rx0*2) float sy = sqr(ry0) * (3 - ry0*2) float a = u1 + sx*(v1-u1) float b = u2 + sx*(v2-u2) sum = sum + (a + sy*(b-a))*freq*fade2 freq = freq * @step p = p * r2 / @step i = i - 1 ENDWHILE final: #index = (sum + @coloroffset) * 0.5 default: title = "Fractional Brownian Motion II" helpfile = "dmj-pub\dmj-pub-uf-fbm.htm" param fade caption = "Fade Amount" default = 16.0 hint = "Indicates the degree to which the texture fades with distance \ from the origin. Larger values give less fading." endparam param coloroffset caption = "Coloring Offset" default = 1.0 hint = "Sets the coloring offset from zero. Increase this if you see \ large areas of solid color." endparam param offset caption = "Offset" default = (0,0) hint = "This is the offset of the pattern. You can use this to shift \ the pattern around on the complex plane." endparam param scale caption = "Scale" default = 1.0 hint = "This is the overall scale of the noise." endparam param angle caption = "Rotation" default = 0.0 hint = "This is the angle, in degrees, of the noise." endparam param step caption = "Scale Step" default = 0.5 hint = "This is the step in scale between noise iterations." endparam param anglestep caption = "Rotation Step" default = 37.0 hint = "This is the angle, in degrees, to rotate between noise \ iterations." endparam param octaves caption = "Octaves" default = 1 min = 1 hint = "This is the number of iterations of the noise formula." endparam param power caption = "Exponent" default = 2.0 hint = "This is the exponent used to scramble numbers." endparam } dmj-FiniteAttractor(INSIDE) { ; ; This coloring method colors based on how quickly the ; orbit falls into a periodic orbit. Since the coloring ; is based on iteration count, it is discrete. ; ; NOTE: This formula uses oldz(). It is thus not suitable ; for very high iteration limits (over 100,000). ; final: int i = 0 bool found = false WHILE (i < #maxit && !found) IF (|oldz(i)-#z| < @threshold) found = true ELSE i = i + 1 ENDIF ENDWHILE #index = i * 0.05 default: title = "Finite Attractor" helpfile = "dmj-pub\dmj-pub-uf-fa.htm" param threshold caption = "Threshold" default = 1e-5 hint = "Specifies how close the orbit must approach the \ finite attractor before being considered to have \ reached it. This is similar to a bailout for \ convergent types." endparam } dmj-Light { ; ; This coloring method expects a surface normal in z ; and performs lighting for that point. This is best ; used with the "Slope" family of formulas, which ; provide exactly the information expected by this ; coloring method. ; ; Surface normals are 3D vectors with a length of 1. ; So while only two components of this vector can be ; passed through a complex number, the third can be ; derived by sqrt(1-sqr(real)-sqr(imag)). The full ; vector gives the orientation of a 3D surface, which ; is needed to calculate lighting. Note that this is ; an *extremely* simple lighting model; to do a more ; accurate lighting model would require more parameters ; from the fractal formula (such as a position in 3D ; space). ; final: float vz = -sqrt(1-|#z|) ; extract implied portion of normal float d2r = #pi/180 ; degrees to radians conversion factor ; create vector for light direction float lx = cos(@angle*d2r) * cos(@elevation*d2r) float ly = sin(@angle*d2r) * cos(@elevation*d2r) float lz = -sin(@elevation*d2r) ; compute cosine of angle between these vectors ; (this is the amount of lighting on the surface) float l = lx*real(#z) + ly*imag(#z) + lz*vz IF (l < @ambient) ; light is below the ambient level l = @ambient ; set it to the ambient level ENDIF IF (@ambient < 0) ; the ambient level is negative l = l + 1 ; offset to prevent clipping at 0 ENDIF #index = l*0.99 ; reduce it just a bit to prevent ; the colors from wrapping default: title = "Lighting" helpfile = "dmj-pub\dmj-pub-uf-lighting.htm" param @angle caption = "Light Rotation" default = 90.0 hint = "Gives the rotation of the light source, in degrees." endparam param @elevation caption = "Light Elevation" default = 30.0 hint = "Gives the elevation of the light source, in degrees." endparam param @ambient caption = "Ambient Light" default = 0.0 min = -1.0 max = 1.0 hint = "Specifies the level of ambient light. Use -1.0 to \ color all surfaces." endparam } dmj-Logistics { ; ; This coloring algorithm should be used with the ; Markus-Lyapunov fractal type. It expects the sum ; of the logarithms at each iteration to be passed ; in real(z). It then performs some simple scaling. ; final: #index = -real(#z)/#numiter default: title = "Markus-Lyapunov" helpfile = "dmj-pub\dmj-pub-uf-mlc.htm" } dmj-Lyapunov { ; ; This algorithm computes the Lyapunov exponent ; for Mandelbrot types. This exponent is usually ; negative for divergent orbits (outside points), ; and positive for convergent orbits (inside ; points). ; ; This might have interesting results when used ; with other fractal types, although the results ; would not be mathematically accurate. ; ; Optimizations suggested by Charles Vassallo. ; init: float oldsum = 0 float sum = 1 float v = 0 float il = 1/log(real(@power)) float lp = log(log(@bailout)/2.0) float f = 0.0 loop: IF (@trackvariable == 0) ; |z| v = cabs(#z) ELSEIF (@trackvariable == 1) ; real(z) v = real(#z) ELSEIF (@trackvariable == 2) ; imag(z) v = imag(#z) ELSEIF (@trackvariable == 3) ; real(z)/imag(z) v = real(#z)/imag(#z) ELSEIF (@trackvariable == 4) ; imag(z)/real(z) v = imag(#z)/real(#z) ELSEIF (@trackvariable == 5) ; arg(z) v = atan2(#z) ELSEIF (@trackvariable == 6) ; 1/real(z) v = 1.0/real(#z) ELSEIF (@trackvariable == 7) ; 1/imag(z) v = 1.0/imag(#z) ENDIF oldsum = sum ; sum = sum + log(abs(2*v)) ; sum the Lyapunov exponent (slow method) sum = sum * (abs(2*v)) ; sum the Lyapunov exponent final: oldsum = log(oldsum) sum = log(sum) IF (@negative == 1) sum = -sum/#numiter oldsum = -oldsum/(#numiter-1) ELSEIF (@negative == 2) sum = abs(sum/#numiter) oldsum = abs(oldsum/(#numiter-1)) ELSE sum = sum/#numiter oldsum = oldsum/(#numiter-1) ENDIF IF (@smooth) f = il*lp - il*log(log(cabs(#z))) #index = oldsum + (sum-oldsum) * (f+1) ELSE #index = sum ENDIF default: title = "Lyapunov" helpfile = "dmj-pub\dmj-pub-uf-lyapunov.htm" param trackvariable caption = "Variable to Track" default = 0 enum = "magnitude of z" "real part of z" "imaginary part of z" \ "real / imag" "imag / real" "angle of z" "1 / real(z)" "1 / imag(z)" hint = "Indicates which variable to measure the Lyapunov exponent for." endparam param negative caption = "Sign" default = 2 enum = "positive" "negative" "absolute value" hint = "Affects the sign of the exponent. 'Negative' and 'absolute \ value' are useful for inside coloring." endparam param power caption = "Exponent" default = 2.0 hint = "This should be set to match the exponent of the \ formula you are using. For Mandelbrot, this is 2." endparam param bailout caption = "Bailout" default = 1e20 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab. Use a very high bailout!" endparam param smooth caption = "Smooth Coloring" default = false hint = "If set, results will be 'smoothed' to hide iteration bands." endparam } dmj-Sampling { ; ; This coloring algorithm samples the iterating fractal for ; the current point. It works best when the fractal formula ; iterates the same sequence for every pixel or large group ; of pixels, as the bifurcation fractals do. But, it will ; actually produce results for all fractal types. ; ; NOTE: Box sampling is not completely accurate if your image ; is rotated. This is the crudest form of sampling anyway; ; if you're interested in accuracy, you probably want point ; sampling. Box sampling is only in here as a rough-cut quick ; render for some other formulas I was working on. ; ; NOTE: Box sampling is RESOLUTION-DEPENDENT. This means if ; you zoom such an image, your samples will be more sparsely ; spaced, and you will need to increase the maximum number ; of iterations to compensate. Other sampling modes are NOT ; dependent on image size. ; ; HINT: Don't use box sampling! :-) ; init: int counter = 0 float iradius = 1/@radius float sradius = sqr(@radius) float xmin = real(#pixel) float xmax = real(#pixel) + 4/#magn/#height float ymin = imag(#pixel) float ymax = imag(#pixel) + 4/#magn/#height float sum = 1e20 IF (@select == 0 || @countall == false) sum = 0.0 ENDIF float d = 0.0 complex oldz = (0,0) complex c = (0,0) IF (@mode >= 2) counter = -1 ; these modes have a built-in filter ENDIF loop: IF (counter >= @filter) ; passed iteration filter, count this IF (@mode == 0) ; box sampling IF (real(#z) >= xmin && real(#z) < xmax && \ imag(#z) >= ymin && imag(#z) < ymax) d = 1.0 ; point is inside box, counts as 1 ELSE d = 0.0 ; point is outside box, counts as 0 ENDIF ELSEIF (@mode == 1) ; point sampling d = |#z-#pixel| ; get distance to point IF (@countall) ; counting all pixels d = sqrt(d) ; use simple linear distance ELSEIF (d < sradius) ; within radius 1 d = 1 - sqrt(d)*iradius ; ranges linearly from 1 to 0 ELSE ; not within radius 1 d = 0 ; does not count ENDIF ELSEIF (@mode == 2) ; line sampling ; note: if you use sum with this mode, the endpoints ; of the line segments will count twice. c = (#pixel-oldz)*conj(#z-oldz)/(|#z-oldz|) ; get distance along line ; from oldz to #z d = real(c) ; normalized distance along line IF (d < 0) ; past oldz end of line d = |oldz-#pixel| ; get distance to point IF (@countall) ; counting all pixels d = sqrt(d) ; use simple linear distance ELSEIF (d < sradius) ; within radius 1 d = 1 - sqrt(d)*iradius ; ranges linearly from 1 to 0 ELSE ; not within radius 1 d = 0 ; does not count ENDIF ELSEIF (d > 1) ; past #z end of line d = |#z-#pixel| ; get distance to point IF (@countall) ; counting all pixels d = sqrt(d) ; use simple linear distance ELSEIF (d < sradius) ; within radius 1 d = 1 - sqrt(d)*iradius ; ranges linearly from 1 to 0 ELSE ; not within radius 1 d = 0 ; does not count ENDIF ELSE ; within line segment boundaries d = abs(imag(c)) * cabs(#z-oldz) ; linear distance from line segment IF (@countall == false) ; need to fix distances IF (d < @radius) ; distance is within range d = 1 - d*iradius ; scale/adjust it ELSE ; distance is out of range d = 0 ; does not count ENDIF ENDIF ENDIF ENDIF ; apply sample transfer function ; a function is not used because these are floats ; and not all functions apply to floats d = d * @prescale IF (@sampfunc == 1) ; log d = log(d) ELSEIF (@sampfunc == 2) ; sqrt d = sqrt(d) ELSEIF (@sampfunc == 3) ; cuberoot d = (d)^(1/3) ELSEIF (@sampfunc == 4) ; exp d = exp(-d) ; NOTE this is different! ELSEIF (@sampfunc == 5) ; sqr d = sqr(d) ELSEIF (@sampfunc == 6) ; cube d = (d)^3 ELSEIF (@sampfunc == 7) ; sin d = sin(d) ELSEIF (@sampfunc == 8) ; cos d = cos(d) ELSEIF (@sampfunc == 9) ; tan d = tan(d) ELSEIF (@sampfunc == 10) ; sinc (modified form) IF (d == 0) d = 1 ELSE d = sin(#pi*d)/d ENDIF ENDIF d = d * @postscale IF (@select == 0) sum = sum + d ELSEIF (@select == 1) IF (d > sum) sum = d ENDIF ENDIF ENDIF counter = counter + 1 IF (@mode >= 2) oldz = #z ENDIF final: IF (@mode == 0) #index = sum/(#maxit-@filter) ELSEIF (@mode >= 1) #index = sum ENDIF default: title = "Sampling" helpfile = "dmj-pub\dmj-pub-uf-sampling.htm" param mode caption = "Sampling Mode" enum = "box" "points" "lines" default = 2 hint = "Sets the sampling mode. 'Box' counts orbits that fall \ within the pixel. 'Points' counts orbits that fall near \ the pixel. 'Lines' counts orbits that fall near lines \ connecting orbit points." endparam param select caption = "Sample Select" enum = "sum" "closest" default = 1 hint = "Sets which sample(s) are selected for plotting. 'Sum' works \ well with point sampling, 'closest' works well with lines." endparam param radius caption = "Sampling Radius" default = 0.1 hint = "Sets the sampling radius. A wider radius will result in \ more points being counted per pixel." endparam param sampfunc caption = "Sampling Transfer" enum = "linear" "log" "sqrt" "cuberoot" "exp" "sqr" "cube" \ "sin" "cos" "tan" "sinc" hint = "This function is applied to orbit distances at every \ iteration. It can be used to apply some special effects \ to sample shapes." endparam param prescale caption = "Distance Pre-scale" default = 1.0 hint = "This is a multiplier applied to the sample distance, \ before the Sampling Transfer function is applied to it." endparam param postscale caption = "Distance Post-scale" default = 1.0 hint = "This is a multiplier applied to the sample distance, \ after the Sampling Transfer function is applied to it." endparam param countall caption = "Ignore Radius" default = false hint = "If enabled, distances will not be restricted to a range \ from 0 to 1." endparam param filter caption = "Filter Iterations" default = 0 hint = "Specifies the number of iterations to skip before \ plotting orbit points. For most Mandelbrot types, you will \ need to set this to at least 1." endparam } dmj-Smooth(OUTSIDE) { ; ; This coloring method provides smooth iteration ; colors for Mandelbrot and other z^2 formula types ; (Phoenix, Julia). Results on other types may be ; unpredictable, but might be interesting. ; ; Thanks to F. Slijkerman for some tweaks. ; Thanks to Linas Vepstas for the math. ; init: complex il = 1/log(@power) ; Inverse log (power). float lp = log(log(@bailout)) ; log(log bailout). final: #index = 0.05 * real(#numiter + il*lp - il*log(log(cabs(#z)))) default: title = "Smoothed Iterations (Mandelbrot)" helpfile = "dmj-pub\dmj-pub-uf-smooth.htm" param power caption = "Exponent" default = (2,0) hint = "This should be set to match the exponent of the \ formula you are using. For Mandelbrot, this is 2." endparam param bailout caption = "Bailout value" default = 128.0 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab." endparam } dmj-SmoothConverge(OUTSIDE) { ; ; This coloring method provides smooth iteration ; colors for Magnet and other convergent/divergent ; formula types. Results on other types may be ; unpredictable, but might be interesting. ; ; The default values have been specifically tweaked ; for Magnet 1, although they are not perfect. ; init: complex il = 1/log(@power) ; Inverse log (power). float lp = log(log(@bailout)) ; log(log bailout). complex zold = (0,0) loop: zold = #z ; save old Z final: IF (|#z - zold| < 0.5) ; convergent bailout. IF (@converge) #index = 0.05 * real(#numiter + @fudge*0.5*(il*lp - il*log(log(cabs(#z-zold))))) ELSE #index = 0 ENDIF ELSE ; divergent bailout. IF (@diverge) #index = 0.05 * real(#numiter + il*lp - il*log(log(cabs(#z)))) ELSE #index = 0 ENDIF ENDIF default: title = "Smoothed Iterations (General)" helpfile = "dmj-pub\dmj-pub-uf-smooth.htm" param power caption = "Exponent" default = (1.414213562,0) hint = "This should be set to match the exponent of the \ formula you are using. For Magnet, this is sqrt(2)." endparam param bailout caption = "Bailout 1" default = 128 hint = "Divergent bailout value; this should be set to match the \ large bailout in the Formula tab." endparam param bailout2 caption = "Bailout 2" default = 0.00001 hint = "Convergent bailout value; this should be set to match the \ small bailout in the Formula tab." endparam param diverge caption = "Color Divergent" default = TRUE hint = "If set, points which escape to infinity will be \ colored." endparam param converge caption = "Color Convergent" default = TRUE hint = "If set, points which collapse to one value will be \ colored." endparam param fudge caption = "Fudge Factor" default = 1.0 hint = "Fractional iteration strength; use 15 for Nova." endparam } dmj-Smooth2 { ; ; This coloring method provides smooth iteration ; colors for all divergent fractal types. It was ; developed by Ron Barnett. It doesn't map ; precisely to iterations, but it's close. ; ; Note: You can achieve exactly this coloring by ; using the Orbit Traps coloring, with a merge ; mode setting of "funky average 4" and a distance ; transfer function of "sqrt". However, this ; formula is probably easier to add. ; init: float sum = 0.0 loop: sum = sum + exp(-cabs(#z)) final: #index = sum default: title = "Exponential Smoothing (Divergent)" helpfile = "dmj-pub\dmj-pub-uf-es.htm" } dmj-Smooth2Converge { ; ; This coloring method provides smooth iteration ; colors for all convergent fractal types. It was ; developed by Ron Barnett. It doesn't map ; precisely to iterations, but it's close. ; init: float sum2 = 0.0 complex zold = (0,0) loop: sum2 = sum2 + exp(-1/cabs(zold-#z)) zold = #z final: #index = sum2 default: title = "Exponential Smoothing (Convergent)" helpfile = "dmj-pub\dmj-pub-uf-es.htm" } dmj-Smooth2General { ; ; This coloring method provides smooth iteration ; colors for all fractal types, convergent or ; divergent (or both). It combines the two methods ; developed by Ron Barnett. It doesn't map ; precisely to iterations, but it's close. ; init: float sum = 0.0 float sum2 = 0.0 complex zold = (0,0) loop: IF (@diverge) sum = sum + exp(-cabs(#z)) ENDIF IF (@converge) sum2 = sum2 + exp(-1/cabs(zold-#z)) ENDIF zold = #z final: IF (|#z - zold| < 0.5) ; convergent bailout. IF (@converge) #index = sum2 ELSE #index = 0 ENDIF ELSE ; divergent bailout. IF (@diverge) #index = sum * @divergescale ELSE #index = 0 ENDIF ENDIF default: title = "Exponential Smoothing (General)" helpfile = "dmj-pub\dmj-pub-uf-es.htm" param diverge caption = "Color Divergent" default = FALSE hint = "If set, points which escape to infinity will be \ colored." endparam param converge caption = "Color Convergent" default = TRUE hint = "If set, points which collapse to one value will be \ colored." endparam param divergescale caption = "Divergent Density" default = 1.0 hint = "Sets the divergent coloring density, relative to the \ convergent coloring. If set to 1.0, they will use \ the same color density." endparam } dmj-Spiral(OUTSIDE) { ; ; This method is an implementation of Ben Weiss's ; decomposition estimator algorithm, as used in ; Kai's Power Tools' Fractal Explorer from ; MetaCreations. Many thanks to Ben for sharing ; this with me. ; ; Note: this version only works correctly with ; the classic Mandelbrot set. Results with other ; fractals will be unpredictable, but could be ; interesting. ; ; Further note: as with KPT, this algorithm is ; not wholly accurate. There are instances (most ; notably on tight spirals) where the radial ; portion of the spiral calculation breaks down, ; and discontinuities result. This is because ; calculating field lines is a difficult task. ; The algorithm used here gives only an approximate ; field line. ; final: float power = 2.0 int i = #numiter complex zn = (0,0) complex p = oldz(i) float d = abs(atan2(p)) / (#pi*2.0) ; initial approximation, in range 0 to 1 ; determine field line (estimated) ; this is the Weiss algorithm WHILE (i > 0) ; still an iteration left to process i = i - 1 ; move to previous iteration zn = oldz(i) ; get orbit value at that iteration p = sqrt(p) ; get positive square root d = d * 0.5 ; all existing angle bits are less significant now IF (real(p)*real(zn)+imag(p)*imag(zn) < 0) ; dot product is negative ; this is an abbreviated calculation ; since we only care about the sign p = -p ; flip sign for next time d = 0.5 - d ; and do some bit mangling ENDIF ENDWHILE IF (imag(#pixel) < 0 ); && @mandel) ; final bit d = 2 - abs(d)*2 ELSE d = d * 2 ENDIF ; determine potential complex il = 1/log(power) ; Inverse log (power). float lp = log(log(@bailout)) ; log(log bailout). float s = real(#numiter + il*lp - il*log(log(cabs(#z)))) ; determine color IF (@potscale < 0) ; reverse spiral #index = @decscale*(2-d) - @potscale*s*0.05 ELSE ; forward spiral #index = @decscale*d + @potscale*s*0.05 ENDIF default: title = "Spiral" helpfile = "dmj-pub\dmj-pub-uf-spiral.htm" param decscale caption = "Radial Scale" default = 1.0 hint = "This is the number of times the gradient should repeat radially; \ this determines the number of 'arms' in the spiral." endparam param potscale caption = "Potential Scale" default = 1.0 hint = "This is the 'tightness' of the spiral, how closely it wraps \ around the fractal." endparam ; param mandel ; caption = "Mandelbrot" ; default = TRUE ; hint = "Check this if you're using a Mandelbrot (not Julia) fractal type." ; endparam ; param power ; caption = "Exponent" ; default = (2,0) ; hint = "This should be set to match the exponent of the \ ; formula you are using. For Mandelbrot, this is 2." ; endparam param bailout caption = "Bailout value" default = 65536.0 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab." endparam } dmj-Trap { ; ; This single type replaces all of the trap types ; in dmj-pub.frm - and includes all the extra UF ; options, so they're more flexible, too. Many ; of the options and shapes in here are built on ; ideas from others, including: ; Cliff Pickover ; Phil Pickard ; Paul Carlson ; Kerry Mitchell ; Luke Plant ; Stephen Ferguson ; Earl Hinrichs ; Janet Preslar ; Mark Townsend ; and probably some others I forgot to credit. ; ; NOTE: SOME THINGS DO NOT WORK YET! ; - Some of the trap modes (sum, average, sign ; average; funky averages) don't look different ; from each other except in distance coloring. ; Unless you're using distance coloring, use ; sum instead of average or sign average, and ; funky average instead of funky average 2. ; init: float closest = 1e38 float closest1 = 1e38 float d = 0.0 float d2 = 0.0 complex point = (0,0) complex point1 = (0,0) complex point2 = (0,0) complex point3 = (0,0) complex z2 = (0,0) BOOL done = false int i = 0 int i1 = 0 int iter = 0 float fiter = @trapstart BOOL trapping = false float diameter2 = sqr(@diameter) float r1 = @angle / 90.0 float r2 = @anglestep / 90.0 float s1 = @skew / 90.0 float s2 = @skewstep / 90.0 complex r = (0,1) ^ r1 complex r0 = (0,0) complex s = (0,1) ^ s1 complex rh = (0,1) ^ (@traporder / 8) ; heart rotation value complex zh = (0,0) complex trapcenter2 = @trapcenter IF (@movetrap) ; Trap Center follows pixel trapcenter2 = #pixel + @trapcenter ENDIF IF (@trapshape >= 17 || @trapshape <= 19) ; ripple shapes diameter2 = #pi / @diameter ENDIF IF (@traptype == 1 || @traptype == 4 || \ @traptype == 5 || @traptype == 7 || \ @traptype == 12 || @traptype == 13 || \ @traptype == 14 || @traptype == 15 || \ @traptype == 16 || @traptype == 17 || \ @traptype == 18) ; last, sum, average, sign average closest = 0.0 ELSEIF (@traptype == 6) ; product closest = 1.0 ELSEIF (@traptype == 9 || @traptype == 11) ; second/two farthest closest = 0.0 closest1 = 0.0 ENDIF BOOL usesolid = true ; assume a solid color loop: iter = iter + 1 ; iteration counter IF (@trapskip != 0) ; we are skipping some iterations fiter = fiter - 1 ; one less to go before we trap WHILE (fiter < 0.0) ; iterations all used up IF (trapping) ; we are currently trapping trapping = false ; so stop fiter = fiter + @trapskip ; skip this many iterations ELSE ; we aren't currently trapping trapping = true ; so start fiter = fiter + @trapiter ; do this many iterations ENDIF ENDWHILE ENDIF IF (@trapskip == 0.0 || trapping) ; if we're checking for traps... ; adjust rotation/skew/trapcenter IF (@anglestep != 0) ; changing rotation angle. r = (0,1) ^ (r1+r2*(iter-1)) ENDIF IF (@skewstep != 0) ; changing skew angle. s = (0,1) ^ (s1+s2*(iter-1)) ENDIF IF (@trapdrift != (0,0)) ; changing trap center trapcenter2 = trapcenter2 + @trapdrift ENDIF IF (@traporbit != (0,0)) ; changing trap orbit trapcenter2 = trapcenter2 * @traporbit ENDIF ; rotate/squash/skew/repeat z IF (@traptype == 18) ; trap only, work on unadulterated pixel z2 = #pixel ELSE ; some other trap mode, use z z2 = #z ENDIF ; 1. radially repeat IF (@gausss > 0 && @gaussr == 0) ; concentrically repeating trap, but not radial z2 = z2 - @gausscenter d = cabs(z2) d = ((d / @gausss) - round(d / @gausss)) * @gausss z2 = z2 * d / cabs(z2) + @gausscenter ENDIF IF (@gaussr > 0.0) ; radially repeating trap z2 = z2 - @gausscenter d = cabs(z2) d2 = atan(imag(z2)/real(z2)) IF (real(z2) < 0) d2 = d2 + #pi ENDIF IF (d2 < 0) d2 = d2 + 2*#pi ENDIF d2 = d2 / (#pi*2) d2 = ((d2 * @gaussr) - round(d2 * @gaussr)) / @gaussr * #pi*2 IF (@gausss > 0) ; concentrically repeating trap d = ((d / @gausss) - round(d / @gausss)) * @gausss ENDIF IF (@radialmode == 0) ; standard radial repeat mode z2 = cos(d2)*d + flip(sin(d2)*d) + @gausscenter ELSEIF (@radialmode == 1) ; unwrap mode z2 = d + flip(d2) ENDIF ENDIF ; 2. grid repeat IF (@gauss > 0) ; repeating trap z2 = z2 - @gaussgcenter z2 = ((z2 / @gauss) - round(z2 / @gauss)) * @gauss z2 = z2 + @gaussgcenter ENDIF ; 3. rotate z2 = (z2 - trapcenter2) * r ; rotate ; 4. apply aspect IF (@aspect != 1.0) z2 = real(z2) + flip(imag(z2) * @aspect) ; apply aspect ENDIF ; 5. skew z2 = real(z2 * s) + flip(imag(z2)) ; apply skew ; determine distance from trap--different for each shape IF (@trapshape == 0) ; point d = cabs(z2) ELSEIF (@trapshape == 1) ; ring d = abs(cabs(z2) - @diameter) ELSEIF (@trapshape == 2) ; ring2 d = abs(|z2| - diameter2) ELSEIF (@trapshape == 3) ; egg d = (cabs(z2-flip(@diameter)*2) + cabs(z2)*@traporder*0.5) * 0.25 ELSEIF (@trapshape == 4) ; hyperbola d = abs(imag(z2) * real(z2) - @diameter) ELSEIF (@trapshape == 5) ; hypercross d = abs(imag(z2) * real(z2)) ELSEIF (@trapshape == 6) ; cross d = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape == 7) ; astroid d = abs(real(z2))^@traporder + abs(imag(z2))^@traporder IF (@traporder < 0) d = 1/d ENDIF ELSEIF (@trapshape == 8) ; diamond d = abs(real(z2)) + abs(imag(z2)) ELSEIF (@trapshape == 9) ; rectangle d = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 > d) d = d2 ENDIF ELSEIF (@trapshape == 10) ; box d = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 > d) d = d2 ENDIF d = abs(d - @diameter) ELSEIF (@trapshape == 11) ; lines d = abs(abs(imag(z2)) - @diameter) ELSEIF (@trapshape == 12) ; waves d = abs(abs(imag(z2) + sin(real(z2)*@trapfreq)*@traporder*0.25) - @diameter) ELSEIF (@trapshape == 13) ; mirrored waves d = abs(abs(imag(z2)) - @diameter + sin(real(z2)*@trapfreq)*@traporder*0.25) ELSEIF (@trapshape == 14) ; mirrored waves 2 d2 = @diameter - sin(real(z2)*@trapfreq)*@traporder*0.25 ; compute wave height d = abs(abs(imag(z2)) - d2) ; distance to each wave d2 = abs(abs(imag(z2)) + d2) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape == 15) ; radial waves d2 = atan2(z2) d = abs(cabs(z2) * (1 - sin(d2*@trapfreq)*@traporder*0.125) - @diameter) ELSEIF (@trapshape == 16) ; radial waves 2 d2 = atan2(z2) d2 = sin(d2*@trapfreq)*@traporder*0.125 d = abs(cabs(z2) * (1 - d2) - @diameter) d2 = abs(cabs(z2) * (1 + d2) - @diameter) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape == 17) ; ring ripples d = cabs(z2) IF (d < @traporder) d = cos(d * diameter2 * @trapfreq) * sqr(1-d/@traporder) ; d = (cos(d * diameter2 * @trapfreq)+1) * sqr(1-d/@traporder) * 0.5 ELSE d = 0 ENDIF ELSEIF (@trapshape == 18) ; grid ripples d = cabs(z2) IF (d < @traporder) d = (cos(real(z2)*diameter2*@trapfreq) + cos(imag(z2)*diameter2*@trapfreq)) * sqr(1-d/@traporder) * 0.5 ELSE d = 0 ENDIF ELSEIF (@trapshape == 19) ; radial ripples d = atan2(z2) d2 = cabs(z2) IF (d2 < @traporder) d = cos(4 * d * @trapfreq) * sqr(1-d2/@traporder) ELSE d = 0 ENDIF ELSEIF (@trapshape == 20) ; pinch d2 = atan2(z2) IF (d2 < 0) d2 = d2 + 2*#pi ENDIF d = sqrt(cabs(z2)) / abs(sin(d2*@traporder*0.5)) ELSEIF (@trapshape == 21) ; spiral d = 1/(cabs(z2)) * @diameter r0 = (0,1) ^ d z2 = z2 * r0 d = atan(abs(imag(z2)/real(z2))) ELSEIF (@trapshape == 22) ; heart zh = real(z2) + flip(abs(imag(z2))) zh = zh*rh * 3 / @diameter d = abs(real(zh) - sqr(imag(zh)) + 3) ENDIF ; apply trap transfer function ; a function is not used because these are floats ; and not all functions apply to floats d = d * @prescale IF (@trapfunc == 1) ; log d = log(d) ELSEIF (@trapfunc == 2) ; sqrt d = sqrt(d) ELSEIF (@trapfunc == 3) ; cuberoot d = (d)^(1/3) ELSEIF (@trapfunc == 4) ; exp d = exp(d) ELSEIF (@trapfunc == 5) ; sqr d = sqr(d) ELSEIF (@trapfunc == 6) ; cube d = (d)^3 ELSEIF (@trapfunc == 7) ; sin d = sin(d) ELSEIF (@trapfunc == 8) ; cos d = cos(d) ELSEIF (@trapfunc == 9) ; tan d = tan(d) ELSEIF (@trapfunc == 10) ; sinc (modified form) IF (d == 0) d = 1 ELSE d = sin(#pi*d)/d ENDIF ENDIF d = d * @postscale IF (@trapabs) ; absolute value only d = abs(d) ENDIF ; now adjust closest/point/i as needed IF (@movetrap == false || iter > 1) IF (@traptype == 0) ; closest IF (d < closest) i = iter point = #z point2 = z2 closest = d ENDIF IF (d < @threshold) usesolid = false ENDIF ELSEIF (@traptype == 1) ; farthest (within threshold) IF (d > closest && d < @threshold) i = iter point = #z point2 = z2 closest = d usesolid = false ENDIF ELSEIF (@traptype == 2) ; first (within threshold) IF (d < @threshold && done == false) i = iter point = #z point2 = z2 closest = d done = true usesolid = false ENDIF ELSEIF (@traptype == 3) ; last (within threshold) IF (d < @threshold) i = iter point = #z point2 = z2 closest = d done = true usesolid = false ENDIF ELSEIF (@traptype == 4) ; sum (within threshold) IF (d < @threshold) i = iter point = point + #z point2 = point2 + z2 closest = closest + d usesolid = false ENDIF ELSEIF (@traptype == 5) ; average (within threshold) IF (d < @threshold) i = iter i1 = i1 + 1 point = point + #z point2 = point2 + z2 closest = closest + d usesolid = false ENDIF ELSEIF (@traptype == 6) ; product (within threshold) IF (d < @threshold) i = iter point = point * #z / @threshold point2 = point2 * z2 / @threshold closest = closest * d / @threshold usesolid = false ENDIF ELSEIF (@traptype == 7) ; sign average IF (d < d2) i = i + 1 point = point + #z point2 = point2 + z2 closest = closest + 1 usesolid = false ELSE i = i - 1 ENDIF d2 = d ELSEIF (@traptype == 8 || @traptype == 10) ; second/two closest IF (d < closest) i1 = i point1 = point point3 = point2 closest1 = closest i = iter point = #z point2 = z2 closest = d ELSEIF (d < closest1) i1 = iter point1 = #z point3 = z2 closest1 = d ENDIF IF (d < @threshold) usesolid = false ENDIF ELSEIF (@traptype == 9 || @traptype == 11) ; second/two farthest IF (d > closest && d < @threshold) i1 = i point1 = point point3 = point2 closest1 = closest i = iter point = #z point2 = z2 closest = d usesolid = false ELSEIF (d > closest1 && d < @threshold) i1 = iter point1 = #z point3 = z2 closest1 = d usesolid = false ENDIF ELSEIF (@traptype == 12) ; funky average IF (d < @threshold) i = i + 1 point = #z - point point2 = z2 - point2 closest = @threshold - abs(closest - d) usesolid = false ENDIF ELSEIF (@traptype == 13) ; funky average 2 IF (d < @threshold) i = i + 1 point = #z - point point2 = z2 - point2 closest = abs(d - @threshold + closest) usesolid = false ENDIF ELSEIF (@traptype == 14) ; funky average 3 (Luke Plant) IF (d < @threshold) i = i + 1 d2 = d/@threshold point = #z + (point-#z) * d2 point2 = z2 + (point2-z2) * d2 closest = closest + d usesolid = false ENDIF ELSEIF (@traptype == 15) ; funky average 4 (exponential average) IF (d < @threshold) i = i + 1 point = #z - point point2 = z2 - point2 closest = closest + exp(-d) usesolid = false ENDIF ELSEIF (@traptype == 16) ; funky average 5 (average distance change) IF (d < d2) point = point + #z point2 = point2 + z2 closest = closest + d2-d usesolid = false ENDIF d2 = d ELSEIF (@traptype == 17) ; funky average 6 (Luke Plant, 1/squared) IF (d < @threshold) i = i + 1 usesolid = false ENDIF d2 = sqr(d/@threshold) point = #z + (point-#z) * d2 point2 = z2 + (point2-z2) * d2 closest = closest + 1/d2 ELSEIF (@traptype == 18) ; trap only, do first iteration IF (iter == 1) point = #z point2 = z2 closest = d/@threshold IF (d < @threshold) usesolid = false ENDIF ENDIF ENDIF ENDIF ENDIF final: ; un-fudge anything that was fudged IF (@traptype == 5) ; traptype average point = point / i1 point2 = point2 / i1 closest = closest / i1 ELSEIF (@traptype == 6) ; traptype product closest = abs(closest) ELSEIF (@traptype == 7) ; traptype sign average point = point / iter point2 = point2 / iter closest = closest / iter ELSEIF (@traptype == 8 || @traptype == 9) ; second closest or farthest i = i - i1 point = point - point1 point2 = point2 - point3 closest = closest - closest1 ELSEIF (@traptype == 10 || @traptype == 11) ; two closest or farthest i = round((i + i1) / 2) point = (point + point1) / 2 point2 = (point2 + point3) / 2 closest = (closest + closest1) / 2 ELSEIF (@traptype == 14) ; funky average 3 closest = @threshold * i - closest ENDIF ; choose coloring based on method IF (!@solidcolor) ; solid color not allowed #solid = false ELSE #solid = usesolid ENDIF IF (@trapcolor == 0) ; distance IF (@traptype == 2 || @traptype == 3) ; first or last type #index = closest / @threshold ELSE ; any other trap type #index = closest ENDIF ELSEIF (@trapcolor == 1) ; magnitude #index = cabs(point2) ELSEIF (@trapcolor == 2) ; real #index = abs(real(point2)) ELSEIF (@trapcolor == 3) ; imaginary #index = abs(imag(point2)) ELSEIF (@trapcolor == 4) ; angle to trap d = atan2(point2) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 5) ; angle to trap 2 (no aspect) point = point - @trapcenter d = atan2(point) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 6) ; angle to origin d = atan2(point) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 7) ; angle to origin 2 (old ReallyCool) #index = 0.02 * abs(atan(imag(point) / real(point)) * 180/#pi) ELSEIF (@trapcolor == 8) ; iteration d = i #index = d / #maxiter ENDIF default: title = "Orbit Traps" helpfile = "dmj-pub\dmj-pub-uf-ot.htm" param trapshape caption = "Trap Shape" default = 0 enum = "point" "ring" "ring 2" "egg" "hyperbola" "hypercross" \ "cross" "astroid" "diamond" "rectangle" "box" "lines" \ "waves" "mirrored waves" "mirrored waves 2" \ "radial waves" "radial waves 2" "ring ripples" \ "grid ripples" "radial ripples" "pinch" "spiral" "heart" ; old ordering: ; enum = "point" "ring" "ring 2" "hyperbola" "hypercross" "cross" \ ; "rectangle" "box" "lines" "spiral" "diamond" "pinch" \ ; "astroid" "heart" "ring ripples" "grid ripples" "egg" \ ; "waves" "mirrored waves" "mirrored waves 2" hint = "This is the shape of the orbit trap." endparam param trapcolor caption = "Trap Coloring" default = 0 enum = "distance" "magnitude" "real" "imaginary" "angle to trap" \ "angle to trap 2" "angle to origin" "angle to origin 2" "iteration" ; old ordering: ; enum = "distance" "magnitude" "real" "imaginary" "angle to trap" \ ; "angle to origin" "iteration" "angle to trap 2" hint = "This is the information used to produce a color." endparam param traptype caption = "Trap Mode" default = 0 enum = "closest" "farthest" "first" "last" "sum" "average" "product" \ "sign average" "second closest" "second farthest" "two closest" \ "two farthest" "funky average" "funky average 2" "funky average 3" \ "funky average 4" "funky average 5" "funky average 6" \ "trap only" hint = "This is how points will be chosen to use for coloring." endparam param traporder caption = "Trap Order" default = 4.0 hint = "Number of leaves for the pinch trap shape, the \ exponent to use for astroid curves (try 0.66667), \ 'egginess', or the height of waves." endparam param trapfreq caption = "Trap Frequency" default = 1.0 hint = "The frequency of ripples or waves." endparam param trapcenter caption = "Trap Center" default = (0,0) hint = "This is the location of the trap in the complex plane." endparam param trapdrift caption = "Trap Drift" default = (0,0) hint = "This is the amount the trap center will move with each \ iteration." endparam param traporbit caption = "Trap Orbit" default = (0,0) hint = "The trap center location is multiplied by this value \ with each iteration." endparam param movetrap caption = "Trap Center Moves" default = false hint = "If this is enabled, the trap center will match the \ pixel location instead of being fixed. This overrides \ Trap Center." endparam param aspect caption = "Aspect Ratio" default = 1.0 min = 0.0000000001 hint = "This is how square the trap is. You can distort the \ trap by using a value other than 1.0." endparam param threshold caption = "Threshold" default = 0.25 min = 0 hint = "This is the width of the trap area, used for most \ trap modes." endparam param diameter caption = "Diameter" default = 1.0 ; min = 0 hint = "This is the diameter of the trap (for ring, box, and \ line shapes)." endparam param angle caption = "Rotation" default = 0.0 hint = "This is the angle, in degrees, that the trap should \ be rotated." endparam param anglestep caption = "Rotation Step" default = 0.0 hint = "This is the angle, in degrees, that the trap will \ be rotated with each iteration." endparam param skew caption = "Skew" default = 0.0 hint = "This is the angle, in degrees, to skew the vertical \ axis of the trap shape." endparam param skewstep caption = "Skew Step" default = 0.0 hint = "This is the angle, in degrees, that the skew angle \ should be increased by with each iteration." endparam param trapfunc caption = "Distance Transfer" default = 0 enum = "linear" "log" "sqrt" "cuberoot" "exp" "sqr" "cube" \ "sin" "cos" "tan" "sinc" hint = "This function is applied to trap distances at every \ iteration. It can be used to apply some special effects \ to trap shapes." endparam param prescale caption = "Distance Pre-scale" default = 1.0 hint = "This is a multiplier applied to the trap distance, \ before the Distance Transfer function is applied to it." endparam param postscale caption = "Distance Post-scale" default = 1.0 hint = "This is a multiplier applied to the trap distance, \ after the Distance Transfer function is applied to it." endparam param trapabs caption = "No Negative Distance" default = false hint = "If enabled, eliminates 'negative' distances by taking \ the absolute value of the distance." endparam param trapstart caption = "Start Iteration" default = 0.0 hint = "This is the iteration at which to start watching for \ orbit traps." endparam param trapiter caption = "Trap Iterations" default = 10000.0 hint = "This is the number of iterations to watch for traps \ in a single block." endparam param trapskip caption = "Skip Iterations" default = 0.0 hint = "This is the number of iterations to skip watching for \ traps, after a block of watched iterations." endparam param gauss caption = "Repeat Spacing" default = 0.0 min = 0.0 hint = "This is the spacing at which the trap will repeat, \ all across the complex plane. If zero, trap does \ not repeat." endparam param gaussgcenter caption = "Repeat Center" default = (0,0) hint = "This is the center point of grid trap repetition." endparam param gaussr caption = "Radial Spacing" default = 0.0 min = 0 hint = "This is the number of times the trap will repeat, \ radially around the radial center point. If zero, \ trap does not repeat radially." endparam param gausss caption = "Concentric Spacing" default = 0.0 min = 0 hint = "This is the distance at which the trap will repeat, \ radially out from the radial center point. If zero, \ trap does not repeat radially out." endparam param gausscenter caption = "Radial Center" default = (0,0) hint = "This is the center point of radial trap repetition." endparam param radialmode caption = "Radial Mode" default = 0 enum = "kaleidoscope" "polar coordinates" hint = "Indicates the type of radial repetition to use. The \ two styles produce very different results." endparam param solidcolor caption = "Use Solid Color" default = false hint = "If enabled, areas 'outside' the trap area will be colored \ with the 'solid color' on the coloring tab." endparam } dmj-Triangle { ; ; This coloring method implements a variation of Kerry ; Mitchell's Triangle Inequality Average coloring ; method. Because the smoothing used in this formula ; is based on the dmj-Smooth formula, which only works ; for z^n+c and derivates, the smoothing here will only ; work for those types as well. ; init: float sum = 0.0 float sum2 = 0.0 float ac = cabs(#pixel) float il = 1/log(@power) float lp = log(log(@bailout)/2.0) float az2 = 0.0 float lowbound = 0.0 float f = 0.0 BOOL first = true float ipower = 1/@apower loop: sum2 = sum IF (!first) az2 = cabs(#z - #pixel) lowbound = abs(az2 - ac) IF (@aflavor == 0) sum = sum + ((cabs(#z) - lowbound) / (az2+ac - lowbound))^@apower ELSEIF (@aflavor == 1) sum = sum + 1-(1-(cabs(#z) - lowbound) / (az2+ac - lowbound))^ipower ENDIF ELSE first = false ENDIF final: sum = sum / (#numiter) sum2 = sum2 / (#numiter-1) f = il*lp - il*log(log(cabs(#z))) #index = sum2 + (sum-sum2) * (f+1) default: title = "Triangle Inequality Average" helpfile = "dmj-pub\dmj-pub-uf-tia.htm" param apower caption = "Average Exponent" default = 1.0 hint = "This skews the values averaged by raising them to \ this power. Use 1.0 for the classic coloring." endparam param aflavor caption = "Average Flavor" default = 0 enum = "normal" "reversed" hint = "Controls whether values are reversed before being \ raised to a power. Has no effect if Average Exponent \ is 1.0." endparam param power caption = "Exponent" default = 2.0 hint = "This should be set to match the exponent of the \ formula you are using. For Mandelbrot, this is 2." endparam param bailout caption = "Bailout" default = 1e20 min = 1 hint = "This should be set to match the bailout value in \ the Formula tab. Use a very high bailout!" endparam } dmj-TwinTrap { ; ; This formula is based in large part on the ; ideas of Kerry Mitchell, in his dual-trap ; formula. What it does is set up *two* trap ; shapes, and color based on the interaction ; between these shapes and the fractal orbit. ; ; The idea is to position two trap shapes, ; and measure two trap distances at each ; iteration. From these two trap distances, a ; final trap distance is computed; this is ; then used for the regular trap modes. ; init: float closest = 1e38 float closest1 = 1e38 float d = 0.0 float d1 = 0.0 float d2 = 0.0 complex point = (0,0) complex point1 = (0,0) complex point2 = (0,0) complex point3 = (0,0) complex point4 = (0,0) complex point5 = (0,0) complex point6 = (0,0) complex point7 = (0,0) complex z1 = (0,0) complex z2 = (0,0) complex z3 = (0,0) BOOL done = false int i = 0 int i1 = 0 int iter = 0 float fiter = @trapstart BOOL trapping = false float diametersq1 = sqr(@diameter1) float diametersq2 = sqr(@diameter2) complex r1 = (0,1) ^ (@angle1 / 90.0) complex r2 = (0,1) ^ (@angle2 / 90.0) complex r0 = (0,0) complex rh1 = (0,1) ^ (@traporder1 / 8) ; heart rotation value complex rh2 = (0,1) ^ (@traporder2 / 8) ; heart rotation value complex zh = (0,0) BOOL ab = false IF (@trapshape1 >= 17 || @trapshape1 <= 19) ; ripple shapes diametersq1 = #pi / @diameter1 ENDIF IF (@trapshape2 >= 17 || @trapshape2 <= 19) ; ripple shapes diametersq2 = #pi / @diameter2 ENDIF IF (@traptype == 1 || @traptype == 4 || \ @traptype == 5 || @traptype == 7 || \ @traptype == 12 || @traptype == 13 || \ @traptype == 14 || @traptype == 15 || \ @traptype == 16 || @traptype == 17 || \ @traptype == 18) ; last, sum, average, sign average closest = 0.0 ELSEIF (@traptype == 6) ; product closest = 1.0 ELSEIF (@traptype == 9 || @traptype == 11) ; second/two farthest closest = 0.0 closest1 = 0.0 ENDIF loop: iter = iter + 1 ; iteration counter IF (@trapskip != 0) ; we are skipping some iterations fiter = fiter - 1 ; one less to go before we trap WHILE (fiter < 0.0) ; iterations all used up IF (trapping) ; we are currently trapping trapping = false ; so stop fiter = fiter + @trapskip ; skip this many iterations ELSE ; we aren't currently trapping trapping = true ; so start fiter = fiter + @trapiter ; do this many iterations ENDIF ENDWHILE ENDIF IF (@trapskip == 0.0 || trapping) ; if we're checking for traps... ; adjust rotation/skew/trapcenter ; (we don't do this for this trap formula) ; rotate/squash/skew/repeat z IF (@traptype == 18) ; trap only, work on unadulterated pixel z1 = #pixel z2 = #pixel ELSE ; some other trap mode, use z z1 = #z z2 = #z ENDIF ; 1. radially repeat ; (we don't do this for this trap formula) ; 2. grid repeat ; (we don't do this for this trap formula) ; 3. rotate z1 = (z1 - @trapcenter1) * r1 z2 = (z2 - @trapcenter2) * r2 ; rotate ; 4. apply aspect IF (@aspect1 != 1.0) z1 = real(z1) + flip(imag(z1) * @aspect1) ENDIF IF (@aspect2 != 1.0) z2 = real(z2) + flip(imag(z2) * @aspect2) ; apply aspect ENDIF ; 5. skew ; (we don't do this for this trap formula) ; determine distance from trap--different for each shape ; first trap shape IF (@trapshape1 == 0) ; point d = cabs(z1) ELSEIF (@trapshape1 == 1) ; ring d = abs(cabs(z1) - @diameter1) ELSEIF (@trapshape1 == 2) ; ring2 d = abs(|z1| - diametersq1) ELSEIF (@trapshape1 == 3) ; egg d = (cabs(z1-flip(@diameter1)*2) + cabs(z1)*@traporder1*0.5) * 0.25 ELSEIF (@trapshape1 == 4) ; hyperbola d = abs(imag(z1) * real(z1) - @diameter1) ELSEIF (@trapshape1 == 5) ; hypercross d = abs(imag(z1) * real(z1)) ELSEIF (@trapshape1 == 6) ; cross d = abs(real(z1)) d2 = abs(imag(z1)) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape1 == 7) ; astroid d = abs(real(z1))^@traporder1 + abs(imag(z1))^@traporder1 IF (@traporder1 < 0) d = 1/d ENDIF ELSEIF (@trapshape1 == 8) ; diamond d = abs(real(z1)) + abs(imag(z1)) ELSEIF (@trapshape1 == 9) ; rectangle d = abs(real(z1)) d2 = abs(imag(z1)) IF (d2 > d) d = d2 ENDIF ELSEIF (@trapshape1 == 10) ; box d = abs(real(z1)) d2 = abs(imag(z1)) IF (d2 > d) d = d2 ENDIF d = abs(d - @diameter1) ELSEIF (@trapshape1 == 11) ; lines d = abs(abs(imag(z1)) - @diameter1) ELSEIF (@trapshape1 == 12) ; waves d = abs(abs(imag(z1) + sin(real(z1)*@trapfreq1)*@traporder1*0.25) - @diameter1) ELSEIF (@trapshape1 == 13) ; mirrored waves d = abs(abs(imag(z1)) - @diameter1 + sin(real(z1)*@trapfreq1)*@traporder1*0.25) ELSEIF (@trapshape1 == 14) ; mirrored waves 2 d2 = @diameter1 - sin(real(z1)*@trapfreq1)*@traporder1*0.25 ; compute wave height d = abs(abs(imag(z1)) - d2) ; distance to each wave d2 = abs(abs(imag(z1)) + d2) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape1 == 15) ; radial waves d2 = atan2(z1) d = abs(cabs(z1) * (1 - sin(d2*@trapfreq1)*@traporder1*0.125) - @diameter1) ELSEIF (@trapshape1 == 16) ; radial waves 2 d2 = atan2(z1) d2 = sin(d2*@trapfreq1)*@traporder1*0.125 d = abs(cabs(z1) * (1 - d2) - @diameter1) d2 = abs(cabs(z1) * (1 + d2) - @diameter1) IF (d2 < d) d = d2 ENDIF ELSEIF (@trapshape1 == 17) ; ring ripples d = cabs(z1) IF (d < @traporder1) d = cos(d * diametersq1 * @trapfreq1) * sqr(1-d/@traporder1) ELSE d = 0 ENDIF ELSEIF (@trapshape1 == 18) ; grid ripples d = cabs(z1) IF (d < @traporder1) d = (cos(real(z1)*diametersq1*@trapfreq1) + cos(imag(z1)*diametersq1*@trapfreq1)) * sqr(1-d/@traporder1) * 0.5 ELSE d = 0 ENDIF ELSEIF (@trapshape1 == 19) ; radial ripples d = atan2(z1) d2 = cabs(z1) IF (d2 < @traporder1) d = cos(4 * d * @trapfreq1) * sqr(1-d2/@traporder1) ELSE d = 0 ENDIF ELSEIF (@trapshape1 == 20) ; pinch d2 = atan2(z1) IF (d2 < 0) d2 = d2 + 2*#pi ENDIF d = sqrt(cabs(z1)) / abs(sin(d2*@traporder1*0.5)) ELSEIF (@trapshape1 == 21) ; spiral d = 1/(cabs(z1)) * @diameter1 r0 = (0,1) ^ d z1 = z1 * r0 d = atan(abs(imag(z1)/real(z1))) ELSEIF (@trapshape1 == 22) ; heart zh = real(z1) + flip(abs(imag(z1))) zh = zh*rh1 * 3 / @diameter1 d = abs(real(zh) - sqr(imag(zh)) + 3) ENDIF ; second trap shape IF (@trapshape2 == 0) ; point d1 = cabs(z2) ELSEIF (@trapshape2 == 1) ; ring d1 = abs(cabs(z2) - @diameter2) ELSEIF (@trapshape2 == 2) ; ring2 d1 = abs(|z2| - diametersq2) ELSEIF (@trapshape2 == 3) ; egg d1 = (cabs(z2-flip(@diameter2)*2) + cabs(z2)*@traporder2*0.5) * 0.25 ELSEIF (@trapshape2 == 4) ; hyperbola d1 = abs(imag(z2) * real(z2) - @diameter2) ELSEIF (@trapshape2 == 5) ; hypercross d1 = abs(imag(z2) * real(z2)) ELSEIF (@trapshape2 == 6) ; cross d1 = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 < d1) d1 = d2 ENDIF ELSEIF (@trapshape2 == 7) ; astroid d1 = abs(real(z2))^@traporder2 + abs(imag(z2))^@traporder2 IF (@traporder2 < 0) d1 = 1/d1 ENDIF ELSEIF (@trapshape2 == 8) ; diamond d1 = abs(real(z2)) + abs(imag(z2)) ELSEIF (@trapshape2 == 9) ; rectangle d1 = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 > d1) d1 = d2 ENDIF ELSEIF (@trapshape2 == 10) ; box d1 = abs(real(z2)) d2 = abs(imag(z2)) IF (d2 > d1) d1 = d2 ENDIF d1 = abs(d1 - @diameter2) ELSEIF (@trapshape2 == 11) ; lines d1 = abs(abs(imag(z2)) - @diameter2) ELSEIF (@trapshape2 == 12) ; waves d1 = abs(abs(imag(z2) + sin(real(z2)*@trapfreq2)*@traporder2*0.25) - @diameter2) ELSEIF (@trapshape2 == 13) ; mirrored waves d1 = abs(abs(imag(z2)) - @diameter2 + sin(real(z2)*@trapfreq2)*@traporder2*0.25) ELSEIF (@trapshape2 == 14) ; mirrored waves 2 d2 = @diameter2 - sin(real(z2)*@trapfreq2)*@traporder2*0.25 ; compute wave height d1 = abs(abs(imag(z2)) - d2) ; distance to each wave d2 = abs(abs(imag(z2)) + d2) IF (d2 < d1) d1 = d2 ENDIF ELSEIF (@trapshape2 == 15) ; radial waves d2 = atan2(z2) d1 = abs(cabs(z2) * (1 - sin(d2*@trapfreq2)*@traporder2*0.125) - @diameter2) ELSEIF (@trapshape2 == 16) ; radial waves 2 d2 = atan2(z2) d2 = sin(d2*@trapfreq2)*@traporder2*0.125 d1 = abs(cabs(z2) * (1 - d2) - @diameter2) d2 = abs(cabs(z2) * (1 + d2) - @diameter2) IF (d2 < d1) d1 = d2 ENDIF ELSEIF (@trapshape2 == 17) ; ring ripples d1 = cabs(z2) IF (d1 < @traporder2) d1 = cos(d1 * diametersq2 * @trapfreq2) * sqr(1-d1/@traporder2) ELSE d1 = 0 ENDIF ELSEIF (@trapshape2 == 18) ; grid ripples d1 = cabs(z2) IF (d1 < @traporder2) d1 = (cos(real(z2)*diametersq2*@trapfreq2) + cos(imag(z2)*diametersq2*@trapfreq2)) * sqr(1-d1/@traporder2) * 0.5 ELSE d1 = 0 ENDIF ELSEIF (@trapshape2 == 19) ; radial ripples d1 = atan2(z2) d2 = cabs(z2) IF (d2 < @traporder2) d1 = cos(4 * d1 * @trapfreq2) * sqr(1-d2/@traporder2) ELSE d1 = 0 ENDIF ELSEIF (@trapshape2 == 20) ; pinch d2 = atan2(z2) IF (d2 < 0) d2 = d2 + 2*#pi ENDIF d1 = sqrt(cabs(z2)) / abs(sin(d2*@traporder2*0.5)) ELSEIF (@trapshape2 == 21) ; spiral d1 = 1/(cabs(z2)) * @diameter2 r0 = (0,1) ^ d1 z2 = z2 * r0 d1 = atan(abs(imag(z2)/real(z2))) ELSEIF (@trapshape2 == 22) ; heart zh = real(z2) + flip(abs(imag(z2))) zh = zh*rh2 * 3 / @diameter2 d1 = abs(real(zh) - sqr(imag(zh)) + 3) ENDIF ; now combine the two trap distances in some fashion. d = d * @trapweight1 ; apply trap weights d1 = d1 * @trapweight2 z3 = d + flip(d1) ; make synthetic point IF (@trapmerge == 0) ; closest IF (d1 < d) d = d1 ENDIF ELSEIF (@trapmerge == 1) ; farthest IF (d1 > d) d = d1 ENDIF ELSEIF (@trapmerge == 2) ; A+B d = d + d1 ELSEIF (@trapmerge == 3) ; A-B d = d - d1 IF (d < 0) d = 0 ENDIF ELSEIF (@trapmerge == 4) ; B-A d = d1 - d IF (d < 0) d = 0 ENDIF ELSEIF (@trapmerge == 5) ; |A-B| d = abs(d - d1) ELSEIF (@trapmerge == 6) ; A*A+B*B d = sqr(d)+sqr(d1) ELSEIF (@trapmerge == 7) ; A*A-B*B d = sqr(d)-sqr(d1) IF (d < 0) d = 0 ENDIF ELSEIF (@trapmerge == 8) ; B*B-A*A d = sqr(d1)-sqr(d) IF (d < 0) d = 0 ENDIF ELSEIF (@trapmerge == 9) ; |A*A-B*B| d = abs(sqr(d)-sqr(d1)) ELSEIF (@trapmerge == 10) ; A*A*A+B*B*B d = d^3+d1^3 ELSEIF (@trapmerge == 11) ; A*A*A-B*B*B d = d^3-d1^3 ELSEIF (@trapmerge == 12) ; B*B*B-A*A*A d = d1^3-d^3 ELSEIF (@trapmerge == 13) ; |A*A*A-B*B*B| d = abs(d^3-d1^3) ELSEIF (@trapmerge == 14) ; A*B d = d*d1 ELSEIF (@trapmerge == 15) ; A/B d = d/d1 ELSEIF (@trapmerge == 16) ; B/A d = d1/d ELSEIF (@trapmerge == 17) ; 1/(A*B) d = 1/(d*d1) ELSEIF (@trapmerge == 19) ; B only d = d1 ELSEIF (@trapmerge == 20) ; A, then B IF (ab) d = d1 ENDIF ab = !ab ENDIF ; now adjust closest/point/i as needed IF (@traptype == 0) ; closest IF (d < closest) i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d ENDIF ELSEIF (@traptype == 1) ; farthest (within threshold) IF (d > closest && d < @threshold) i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d ENDIF ELSEIF (@traptype == 2) ; first (within threshold) IF (d < @threshold && done == false) i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d done = true ENDIF ELSEIF (@traptype == 3) ; last (within threshold) IF (d < @threshold) i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d done = true ENDIF ELSEIF (@traptype == 4) ; sum (within threshold) IF (d < @threshold) i = iter point = point + #z point2 = point2 + z1 point4 = point4 + z2 point6 = point6 + z3 closest = closest + d ENDIF ELSEIF (@traptype == 5) ; average (within threshold) IF (d < @threshold) i = iter i1 = i1 + 1 point = point + #z point2 = point2 + z1 point4 = point4 + z2 point6 = point6 + z3 closest = closest + d ENDIF ELSEIF (@traptype == 6) ; product (within threshold) IF (d < @threshold) i = iter point = point * #z / @threshold point2 = point2 * z1 / @threshold point4 = point4 * z2 / @threshold point6 = point6 * z3 / @threshold closest = closest * d / @threshold ENDIF ELSEIF (@traptype == 7) ; sign average IF (d < d2) i = i + 1 point = point + #z point2 = point2 + z1 point4 = point4 + z2 point6 = point6 + z3 closest = closest + 1 ELSE i = i - 1 ENDIF d2 = d ELSEIF (@traptype == 8 || @traptype == 10) ; second/two closest IF (d < closest) i1 = i point1 = point point3 = point2 point5 = point4 point7 = point6 closest1 = closest i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d ELSEIF (d < closest1) i1 = iter point1 = #z point3 = z1 point5 = z2 point7 = z3 closest1 = d ENDIF ELSEIF (@traptype == 9 || @traptype == 11) ; second/two farthest IF (d > closest && d < @threshold) i1 = i point1 = point point3 = point2 point5 = point4 point7 = point6 closest1 = closest i = iter point = #z point2 = z1 point4 = z2 point6 = z3 closest = d ELSEIF (d > closest1 && d < @threshold) i1 = iter point1 = #z point3 = z1 point5 = z2 point7 = z3 closest1 = d ENDIF ELSEIF (@traptype == 12) ; funky average IF (d < @threshold) i = i + 1 point = #z - point point2 = z1 - point2 point4 = z2 - point4 point6 = z3 - point6 closest = @threshold - abs(closest - d) ENDIF ELSEIF (@traptype == 13) ; funky average 2 IF (d < @threshold) i = i + 1 point = #z - point point2 = z1 - point2 point4 = z2 - point4 point6 = z3 - point6 closest = abs(d - @threshold + closest) ENDIF ELSEIF (@traptype == 14) ; funky average 3 (Luke Plant) IF (d < @threshold) i = i + 1 d1 = d/@threshold point = #z + (point-#z) * d1 point2 = z1 + (point2-z1) * d1 point4 = z2 + (point4-z2) * d1 point6 = z3 + (point6-z3) * d1 closest = closest + d ENDIF ELSEIF (@traptype == 15) ; funky average 4 (exponential average) IF (d < @threshold) i = i + 1 point = #z - point point2 = z1 - point2 point4 = z2 - point4 point6 = z3 - point6 closest = closest + exp(-d) ENDIF ELSEIF (@traptype == 16) ; funky average 5 (average distance change) IF (d < d2) point = point + #z point2 = point2 + z1 point4 = point4 + z2 point6 = point6 + z3 closest = closest + d2-d ENDIF d2 = d ELSEIF (@traptype == 17) ; funky average 6 (Luke Plant, 1/squared) IF (d < @threshold) i = i + 1 ENDIF d1 = sqr(d/@threshold) point = #z + (point-#z) * d1 point2 = z1 + (point2-z1) * d1 point4 = z2 + (point4-z2) * d1 point6 = z3 + (point6-z3) * d1 closest = closest + 1/d1 ELSEIF (@traptype == 18) ; trap only, do first iteration IF (iter == 1) point = #z point2 = z1 point4 = z2 point6 = z3 closest = d/@threshold ENDIF ENDIF ENDIF final: ; un-fudge anything that was fudged IF (@traptype == 5) ; traptype average point = point / i1 point2 = point2 / i1 point4 = point4 / i1 point6 = point6 / i1 closest = closest / i1 ELSEIF (@traptype == 6) ; traptype product closest = abs(closest) ELSEIF (@traptype == 7) ; traptype sign average point = point / iter point2 = point2 / iter point4 = point4 / iter point6 = point6 / iter closest = closest / iter ELSEIF (@traptype == 8 || @traptype == 9) ; second closest or farthest i = i - i1 point = point - point1 point2 = point2 - point3 point4 = point4 - point5 point6 = point6 - point7 closest = closest - closest1 ELSEIF (@traptype == 10 || @traptype == 11) ; two closest or farthest i = round((i + i1) / 2) point = (point + point1) / 2 point2 = (point2 + point3) / 2 point4 = (point4 + point5) / 2 point6 = (point6 + point7) / 2 closest = (closest + closest1) / 2 ELSEIF (@traptype == 14) ; funky average 3 closest = @threshold * i - closest ENDIF ; choose coloring based on method IF (@trapcolor == 0) ; distance IF (@traptype == 2 || @traptype == 3) ; first or last type #index = closest / @threshold ELSE ; any other trap type #index = closest ENDIF ELSEIF (@trapcolor == 1) ; magnitude #index = cabs(point) ELSEIF (@trapcolor == 2) ; real #index = abs(real(point)) ELSEIF (@trapcolor == 3) ; imaginary #index = abs(imag(point)) ELSEIF (@trapcolor == 4) ; angle to origin d = atan2(point) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 5) ; angle to origin 2 (old ReallyCool) #index = 0.02 * abs(atan(imag(point) / real(point)) * 180/#pi) ELSEIF (@trapcolor == 6) ; distance A IF (@traptype == 2 || @traptype == 3) ; first or last type #index = real(point6) / @threshold ELSE ; any other trap type #index = real(point6) ENDIF ELSEIF (@trapcolor == 7) ; magnitude A #index = cabs(point2) ELSEIF (@trapcolor == 8) ; real A #index = abs(real(point2)) ELSEIF (@trapcolor == 9) ; imaginary A #index = abs(imag(point2)) ELSEIF (@trapcolor == 10) ; angle to trap A d = atan2(point2) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 11) ; angle to trap A 2 (no aspect) point = point - @trapcenter1 d = atan2(point) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 12) ; distance B IF (@traptype == 2 || @traptype == 3) ; first or last type #index = imag(point6) / @threshold ELSE ; any other trap type #index = imag(point6) ENDIF ELSEIF (@trapcolor == 13) ; magnitude B #index = cabs(point4) ELSEIF (@trapcolor == 14) ; real B #index = abs(real(point4)) ELSEIF (@trapcolor == 15) ; imaginary B #index = abs(imag(point4)) ELSEIF (@trapcolor == 16) ; angle to trap B d = atan2(point4) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 17) ; angle to trap B 2 (no aspect) point = point - @trapcenter2 d = atan2(point) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 18) ; synthetic magnitude #index = cabs(point6) ELSEIF (@trapcolor == 19) ; synthetic angle d = atan2(point6) IF (d < 0) d = d + #pi * 2 ENDIF #index = d / (#pi * 2) ELSEIF (@trapcolor == 20) ; iteration d = i #index = d / #maxiter ENDIF default: title = "Twin Traps" helpfile = "dmj-pub\dmj-pub-uf-tt.htm" param trapshape1 caption = "Trap A Shape" default = 1 enum = "point" "ring" "ring 2" "egg" "hyperbola" "hypercross" \ "cross" "astroid" "diamond" "rectangle" "box" "lines" \ "waves" "mirrored waves" "mirrored waves 2" \ "radial waves" "radial waves 2" "ring ripples" \ "grid ripples" "radial ripples" "pinch" "spiral" "heart" hint = "This is the shape of the orbit trap." endparam param trapshape2 caption = "Trap B Shape" default = 6 enum = "point" "ring" "ring 2" "egg" "hyperbola" "hypercross" \ "cross" "astroid" "diamond" "rectangle" "box" "lines" \ "waves" "mirrored waves" "mirrored waves 2" \ "radial waves" "radial waves 2" "ring ripples" \ "grid ripples" "radial ripples" "pinch" "spiral" "heart" hint = "This is the shape of the orbit trap." endparam param trapmerge caption = "Trap Merging" default = 0 enum = "closest" "farthest" "A + B" "A - B" "B - A" "| A - B |" \ "A*A + B*B" "A*A - B*B" "B*B - A*A" "| A*A - B*B |" \ "A*A*A + B*B*B" "A*A*A - B*B*B" "B*B*B - A*A*A" "| A*A*A - B*B*B |" \ "A * B" "A / B" "B / A" "1 / (A*B)" "A only" "B only" "alternating" hint = "Sets how the two trap distances (A and B) will be combined \ before being checked by the regular trap modes." endparam param trapweight1 caption = "Trap A Weight" default = 1.0 hint = "This is the relative weight of the distance to trap A, \ prior to applying the Trap Merging parameter." endparam param trapweight2 caption = "Trap B Weight" default = 1.0 hint = "This is the relative weight of the distance to trap B, \ prior to applying the Trap Merging parameter." endparam param trapcolor caption = "Trap Coloring" default = 0 enum = "distance" "magnitude" "real" "imaginary" "angle to origin" \ "angle to origin 2" "distance A" "magnitude A" "real A" \ "imaginary A" "angle to trap A" "angle to trap A 2" \ "distance B" "magnitude B" "real B" "imaginary B" \ "angle to trap B" "angle to trap B 2" "synthetic magnitude" \ "synthetic angle" "iteration" hint = "This is the information used to produce a color." endparam param traptype caption = "Trap Mode" default = 0 enum = "closest" "farthest" "first" "last" "sum" "average" "product" \ "sign average" "second closest" "second farthest" "two closest" \ "two farthest" "funky average" "funky average 2" "funky average 3" \ "funky average 4" "funky average 5" "funky average 6" \ "trap only" hint = "This is how points will be chosen to use for coloring." endparam param traporder1 caption = "Trap A Order" default = 4.0 hint = "Number of leaves for the pinch trap shape, the \ exponent to use for astroid curves (try 0.66667), \ 'egginess', or the height of waves." endparam param traporder2 caption = "Trap B Order" default = 4.0 hint = "Number of leaves for the pinch trap shape, the \ exponent to use for astroid curves (try 0.66667), \ 'egginess', or the height of waves." endparam param trapfreq1 caption = "Trap A Frequency" default = 1.0 hint = "The frequency of ripples or waves." endparam param trapfreq2 caption = "Trap B Frequency" default = 1.0 hint = "The frequency of ripples or waves." endparam param trapcenter1 caption = "Trap A Center" default = (0,0) hint = "This is the location of the trap in the complex plane." endparam param trapcenter2 caption = "Trap B Center" default = (0,0) hint = "This is the location of the trap in the complex plane." endparam param aspect1 caption = "Aspect Ratio A" default = 1.0 min = 0.0000000001 hint = "This is how square the trap is. You can distort the \ trap by using a value other than 1.0." endparam param aspect2 caption = "Aspect Ratio B" default = 1.0 min = 0.0000000001 hint = "This is how square the trap is. You can distort the \ trap by using a value other than 1.0." endparam param threshold caption = "Threshold" default = 0.25 min = 0 hint = "This is the width of the trap area, used for most \ trap modes." endparam param diameter1 caption = "Diameter A" default = 1.0 hint = "This is the diameter of the trap (for ring, box, and \ line shapes)." endparam param diameter2 caption = "Diameter B" default = 1.0 hint = "This is the diameter of the trap (for ring, box, and \ line shapes)." endparam param angle1 caption = "Rotation A" default = 0.0 hint = "This is the angle, in degrees, that the trap should \ be rotated." endparam param angle2 caption = "Rotation B" default = 0.0 hint = "This is the angle, in degrees, that the trap should \ be rotated." endparam ; param anglestep ; caption = "Rotation Step" ; default = 0.0 ; hint = "This is the angle, in degrees, that the trap will \ ; be rotated with each iteration." ; endparam ; param skew ; caption = "Skew" ; default = 0.0 ; hint = "This is the angle, in degrees, to skew the vertical \ ; axis of the trap shape." ; endparam ; param skewstep ; caption = "Skew Step" ; default = 0.0 ; hint = "This is the angle, in degrees, that the skew angle \ ; should be increased by with each iteration." ; endparam param trapstart caption = "Start Iteration" default = 0.0 hint = "This is the iteration at which to start watching for \ orbit traps." endparam param trapiter caption = "Trap Iterations" default = 10000.0 hint = "This is the number of iterations to watch for traps \ in a single block." endparam param trapskip caption = "Skip Iterations" default = 0.0 hint = "This is the number of iterations to skip watching for \ traps, after a block of watched iterations." endparam }