comment { Fractal formulas by Andreas Hoppler \$LastChangedDate: 2007-07-07 17:06:19 +0200 (Sa, 07 Jul 2007) \$ 2004-02-01 added more hints and comments 2004-02-02 Separated inner and outer bailout where applicable 2004-02-02 Added custom orbit start in most M-type formulas 2004-02-08 Added 'frothy_basin' formula 2004-02-15 Some fixes to 'frothy basin' formula 2004-03-21 Added Hénon, Perlin, Poly_3 formulas. Removed version numbers from titles. 2007-07-01 Added 'Circle-Logo' and 'Lacunary' formulas, updated 'hring' 2007-07-03 Re-renamed 'cvalue' parameter back to 'cp' for hring-formulas to keep old params working 2007-07-03 Fixed a bug with hring M-c and stripes 2007-07-10 Editorial changes to 'Lacunary3_J' and '_M' 2021-12-27 Added perturbation calculation to CircleLogo_M and _J 2022-01-02 Removed spurious parameter on CircleLogo_M. Also removed DEBUG defines 2022-02-03 Fixed Circle Logo perturbation bug for power 5 } ##################################################################### ### NEWTON'S METHOD FRACTALS ### ##################################################################### # --- Polynomials of degree 3 --------------------------------------- N_poly3_J { ; newton's method applied to P(z) = (z + 1) * (z - 0.5 - a) * (z - 0.5 + a) ; Complex plane fractal (J-set) for a specific 'a' ; F(z) = z - P(z)/P'(z) ; Attractors: the 3 zeros of P(z): -1, 0.5+a, 0.5-a ; Critical points: 0, as well the 3 zeros of P(z). global: ; the attractors complex attract[3] int n_attract = 3 complex cf1 ; pre-calculate some factors outside the loop complex cf2 if @pchoice=="distance" attract[0] = -1 attract[1] = (1-@a)/2 attract[2] = (1+@a)/2 cf1 = -(3+Sqr(@a))/4 cf2 = (Sqr(@a) - 1)/4 elseif @pchoice=="halfdist" ; for some older param files attract[0] = -1 attract[1] = 0.5 - @a attract[2] = 0.5 + @a cf1 = -0.75 - Sqr(@a) cf2 = Sqr(@a) - 0.25 elseif @pchoice=="square" attract[0] = -1 complex tmp = sqrt(@a)/2 attract[1] = 0.5 - tmp attract[2] = 0.5 + tmp cf1 = -(3+@a)/4 cf2 = (@a - 1)/4 else ; @pChoice=="zero" attract[0] = -1 attract[1] = 1 attract[2] = @a cf1 = -2 * @a cf2 = -@a endif init: z=#pixel complex zprev = 0 int step=0 ; which attractor are we checking ? (with @stripes) int found=-1 ; did we find one? loop: complex dist ; distance to attractor if (@stripes) ; stripes: check only one attractor dist = z-attract[step] if |dist| < @delta found=step else step = step+1 if step>=n_attract step = 0 endif endif else ; no stripes int i=0 ; no stripes: check all attractors repeat dist = z-attract[i] if |dist| < @delta found = i endif i=i+1 until i>=n_attract || found>=0 endif ; stripes if (found<0) && (!@stripes || step==0) ; here we step z (this is the core of the formula) zprev = z ; z[new] = z - P(z)/P'(z) = (z*P'(z) - P(z)) / P'(z) if @pChoice=="zero" z = ((2 * z + cf2) * z^2 + cf2) / (( 3 * z + cf1 ) * z - 1) else z = ( 2 * z^3 + cf2 ) / ( 3 * z^2 + cf1 ) endif endif ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm) if @zmode != "z" if found>=0 if @zmode=="distance" z = dist elseif @zmode=="step" z = z-zprev endif ; zmode endif ; found endif ; zmode bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue found<0 default: title = "Newton-Julia, degree-3 polynomial" ; v1.3 rating = recommended periodicity = 0 ; turn off periodicity checking by default as it interferes with 'stripes' heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="distance" "halfdist" "square" "zero" default=0 hint="Governs the meaning of the parameter 'a'. \ 'distance' means 'a' is the distance between the two zeros. The third zero is -1. \ 'halfdist' means 'a' is half the distance. \ 'square' means sqrt(a) is the distance. \ 'zero' means 'a' is the location of a zero, the other two zeros are -1 and 1." endparam complex param a caption="a" default=flip(sqrt(3)) hint="See 'param interpretation' above for the meaning of this parameter." endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="How close can z get to a fixed point? \ Actually, this is the square of the distance. \ Smaller values mean more iterations." default=1e-10 min=0 endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 3 stripes." default=false endparam param zmode caption = "final z" enum = "z" "distance" "step" hint="'distance' replaces the final z value with the distance from the (known) attractor. \ 'step' replaces the final z value with the distance from the previous z. \ Since Newton's method converges quadratically, these choices go towards zero in a quadratic way, \ therefore you can use a smooth coloring method (Smooth Mandel / Smooth Stripes / ...) \ with power=2." default = 1 endparam switch: type = "N_poly3_M" delta = delta pchoice = pchoice stripes = stripes zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_poly3_M { ; newton's method applied to P(z) = (z + 1) * (z - 0.5 - a) * (z - 0.5 + a) ; Parameter-plane (M-set) fractal ; F(z) = z - P(z)/P'(z) ; Attractors: the 3 zeros of P(z): -1, 0.5+a, 0.5-a ; Critical points: 0, as well the 3 zeros of P(z). init: ; the attractors complex attract[3] int n_attract = 3 complex cf1 ; pre-calculate some factors outside the loop complex cf2 z = @z_init complex zprev = @z_init ; the previous z if @pchoice=="distance" attract[0] = -1 attract[1] = (1-#pixel)/2 attract[2] = (1+#pixel)/2 cf1 = -(3+Sqr(#pixel))/4 cf2 = (Sqr(#pixel) - 1)/4 elseif @pchoice=="halfdist" attract[0] = -1 attract[1] = 0.5 - #pixel attract[2] = 0.5 + #pixel cf1 = -0.75 - Sqr(#pixel) cf2 = Sqr(#pixel) - 0.25 elseif @pchoice=="square" attract[0] = -1 complex tmp = sqrt(#pixel)/2 attract[1] = 0.5 - tmp attract[2] = 0.5 + tmp cf1 = -(3+#pixel)/4 cf2 = (#pixel - 1)/4 else ; @pChoice=="zero" attract[0] = -1 attract[1] = 1 attract[2] = #pixel cf1 = -2 * #pixel cf2 = -#pixel z = #pixel/3 endif int step=0 ; which attractor are we checking ? (with @stripes) int found=-1 ; did we find one? loop: complex dist ; distance to attractor if (@stripes) ; stripes: check only one attractor dist = z-attract[step] if |dist| < @delta found=step else step = step+1 if step>=n_attract step = 0 endif endif else ; no stripes: check all attractors int i=0 repeat dist = z-attract[i] if |dist| < @delta found = i endif i=i+1 until i>=n_attract || found>=0 endif ; stripes if (found<0) && (!@stripes || step==0) ; here we step z (this is the core of the formula) zprev = z ; z[new] = z - P(z)/P'(z) = (z*P'(z) - P(z)) / P'(z) if @pChoice=="zero" z = ((2 * z + cf2) * z^2 + cf2) / (( 3 * z + cf1 ) * z - 1) else z = ( 2 * z^3 + cf2 ) / ( 3 * z^2 + cf1 ) endif endif ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm) if @zmode != "z" if found>=0 if @zmode=="distance" z = dist elseif @zmode=="step" z = z-zprev endif ; zmode endif ; found endif ; zmode bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue found<0 default: title = "Newton-Mandel, degree-3 polynomial" ; v1.3 rating = recommended periodicity = 0 ; turn off periodicity checking by default as it interferes with 'stripes' magn = 0.5 ; default zoom factor. heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="distance" "halfdist" "square" "zero" hint="Choice of a parameter plane. \ 'distance' means #pixel is the distance between the two zeros. The third zero is -1. \ 'halfdist' means #pixel is half the distance. \ 'square' means sqrt(#pixel) is the distance. \ 'zero' means #pixel is the location of a zero, the other two zeros are -1 and 1." default=0 endparam complex param z_init caption="initial z" default = 0 hint = "Where to start the iteration. The default is to start at the critical point 0+0i." endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="How close can z get to a fixed point? \ Actually, this is the square of the distance. \ Smaller values mean more iterations." default=1e-16 min=0 endparam heading caption = "Coloring" endheading bool param stripes default=false hint="Checking this box means you can use one of the striped coloring methods. Use 3 stripes." endparam param zmode caption = "final z" enum = "z" "distance" "step" default = 1 hint="'distance' replaces the final z value with the distance from the (known) attractor. \ 'step' replaces the final z value with the distance from the previous z. \ Since Newton's method converges quadratically, these choices go towards zero in a quadratic way, \ therefore you can use a smooth coloring method (Smooth Mandel / Smooth Stripes / ...) \ with power=2." endparam switch: type = "N_poly3_J" delta = delta pchoice = pchoice a = #pixel stripes = stripes zmode = zmode } # --- Polynomials of degree 4 --------------------------------------- N_Poly4_J { ; newton's method applied to P(z) = (z + 1) * (z - 1) * (z - b + a) * (z - b - a) ; Complex plane (J-set) for specific 'a' and 'b' ; F(z) = z - P(z)/P'(z) ; F(z) = (3*z^4 - 4*b*z^3 + (b^2-a^2-1)*z^2 + (b^2-a^2)) / ( 4*z^3 - 6*b*z^2 + 2*(b^2-a^2-1)*z + 2*b) ; Attractors: the 4 zeros of P(z): 1,-1, b+a, b-a ; Critical points: the 4 zeros, as well as ; b/2 +/- sqrt(3)*sqrt(b^2 + 2*a^2 + 2)/6 ; ; For b=0, this reverts to the 'Newt_fang_xxx' formulas from pwc_convert.ufm global: ; precalculate some factors outside the loop complex fac0 complex fac2 ; turning 'reorder' on will switch attractors 1 and 2 (of 0..3) int ixA int ixB if @reorder ixA = 1 ixB = 2 else ixA = 2 ixB = 1 endif ; the attractors complex attract[4] int n_attract = 4 attract[0] = -1 attract[ixA] = 1 if @pchoice=="halfdist" attract[ixB] = @b-@a attract[3] = @b+@a fac0 = @b^2 - @a^2 fac2 = fac0 - 1 else ; @pchoice=="square" complex tmp = sqrt(-@a) attract[ixB] = @b + tmp attract[3] = @b - tmp fac0 = @b^2 + @a fac2 = fac0 - 1 endif init: z=#pixel complex zprev = z int step=0 ; which attractor are we checking ? (with @stripes) int found=-1 ; did we find one? loop: complex dist ; distance to attractor if (@stripes) ; stripes: check only one attractor dist = z-attract[step] if |dist| < @delta found=step else step = step+1 if step>=n_attract step = 0 endif endif else ; no stripes: check all attractors int i=0 repeat dist = z-attract[i] if |dist| < @delta found = i endif i=i+1 until i>=n_attract || found>=0 endif ; stripes if (found<0) && (!@stripes || step==0) ; here we step z (this is the core of the formula) zprev = z z = (3*z^4 - 4*@b*z^3 + fac2*z^2 + fac0) / ( 4*z^3 - 6*@b*z^2 + 2*fac2*z + 2*@b) endif ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm) if @zmode != "z" if found>=0 if @zmode=="distance" z = dist elseif @zmode=="step" z = z-zprev endif ; zmode endif ; found endif ; zmode bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue found<0 default: title = "Newton-Julia, degree-4 polynomial" ; v1.0 periodicity = 0 ; turn off periodicity checking by default as it interferes with 'stripes' heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="halfdist" "square" default=0 endparam complex param a caption="a" default=exp( (0,1) * #pi * (sqrt(5) - 1) ) endparam complex param b caption="b" default=0 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="How close can z get to a fixed point? \ Actually, this is the square of the distance. \ Smaller values mean more iterations." default=1e-10 min=0 endparam heading caption = "Coloring" endheading bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 4 stripes." default=false endparam bool param reorder ; enabled = @stripes default = false caption = "Reorder attractors" hint = "Exchange the attractors 1 and 2 (of 0..3). \ Matters if you use staggered iteration, or in the rare case when two attractors coincide." endparam param zmode caption = "final z" enum = "z" "distance" "step" default = 1 endparam switch: type = "N_Poly4_Ma" pchoice = pchoice ; a = a b = b delta = delta reorder = reorder stripes = stripes zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_Poly4_Ma { ; newton's method applied to P(z) = (z + 1) * (z - 1) * (z - b + a) * (z - b - a) ; Parameter plane (M-set) for all 'a' and a specific 'b' ; F(z) = z - P(z)/P'(z) ; F(z) = (3*z^4 - 4*b*z^3 + (b^2-a^2-1)*z^2 + (b^2-a^2)) / ( 4*z^3 - 6*b*z^2 + 2*(b^2-a^2-1)*z + 2*b) ; Attractors: the 4 zeros of P(z): 1,-1, b+a, b-a ; Critical points: the 4 zeros, as well as ; b/2 +/- sqrt(3)*sqrt(b^2 + 2*a^2 + 2)/6 ; For b=0, this reverts to the 'Newt_fang_xxx' formulas from pwc_convert.ufm global: int n_attract = 4 ; number of attractors int n_orbits ; number of orbits that have to be checked if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif ; turning 'reorder' on will switch attractors 1 and 2 (of 0..3) int ixA int ixB if @reorder ixA = 1 ixB = 2 else ixA = 2 ixB = 1 endif init: complex a = #pixel complex attract[4] ; n_attract complex orbit[2] ; n_orbits attract[0] = -1 attract[ixA] = 1 ; pre-calculate some factors outside the loop complex tmps complex fac0 complex fac2 if @pchoice=="halfdist" attract[ixB] = @b-a attract[3] = @b+a fac0 = @b^2 - a^2 fac2 = fac0 - 1 tmps = @b^2 + 2*a^2 + 2 else ; @pchoice=="square" complex tmp = sqrt(-a) attract[ixB] = @b + tmp attract[3] = @b - tmp fac0 = @b^2 + a fac2 = fac0 - 1 tmps = @b^2 - 2*a + 2 endif if @giveup == "custom orbit start" orbit[0] = @z_init else tmps = sqrt(tmps/3) / 2 \$ifdef UNUSED if real(tmps)*imag(@b) > -imag(tmps)*real(@b) tmps = - tmps endif \$endif orbit[0] = @b/2 - tmps orbit[1] = @b/2 + tmps if @giveup == "other" orbit[0] = orbit[1] ; orbit[1] = @b/2 - tmps ; not used in this case endif endif z = orbit[0]-attract[0] complex zprev = z ; give it some value to avoid compiler warnings complex dist = 0 ; ditto bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live ; how many orbits can bail out until we're done? if @giveup=="all done" n_live = n_orbits else n_live = 1 endif int step=0 ; which attractor are we checking ? (with @stripes) loop: int ix = 0 ; count through the orbits repeat int found=-1 if continue[ix] ; orbit still live ? z = orbit[ix] if (@stripes) ; stripes: check only one attractor dist = z-attract[step] if |dist| < @delta found=step endif else ; no stripes: check all attractors int j=0 ; loop through the attractors repeat dist = z-attract[j] if |dist| < @delta found = j endif j=j+1 until j>=n_attract || found>=0 endif ; stripes if (found>=0) continue[ix] = false n_live = n_live - 1 elseif (!@stripes || step==n_attract-1) ; here we step z (this is the core of the formula) zprev = z z = (3*z^4 - 4*@b*z^3 + fac2*z^2 + fac0) / ( 4*z^3 - 6*@b*z^2 + 2*fac2*z + 2*@b) orbit[ix] = z endif ; found endif ; continue ix = ix + 1 until ix >= n_orbits || n_live <= 0 if @stripes step = step+1 if step >= n_attract step = 0 endif endif ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm) if @zmode != "z" if n_live <= 0 if @zmode=="distance" z = dist elseif @zmode=="step" z = z-zprev endif ; zmode endif ; found endif ; zmode bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Newton-Mandel(a), degree-4 polynomial" ; v1.0 periodicity = 0 ; turn off periodicity checking by default as it interferes with 'stripes' heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="halfdist" "square" default=0 endparam complex param b caption="b" default=0 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="How close can z get to a fixed point? \ Actually, this is the square of the distance. \ Smaller values mean more iterations." default=1e-10 min=0 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = 0 endparam heading caption = "Coloring" endheading bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 4 stripes." default=false endparam bool param reorder ; enabled = @stripes default = false caption = "Reorder attractors" hint = "Exchange the attractors 1 and 2 (of 0..3). \ Matters if you use staggered iteration, or in the rare case when two attractors coincide." endparam param zmode caption = "final z" enum = "z" "distance" "step" default = 1 ; "distance" endparam switch: type = "N_Poly4_J" pchoice = pchoice a = #pixel b = b delta = delta reorder = reorder stripes = stripes zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_Poly4_Mb { ; newton's method applied to P(z) = (z + 1) * (z - 1) * (z - b + a) * (z - b - a) ; Parameter plane (M-set) for a specific 'a' and all 'b' ; F(z) = z - P(z)/P'(z) ; F(z) = (3*z^4 - 4*b*z^3 + (b^2-a^2-1)*z^2 + (b^2-a^2)) / ( 4*z^3 - 6*b*z^2 + 2*(b^2-a^2-1)*z + 2*b) ; Attractors: the 4 zeros of P(z): 1,-1, b+a, b-a ; Critical points: the 4 zeros, as well as ; b/2 +/- sqrt(3)*sqrt(b^2 + 2*a^2 + 2)/6 ; v1.1 2004-01-16 Added yipe_bug param global: int n_attract = 4 ; number of attractors int n_orbits ; number of orbits that have to be checked if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif ; pre-calculate some intermediate outside the loop for efficiency complex pre_fac0 complex pre_tmps complex plusminus if @pchoice=="halfdist" plusminus = @a pre_fac0 = - (@a^2) pre_tmps = 2 + 2*@a^2 else ; @pchoice=="square" plusminus = sqrt(-@a) pre_fac0 = @a pre_tmps = 2 - 2*@a endif ; turning 'reorder' on will switch attractors 1 and 2 (of 0..3) int ixA int ixB if @reorder ixA = 1 ixB = 2 else ixA = 2 ixB = 1 endif init: complex b = #pixel complex attract[4] ; n_attract complex orbit[2] ; n_orbits attract[0] = -1 attract[ixA] = 1 attract[ixB] = b-plusminus attract[3] = b+plusminus ; pre-calculate some factors outside the loop complex fac0 complex fac2 complex tmps fac0 = b^2 + pre_fac0 fac2 = fac0 - 1 if @giveup == "custom orbit start" orbit[0] = @z_init else tmps = sqrt( (pre_tmps + b^2) / 3 ) / 2 if @yipe_bug ; v1.1 if real(tmps)*imag(b) > imag(tmps)*real(b) tmps = - tmps endif else ; bug-less if real(tmps)*real(b) > -imag(tmps)*imag(b) tmps = - tmps endif endif orbit[0] = b/2 - tmps orbit[1] = b/2 + tmps if @giveup == "other" orbit[0] = orbit[1] ; orbit[1] = b/2 - tmps ; not used in this case endif endif z = orbit[0]-attract[0] complex zprev = z ; avoid warnings complex dist = 0 ; ditto bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live ; how many orbits can bail out until we're done? if @giveup=="all done" n_live = n_orbits else n_live = 1 endif int step=0 ; which attractor are we checking ? (with @stripes) loop: int ix = 0 ; count through the orbits repeat int found=-1 if continue[ix] ; orbit still live? z = orbit[ix] if (@stripes) ; stripes: check only one attractor dist = z-attract[step] if |dist| < @delta found=step endif ; dist else ; no stripes: check all attractors int j=0 ; loop through the attractors repeat dist = z-attract[j] if |dist| < @delta found = j endif j=j+1 until j>=n_attract || found>=0 endif ; stripes if (found>=0) continue[ix] = false n_live = n_live - 1 elseif (!@stripes || step==n_attract-1) ; here we step z (this is the core of the formula) zprev = z z = (3*z^4 - 4*b*z^3 + fac2*z^2 + fac0) / ( 4*z^3 - 6*b*z^2 + 2*fac2*z + 2*b) orbit[ix] = z endif ; found endif ; continue ix = ix + 1 until ix >= n_orbits || n_live <= 0 if @stripes step = step+1 if step >= n_attract step = 0 endif endif ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm) if @zmode != "z" if n_live <= 0 if @zmode=="distance" z = dist elseif @zmode=="step" z = z-zprev endif ; zmode endif ; found endif ; zmode bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Newton-Mandel(b), degree-4 polynomial" ; v1.1 periodicity = 0 ; turn off periodicity checking by default as it interferes with 'stripes' heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="halfdist" "square" default=0 endparam complex param a caption="a" default=exp( (0,1) * #pi * (sqrt(5) - 1) ) endparam bool param yipe_bug caption="reproduce 'yipe' bug?" hint="checking this options reproduces a bug which gives a funny appearance \ for a = 1 [square], \ respectively a = 0+1i [halfdist]" default = false endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="How close can z get to a fixed point? \ Actually, this is the square of the distance. \ Smaller values mean more iterations." default=1e-10 min=0 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = 0 endparam heading caption = "Coloring" endheading bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 4 stripes." default=false endparam bool param reorder ; enabled = @stripes default = false caption = "Reorder attractors" hint = "Exchange the attractors 1 and 2 (of 0..3). \ Matters if you use staggered iteration, or in the rare case when two attractors coincide." endparam param zmode caption = "final z" enum = "z" "distance" "step" default = 1 ; "distance" endparam switch: type = "N_Poly4_J" pchoice = pchoice a = a b = #pixel delta = delta reorder = reorder stripes = stripes zmode = zmode } # --- Polynomials of degree k (k arbitrary) ------------------------- N_polyK_J { ; Newton's method applied to (a subset of the) polynomials of degree k ; ; newton's method for P(z) = z^k + p*z + q ; Complex plane (J-set) for specific 'p' and 'q' ; F(z) = z - P(z)/P'(z) ; F(z) = ((k-1)*z^k - q)/(k*z^(k-1) + p) ; Attractors: the zeroes of P(z) ; Critical points: the zeroes of P(z), as well as the origin (z=0) init: z=#pixel complex delta_z bool continue = true ; bailout loop: if @zdirect ; the result should be the same with or without, except for rounding error complex zprev = z z = ( (@k-1)*(z^@k) - @q ) / ( @k*z^(@k-1) + @p ) delta_z = z-zprev else delta_z = - (z^@k + @p*z + @q) / (@k*z^(@k-1) + @p) z = z + delta_z endif if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z elseif @zmode=="f" z = z^@k + @p*z + @q endif ; zmode endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Julia, degree-k polynomial" ; v1.0 periodicity = 0 ; turn off periodicity checking by default magn = 0.25 ; default zoom factor. heading caption = "Fractal shape" endheading int param k caption="degree k" hint="For positive k, this is the degree (or order) of the polynomial P(z) whose zeroes we're seeking \ with Newton's method. \ 'k' should be greater than 2 or smaller than -1. \ Actually, 2 and -1 work, but are trivial." default=5 endparam complex param p caption="p" hint="The coefficient p in the function P(z) := z^k + p*z + q (We're seeking the zeroes of this function using Newton iteration.)" default=-1 endparam complex param q caption="q" hint="The coefficient q in the function P(z) := z^k + p*z + q (We're seeking the zeroes of this function using Newton iteration.)" default=1 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="If (the square of) the change in z is smaller than this value, bail out. \ Smaller values mean more iterations." default=1e-8 min=0 endparam bool param zdirect caption = "direct" hint = "If this checkbox does anything, then only because of roundoff errors." default = true endparam heading caption = "Coloring" endheading param zmode caption = "final z" enum = "z" "step" "f" default = 1 hint="'step' replaces the final z value with 'z-z[prev]'. 'f' replaces the final z value with f(z). \ Since Newton's method converges quadratically, these choices go towards zero in a quadratic way, \ therefore you can use a smooth coloring method (Smooth Mandel / Smooth Stripes / ...)" endparam switch: type = "N_polyK_Mp" k = k p = p q = q delta = delta zdirect = zdirect zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_polyK_Mq { ; Newton's method applied to (a subset of the) polynomials of degree k ; ; newton's method for P(z) = z^k + p*z + q ; Parameter plane (M-set) for a specific 'p' and all 'q' ; F(z) = z - P(z)/P'(z) ; F(z) = ((k-1)*z^k - q)/(k*z^(k-1) + p) ; Attractors: the zeroes of P(z) ; Critical points: the zeroes of P(z), as well as the origin (z=0) init: z=0 complex delta_z bool continue = true loop: if @zdirect ; the result should be the same with or without, except for rounding error complex zprev = z z = ( (@k-1)*(z^@k) - #pixel ) / ( @k*z^(@k-1) + @p ) delta_z = z-zprev else delta_z = - (z^@k + @p*z + #pixel) / ( @k*z^(@k-1) + @p ) z = z + delta_z endif if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z elseif @zmode=="f" z = z^@k + @p*z + #pixel endif ; zmode endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Mandel(q), degree-k polynomial" ; v1.0 periodicity = 0 ; turn off periodicity checking by default heading caption = "Fractal shape" endheading int param k caption="degree k" hint="For positive k, this is the degree (or order) of the polynomial P(z) whose zeroes we're seeking \ with Newton's method. \ 'k' should be greater than 2 or smaller than -1." default=5 endparam complex param p caption="p" hint="The coefficient p in the function P(z) := z^k + p*z + q (We're seeking the zeroes of this function using Newton iteration.)" default=-1 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="If (the square of) the change in z is smaller than this value, bail out. \ Smaller values mean more iterations." default=1e-8 min=0 endparam bool param zdirect caption = "direct" hint = "If this checkbox does anything, then only because of roundoff errors." default = true endparam heading caption = "Coloring" endheading param zmode caption = "final z" enum = "z" "step" "f" default = 1 hint="'step' replaces the final z value with 'z-z[prev]'. 'f' replaces the final z value with f(z). \ Since Newton's method converges quadratically, these choices go towards zero in a quadratic way, \ therefore you can use a smooth coloring method (Smooth Mandel / Smooth Stripes / ...)" endparam switch: type = "N_polyK_J" k = k p = p q = #pixel delta = delta zdirect = zdirect zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_polyK_Mp { ; Newton's method applied to (a subset of the) polynomials of degree k ; ; newton's method for P(z) = z^k + p*z + q ; Parameter plane (M-set) for all 'p' and a specific'q' ; F(z) = z - P(z)/P'(z) ; F(z) = ((k-1)*z^k - q)/(k*z^(k-1) + p) ; Attractors: the zeroes of P(z) ; Critical points: the zeroes of P(z), as well as the origin (z=0) init: z=0 complex delta_z bool continue = true loop: if @zdirect ; the result should be the same with or without, except for rounding error complex zprev = z z = ( (@k-1)*(z^@k) - @q ) / ( @k*z^(@k-1) + #pixel ) delta_z = z-zprev else delta_z = - (z^@k + #pixel*z + @q) / (@k*z^(@k-1) + #pixel) z = z + delta_z endif if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z elseif @zmode=="f" z = z^@k + #pixel*z + @q endif ; zmode endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Mandel(p), degree-k polynomial" ; v1.0 periodicity = 0 heading caption = "Fractal shape" endheading int param k caption="degree k" hint="For positive k, this is the degree (or order) of the polynomial P(z) whose zeroes we're seeking \ with Newton's method. \ 'k' should be greater than 2 or smaller than -1." default=5 endparam complex param q caption="q" hint="The coefficient q in the function P(z) := z^k + p*z + q (We're seeking the zeroes of this function using Newton iteration.)" default=1 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="If (the square of) the change in z is smaller than this value, bail out. \ Smaller values mean more iterations." default=1e-8 min=0 endparam bool param zdirect caption = "direct" hint = "If this checkbox does anything, then only because of roundoff errors." default = true endparam heading caption = "Coloring" endheading param zmode caption = "final z" enum = "z" "step" "f" default = 1 hint="'step' replaces the final z value with 'z-z[prev]'. 'f' replaces the final z value with f(z). \ Since Newton's method converges quadratically, these choices go towards zero in a quadratic way, \ therefore you can use a smooth coloring method (Smooth Mandel / Smooth Stripes / ...)" endparam switch: type = "N_polyK_J" k = k p = #pixel q = q delta = delta zdirect = zdirect zmode = zmode } # --- Sine ---------------------------------------------------------- N_sin_J { ; newton's method applied to f(z) = sin(z) - c ; Q(z) = z - f(z)/f'(z) = z - ((sin(z)-c)/cos(z)) global: init: z = #pixel complex delta_z = 0 bool continue = true loop: delta_z = ( sin(z) - @cv ) / cos(z) if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z endif ; zmode else z = z - delta_z endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Julia, sin(z)-c" ; v1.0 rating = recommended magn = 0.25 periodicity = 0 heading caption = "Fractal shape" endheading complex param cv caption="c" default = (0, 0.5) endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="Square of max step. Use small values." default=1e-8 min=0 endparam param zmode caption = "final z" enum = "z" "step" default = 1 endparam switch: type = "N_sin_M" delta = delta zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_sin_M { ; newton's method applied to f(z) = sin(z) - c ; Q(z) = z - f(z)/f'(z) = z - ((sin(z)-c)/cos(z)) global: init: z = 0 complex delta_z = 0 complex c = #pixel bool continue = true loop: delta_z = ( sin(z) - c ) / cos(z) if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z endif ; zmode else z = z - delta_z endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-M, sin(z)-c" ; [v1.0] rating = recommended magn = 0.25 periodicity = 0 heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="Square of max step. Use small values." default=1e-8 min=0 endparam param zmode caption = "final z" enum = "z" "step" default = 1 endparam switch: type = "N_sin_J" cv = #pixel delta = delta zmode = zmode } # --- tan ----------------------------------------------------------- N_tan_J { ; newton's method applied to f(z) = tan(z) - c ; Q(z) = z - f(z)/f'(z) = z - cos(z) ( c*cos(z) - sin(z) ) global: init: z = #pixel complex delta_z = 0 bool continue = true loop: complex cosz = cos(z) delta_z = cosz * ( @cv * cosz - sin(z) ) if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z endif ; zmode else z = z + delta_z endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Julia, tan(z)-c" ; [v1.0] magn = 0.25 periodicity = 0 heading caption = "Fractal shape" endheading complex param cv caption="c" default = (0, 0.5) endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="Square of max step. Use small values." default=1e-8 min=0 endparam param zmode caption = "final z" enum = "z" "step" default = 1 endparam switch: type = "N_tan_M" delta = delta zmode = zmode } # - - - - - - - - - - - - - - - - - - - - - - - N_tan_M { ; newton's method applied to f(z) = tan(z) - c ; Q(z) = z - f(z)/f'(z) = z - cos(z) * ( c*cos(z) - sin(z) ) global: init: complex delta_z = 0 complex c = #pixel bool continue = true z = #pi * @crit loop: complex cosz = cos(z) delta_z = cosz * ( c * cosz - sin(z) ) if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z endif ; zmode else z = z + delta_z endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-M, tan(z)-c" ; v1.0 magn = 0.25 periodicity = 0 int param crit default = 0 endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="Square of max step. Use small values." default=1e-8 min=0 endparam param zmode caption = "final z" enum = "z" "step" default = 1 endparam switch: type = "N_tan_J" cv = #pixel delta = delta zmode = zmode } # --- Exp ----------------------------------------------------------- N_exp_J { ; newton's method applied to f(z) = exp(z) - c ; Q(z) = z - f(z)/f'(z) = z - ((exp(z)-c)/exp(z)) = z - 1 + c * exp(-z) global: init: z = #pixel complex delta_z = 0 bool continue = true loop: delta_z = @cv * exp(-z) - 1 if |delta_z| < @delta continue = false if @zmode=="step" z = delta_z endif ; zmode else z = z + delta_z endif ; < @delta bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Newton-Julia of exp(z)-c v1.0" magn = 0.25 periodicity = 0 heading caption = "Fractal shape" endheading complex param cv caption="c" default = (0, 0.5) endparam heading caption = "Bailout" endheading float param delta caption="Bailout delta" hint="Square of max step. Use small values." default=1e-8 min=0 endparam param zmode caption = "final z" enum = "z" "step" default = 1 endparam } ##################################################################### ### Analytic functions (not Newton iteration) ### ##################################################################### # --- hring --------------------------------------------------------- HRing_J { ; Iteration of f(z) = a * z^@power * (z - c) / (1 - c*z) and related formulas ; ; This thing is capable of generating Herman rings for alpha ; equal to exp(2*pi*i*a), a irrational, power k=2, and large c (try 4.0) ; ; ; Sorry, lost the author of the original formula (picked it up ages ago in a newsgroup) ; the formula is the same as HRING.FRM in the orgform collection. ; ; Adapted for UltraFractal by A. Hoppler ; ; 2004-02-01 v1.4 added explicit outer bailout ; 2007-06-27 v1.5 added variants, tweaked calculations global: ; bailout radii float rad_inner = @delta float rad_outer = @bailout if @bailout <= 0 if @delta <= 0 rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif @delta <= 0 rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif ; parameters complex cc = @cp complex dd = @cp if @conjugate dd = conj(dd) endif complex lambda if @pchoice=="angle" lambda = exp( #pi * (0,2) * @alpha ) else lambda = @alpha endif \$ifdef DEBUG print( "lambda=", lambda) \$endif init: #z=#pixel bool continue = true int which = 0 ; Which bailout to check (if @stripes) float rad = 1 ; initialize, so the compiler doesn't complain loop: if (!@stripes) || which==0 ; calculate next z if @power==2 #z = lambda * sqr(#z) * (#z - cc) / (1 - dd*#z) elseif @hipower #z = lambda * #z^@power * ((#z - cc) / (1 - dd*#z))^(@power-1) else #z = lambda * #z^@power * (#z - cc) / (1 - dd*#z) endif rad = |#z| endif if @stripes ; check one bailout which = which+1 if which==1 continue = rad >= rad_inner else which = 0 continue = rad <= rad_outer endif else ; check both bailouts continue = rad >= rad_inner && rad < rad_outer endif bailout: ; false = bail; true = continue continue default: title = "Herman ring J" ; v1.5 rating = recommended magn = 0.5 center = (0,0) periodicity = 0 heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" hint="lambda = exp(2*#pi*i*alpha)" enum="angle" "lambda" default=0 endparam complex param alpha caption="alpha" hint="Angle alpha or factor lambda, depending on param interpretation" default = (Sqrt(5)-1)/2 endparam complex param cp caption="c" default= (4, 0) endparam bool param conjugate default=false endparam int param power caption="Power k" default=2 endparam bool param hipower caption = "High-power variant" endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "HRing_Mc" alpha = alpha pchoice = pchoice center = cp conjugate = conjugate power = power hipower = hipower delta = delta bailout = bailout stripes = stripes } # - - - - - - - - - - - - - - - - - - - - - - - HRing_Mc { ; Parameter-Plane (M-Type) fractal of the hring formula ; ; Iteration of f(z) = a * z^@power * (z - c) / (1 - c*z), with c=#pixel ; or a related formula ; ; Sorry, lost the author of the original formula (picked it up ages ago in a newsgroup) ; the formula is the same as HRING.FRM in the orgform collection. ; ; Adapted for UltraFractal by A. Hoppler ; ;; 2004-02-01 v1.4 added explicit outer bailout, added custom orbit ;; 2006-06-27 v1.5 added hipower variant, tweaked calculations global: ; bailout radii float rad_inner = @delta float rad_outer = @bailout if @bailout <= 0 if @delta <= 0 rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif @delta <= 0 rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 ; calculate all two orbits else n_orbits = 1 ; one orbit only endif complex lambda if @pchoice=="angle" lambda = exp( #pi * (0,2) * @alpha ) else lambda = @alpha endif init: complex c = #pixel complex d = c if (@conjugate) d = conj(c) endif complex orbit[2] ; n_orbits if @giveup == "custom orbit start" orbit[0] = @z_init else if @hipower ; calculate some intermediate values complex t1cd = 1-c*d complex t22 = t1cd*(4*@power^2 - 4*@power + t1cd) complex t1 = 2*@power - t1cd complex t2 = sqrt(t22) complex t3 = 2*@power*d else ; calculate some intermediate values complex dc_p_1 = d*c + 1 complex dc_m_1 = d*c - 1 complex dc2 = sqr(dc_m_1) complex t2 = sqrt( dc2*(@power^2+1) -2*(dc_p_1)*(dc_m_1)*@power ) complex t1 = dc_p_1*@power - dc_m_1 complex t3 = 2*@power*d endif orbit[0] = ( t1+t2 ) / t3 ; critical point orbit[1] = ( t1-t2 ) / t3 ; other critical point ; third critical point: 0 - always superattractive ; fourth critical point: infinity - always superattractive if @giveup == "other" orbit[0] = orbit[1] elseif @giveup == "adjusted" || @giveup == "adjusted 2" bool flip = @giveup != "adjusted" if imag(t1) < 0 flip = ! flip endif if imag(t2) > 0 flip = ! flip endif if flip orbit[0] = orbit[1] endif endif endif ; set to initial values so the compiler doesn't complain #z = orbit[0] float rad = |#z| bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live ; how many orbits must die until we bail out? if @giveup=="all done" n_live = 2 else n_live = 1 endif bool test_in = false ; test in or out (for stripes) loop: int ix = 0 repeat if continue[ix] #z = orbit[ix] float rad = |#z| bool bail if @stripes ; stripes: check only one attractor if test_in bail = (rad < rad_inner) else bail = (rad > rad_outer) endif else ; no stripes: check inner and outer bail = (rad < rad_inner) || (rad > rad_outer) endif ; @stripes if bail continue[ix] = false n_live = n_live - 1 else if (! @stripes) || test_in if @power==2 orbit[ix] = lambda * sqr(#z) * (#z - c) / (1 - d*#z) elseif @hipower orbit[ix] = lambda * #z^@power * ((#z - c) / (1 - d*#z))^(@power-1) else orbit[ix] = lambda * #z^@power * (#z - c) / (1 - d*#z) endif endif endif ; bail endif ; continue ix = ix+1 until ix >= n_orbits || n_live <= 0 if (@stripes) test_in = ! test_in endif bailout: ; false = bail; true = continue n_live > 0 default: title = "Herman ring M(c)" ; v1.4 rating = recommended center = (0, 0) magn = 0.4 periodicity = 0 maxiter = 1000 heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="angle" "lambda" default=0 endparam complex param alpha caption="alpha" hint="Angle alpha. Should be irrational." default = (Sqrt(5)-1)/2 endparam bool param conjugate default=false endparam int param power caption="Power k" default=2 endparam bool param hipower caption = "High-power variant" endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "adjusted" "adjusted 2" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = (1.0, 0.0) endparam bool param stripes default = false endparam switch: type = "HRing_J" alpha = alpha pchoice = pchoice cp = #pixel conjugate = conjugate power = power hipower = hipower delta = delta bailout = bailout stripes = stripes } # - - - - - - - - - - - - - - - - - - - - - - - HRing_Ma { ; Parameter-plane (M-Type) fractal of the hring formula. ; ; Iteration of f(z) = a * z^@power * (z - c) / (1 - c*z), with a=#pixel ; or a similar formula ; ; Sorry, lost the author of the original formula (picked it up ages ago in a newsgroup) ; the formula is the same as HRING.FRM in the orgform collection. ; ; Adapted for UltraFractal by A. Hoppler ; ;; 2004-02-01 v1.4 added explicit outer bailout, added custom orbit ;; 2006-06-27 v1.5 added hipower variant, tweaked calculations global: ; bailout radii float rad_inner = @delta float rad_outer = @bailout if @bailout <= 0 if @delta <= 0 rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif @delta <= 0 rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 ; calculate all two orbits else n_orbits = 1 ; one orbit only endif complex c = @cp complex d = @cp if (@conjugate) d = conj(d) endif complex orbit_init[2] ; n_orbits if @giveup == "custom orbit start" orbit_init[0] = orbit_init[1] = @z_init else if @hipower ; calculate some intermediate values complex t1cd = 1-c*d complex t22 = t1cd*(4*@power^2 - 4*@power + t1cd) complex t1 = 2*@power - t1cd complex t2 = sqrt(t22) complex t3 = 2*@power*d else ; calculate some intermediate values complex dc_p_1 = d*c + 1 complex dc_m_1 = d*c - 1 complex dc2 = sqr(dc_m_1) complex t2 = sqrt( dc2*(@power^2+1) -2*(dc_p_1)*(dc_m_1)*@power ) complex t1 = dc_p_1*@power - dc_m_1 complex t3 = 2*@power*d endif orbit_init[0] = ( t1+t2 ) / t3 ; critical point orbit_init[1] = ( t1-t2 ) / t3 ; other critical point ; third critical point: 0 - always superattractive ; fourth critical point: c - only if @hipower - always mapped to 0 if @giveup == "other" orbit_init[0] = orbit_init[1] endif endif init: complex lambda if @pchoice=="angle" lambda = exp( #pi * (0,2) * #pixel ) else lambda = #pixel endif complex orbit[2] ; n_orbits orbit = orbit_init #z = orbit[0] bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live if @giveup==0 n_live = n_orbits else n_live = 1 endif bool test_in = false ; test in or out (for stripes) loop: int ix = 0 repeat if continue[ix] #z = orbit[ix] float rad = |#z| bool bail if @stripes ; stripes: check only one attractor if test_in bail = (rad < rad_inner) else bail = (rad > rad_outer) endif else ; no stripes: check inner and outer bail = (rad < rad_inner) || (rad > rad_outer) endif ; @stripes if bail continue[ix] = false n_live = n_live - 1 else if (! @stripes) || test_in if @power==2 orbit[ix] = lambda * sqr(#z) * (#z - c) / (1 - d*#z) elseif @hipower orbit[ix] = lambda * #z^@power * ((#z - c) / (1 - d*#z))^(@power-1) else orbit[ix] = lambda * #z^@power * (#z - c) / (1 - d*#z) endif endif endif ; bail endif ; continue ix = ix+1 until ix >= n_orbits || n_live <= 0 if (@stripes) test_in = ! test_in endif bailout: ; false = bail; true = continue n_live > 0 default: title = "Herman ring M(a)" ; v1.5 rating = recommended center = (0, 0) magn = 2 periodicity = 0 maxiter = 1000 heading caption = "Fractal shape" endheading param pchoice caption="param interpretation" enum="angle" "lambda" default=0 endparam complex param cp caption="c" default = (0,9) endparam bool param conjugate default=false endparam int param power caption="Power k" default=2 endparam bool param hipower caption = "High-power variant" endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = (1.0, 0.0) endparam bool param stripes default = false endparam switch: type = "HRing_J" alpha = #pixel pchoice = pchoice cp = cp conjugate = conjugate power = power hipower = hipower bailout = bailout delta = delta stripes = stripes } # --- ringoes (inspired by 'hring') --------------------------------- Ringo23_J { ; ; Iteration of f(z) = z^@k * (z^2 + a*z + b) / (c*z + 1) ; ; inspired by 'hring' formula ; ; 2004-02-02 v1.4 separated inner and outer bailout global: float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif init: z=#pixel bool continue = true int which=0 float rad = 1 ; just so the compiler doesn't complain loop: if (!@stripes) || which==0 z = z^@k * (z^2 + @a*z + @b) / (@c*z + 1) rad = |z| endif if @stripes which = which+1 if which==1 continue = rad >= rad_inner else which = 0 continue = rad <= rad_outer endif else continue = rad >= rad_inner && rad < rad_outer endif bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Ringo 2/3 J" ; v1.4 magn = 0.5 periodicity = 0 heading caption = "Fractal shape" endheading complex param a default= 0 endparam complex param b default= 1 endparam complex param c default= 1 endparam int param k caption="Power" default=2 endparam param plane caption = "switch to plane ..." enum = "a" "b" "c" default = 0 endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "Ringo23_M" a = a b = b c = c k = k plane = plane delta = delta stripes = stripes } # - - - - - - - - - - - - - - - - - - - - - - - Ringo23_M { ; Parameter-Plane (M-Type) fractal ; ; Iteration of f(z) = z^@k * (z^2 + a*z + b) / (c*z + 1) ; ; inspired by 'hring' formula ; ; 2004-02-02 v1.4 separated inner and outer bailout global: float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif complex rt1 = sqrt(3)*(0,-0.5) - 0.5 complex rt2 = sqrt(3)*(0, 0.5) - 0.5 init: complex aa = @a complex bb = @b complex cc = @c if @plane=="a" aa = #pixel elseif @plane=="b" bb = #pixel elseif @plane=="c" cc = #pixel endif complex a3 = cc*(@k+1) complex a2 = (aa*cc+1)*@k + 2 complex a1 = bb*cc*(@k-1) + aa*(@k+1) complex a0 = bb*@k int n_orbits = 3 complex orbit[3] ; n_orbits complex cog complex sqrtD if (a3==0) if (a2==0) n_orbits = 1 orbit[0] = - a0/a1 orbit[1] = orbit[0] orbit[2] = orbit[0] else n_orbits = 2 a1 = a1/(2*a2) a0 = a0/a2 sqrtD = sqrt( a1*a1 - a0 ) orbit[0] = -a1 + sqrtD orbit[1] = -a1 - sqrtD orbit[2] = orbit[1] endif else a2 = a2/a3 a1 = a1/a3 a0 = a0/a3 complex Q = (3*a1 - a2^2)/9 complex R = (9*a2*a1 - 27*a0 - 2*a2^3)/54 complex D = Q^3 + R^2 sqrtD = sqrt(D) complex S = (R+sqrtD)^(1/3) complex T = -Q/S cog = -a2/3 orbit[0] = cog + S + T orbit[1] = cog + rt2*S + rt1*T orbit[2] = cog + rt1*S + rt2*T endif \$ifdef DEBUG if real(#random)<0.05 print( q2, q1, q0 ) int k=0 while k<3 print( k, ":", orbit[k]^3 + q2*orbit[k]^2 + q1*orbit[k] + q0 ) k = k+1 endwhile endif \$endif if @giveup!="all done" && @giveup!="any done" n_orbits = 1 if @giveup == "other" orbit[0] = orbit[1] elseif @giveup == "yet another" orbit[0] = orbit[2] elseif @giveup == "custom orbit start" orbit[0] = @z_init endif endif complex z = orbit[0] bool continue[3] ; n_orbits continue[0] = true continue[1] = true continue[2] = true int n_live if @giveup=="all done" n_live = n_orbits else n_live = 1 endif bool test_in = false loop: int ix = 0 repeat if continue[ix] z = orbit[ix] float rad = |z| bool bail if @stripes if test_in bail = (rad < rad_inner) else bail = (rad > rad_outer) endif else ; not @stripes bail = (rad < rad_inner) || (rad > rad_outer) endif ; @stripes if bail continue[ix] = false n_live = n_live - 1 else if (! @stripes) || test_in orbit[ix] = z^@k * (z^2 + aa*z + bb) / (cc*z + 1) endif endif ; bail / which endif ; continue ix = ix+1 until ix >= n_orbits || n_live <= 0 test_in = ! test_in bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Ringo 2/3 M" ; v1.4 center = (0, 0) magn = 0.4 periodicity = 0 maxiter = 1000 heading caption = "Fractal shape" endheading complex param a default= 0 endparam complex param b default= 1 endparam complex param c default= 1 endparam int param k caption="Power" default=2 endparam param plane caption = "plane" enum = "a" "b" "c" default = 0 endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "yet another" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = (1.0, 0.0) endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "Ringo23_J" plane = plane k = k delta = delta stripes = stripes a = #pixel ; a = a b = b c = c } # - - - - - - - - - - - - - - - - - - - - - - - Ringo_power_J { ; 2004-02-02 v1.4 separated inner and outer bailout global: float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif init: z=#pixel bool continue = true int which = 0 float rad = 1 ; just so the compiler doesn't complain loop: if (!@stripes) || which==0 complex tmp1=@b*z + 1 z = z^@power * (tmp1+@d) / (@b*tmp1) rad = |z| endif if @stripes which = which+1 if which==1 continue = rad >= rad_inner else which = 0 continue = rad <= rad_outer endif else continue = rad >= rad_inner && rad < rad_outer endif bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue continue default: title = "Ringo power J" ; v1.4 magn = 0.5 periodicity = 0 heading caption = "Fractal shape" endheading complex param d caption="d" ; hint="something" default = 1 endparam complex param b caption="b" default= 1 endparam int param power caption="Power k" default=2 endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "ringo_power_Md" b = b ; #center = d power = power delta = delta stripes = stripes } # - - - - - - - - - - - - - - - - - - - - - - - Ringo_power_Mb { ; 2004-02-02 v1.4 separated inner and outer bailout global: float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif complex sol1 = (-@d-2)*@power + @d complex sol2 = SQRT( @d^2*(@power^2 + 1) + (-2*@d^2 - 4*@d)*@power ) init: complex b = #pixel ; now calculate critical points complex orbit[2] ; n_orbits complex sol3 = 2*b*@power if @giveup == "custom orbit start" orbit[0] = orbit[1] = @z_init else orbit[0] = (sol1+sol2) / sol3 ; critical point orbit[1] = (sol1-sol2) / sol3 ; other critical point if @giveup == "other" orbit[0] = orbit[1] elseif @giveup == "adjusted" || @giveup == "adjusted ii" bool test = @giveup != "adjusted" if imag(sol1) < 0 test = !test endif if imag(sol2) < 0 test = !test endif if test orbit[0] = orbit[1] endif endif endif ; @giveup complex z = orbit[0] bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live if @giveup==0 n_live = 2 else n_live = 1 endif bool test_in = false loop: int ix = 0 repeat if continue[ix] z = orbit[ix] float rad = |z| bool bail if @stripes if test_in bail = (rad < rad_inner) else bail = (rad > rad_outer) endif else ; not @stripes bail = (rad < rad_inner) || (rad > rad_outer) endif ; @stripes if bail continue[ix] = false n_live = n_live - 1 else if (! @stripes) || test_in complex tmp1=b*z + 1 orbit[ix] = z^@power * (tmp1+@d) / (b*tmp1) endif endif ; bail / which endif ; continue ix = ix+1 until ix >= n_orbits || n_live <= 0 if (@stripes) test_in = ! test_in endif bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Ringo power M(b)" ; v1.4 center = (0, 0) magn = 0.4 periodicity = 0 maxiter = 1000 heading caption = "Fractal shape" endheading complex param d caption="d" ; hint="something" default = 1 endparam int param power caption="Power k" default=2 endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "adjusted" "adjusted ii" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = (1.0, 0.0) endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "Ringo_power_J" d = d b = #pixel power = power delta = delta stripes = stripes } # - - - - - - - - - - - - - - - - - - - - - - - Ringo_power_Md { ; 2004-02-02 v1.4 separated inner and outer bailout global: float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif complex sol3 = 2*@b*@power init: complex d = #pixel ; now calculate z complex orbit[2] ; n_orbits if @giveup == "custom orbit start" orbit[0] = orbit[1] = @z_init else complex sol1 = (-d-2)*@power + d complex d2 = d^2 complex sol2 = SQRT( d2*(@power^2 + 1) + (-2*d2 - 4*d)*@power ) orbit[0] = (sol1+sol2) / sol3 ; critical point orbit[1] = (sol1-sol2) / sol3 ; other critical point if @giveup == "other" orbit[0] = orbit[1] elseif @giveup == "adjusted" || @giveup == "adjusted ii" bool test = @giveup != "adjusted" if imag(sol1) < 0 test = !test endif if imag(sol2) < 0 test = !test endif if test orbit[0] = orbit[1] endif endif endif complex z = orbit[0] bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live if @giveup==0 n_live = 2 else n_live = 1 endif bool test_in = false loop: int ix = 0 repeat if continue[ix] z = orbit[ix] float rad = |z| bool bail if @stripes if test_in bail = (rad < rad_inner) else bail = (rad > rad_outer) endif else ; not @stripes bail = (rad < rad_inner) || (rad > rad_outer) endif ; @stripes if bail continue[ix] = false n_live = n_live - 1 else if (! @stripes) || test_in complex tmp1=@b*z + 1 orbit[ix] = z^@power * (tmp1+d) / (@b*tmp1) endif endif ; bail / which endif ; continue ix = ix+1 until ix >= n_orbits || n_live <= 0 if (@stripes) test_in = ! test_in endif bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Ringo power M(d)" ; v1.4 center = (0, 0) magn = 0.4 periodicity = 0 maxiter = 1000 heading caption = "Fractal shape" endheading complex param b caption="b" ; hint="something" default = 1 endparam int param power caption="Power k" default=2 endparam heading caption = "Bailout" endheading float param delta caption="small bailout" hint="Orbit trap radius about 0. If <= 0, use reciprocal of large bailout" default=1e-16 endparam float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" "adjusted" "adjusted ii" "custom orbit start" default = 0 endparam complex param z_init caption = "orbit start" enabled = @giveup == "custom orbit start" default = (1.0, 0.0) endparam bool param stripes caption="Staggered iteration" hint="Checking this box means you can use one of the striped coloring methods. Use 2 stripes." default=false endparam switch: type = "Ringo_power_J" d = #pixel b = b power = power delta = delta stripes = stripes } # --- Polynomials of degree 3 --------------------------------------- Poly3_J { ; iteration of f(z) := alpha*z^3 - 3*alpha*z + beta; ; Z-plane (J-set) fractal ; parameters: alpha, beta global: init: z=#pixel loop: z = z * @aa * (z^2 - 3) + @bb; bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue |z| < @bailout default: title = "Polynomial of degree 3" rating = notRecommended periodicity = 0 complex param aa caption="alpha" default=exp( (0,2) * #pi / 3 ) endparam complex param bb caption="beta" default=0 endparam float param bailout caption="Bailout" default=1e10 min=1 endparam switch: type = "Poly3_Ma" bb = bb ; aa = aa bailout = bailout } # - - - - - - - - - - - - - - - - - - - - - - - Poly3_Ma { ; iteration of f(z) := alpha*z^3 - 3*alpha*z + beta; ; alpha-plane (M-set) fractal ; parameters: beta global: int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif init: complex aa = #pixel complex orbit[2] ; n_orbits orbit[0] = -1 orbit[1] = +1 if @giveup == "other" orbit[0] = orbit[1] endif bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live if @giveup=="all done" n_live = n_orbits else n_live = 1 endif loop: int ix = 0 repeat if continue[ix] z = orbit[ix] if (|z|>@bailout) continue[ix] = false n_live = n_live - 1 else z = z * aa * (z^2 - 3) + @bb; orbit[ix] = z endif ; abs endif ; continue ix = ix + 1 until ix >= n_orbits || n_live <= 0 bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Mandel(alpha) of polynomial (deg.3)" ; [v1.0] rating = notRecommended periodicity = 0 complex param bb caption="beta" default=(1,1) endparam heading caption = "Bailout" endheading float param bailout caption="Bailout" default=1e10 min=1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" default = 0 endparam switch: type = "Poly3_J" aa = #pixel bb = bb bailout = bailout } # - - - - - - - - - - - - - - - - - - - - - - - Poly3_Mb { ; iteration of f(z) := alpha*z^3 - 3*alpha*z + beta; ; beta-plane (M-set) fractal ; parameters: alpha global: int n_orbits if @giveup=="all done" || @giveup=="any done" n_orbits = 2 else n_orbits = 1 endif init: complex bb = #pixel complex orbit[2] ; n_orbits orbit[0] = 1 orbit[1] = -1 if @giveup == "other" orbit[0] = orbit[1] endif bool continue[2] ; n_orbits continue[0] = true continue[1] = true int n_live if @giveup=="all done" n_live = n_orbits else n_live = 1 endif loop: int ix = 0 repeat if continue[ix] z = orbit[ix] if (|z|>@bailout) continue[ix] = false n_live = n_live - 1 else z = z * @aa * (z^2 - 3) + bb; orbit[ix] = z endif ; abs endif ; continue ix = ix + 1 until ix >= n_orbits || n_live <= 0 bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue n_live > 0 default: title = "Mandel(beta) of polynomial (deg.3)" ; [v1.0] rating = notRecommended periodicity = 0 complex param aa caption="alpha" default=exp( (0,2) * #pi / 3 ) endparam heading caption = "Bailout" endheading float param bailout caption="Bailout" default=1e10 min=1 endparam param giveup caption = "Stop when" enum = "all done" "any done" "one" "other" default = 0 endparam switch: type = "Poly3_J" bb = #pixel aa = aa bailout = bailout } # --- Lacunary ------------------------------------------------------ Lacunary1 { ; iteration of f(z) = c + z + z^2 + z^4 + ... + z^(2^n) ; Escape-time (Julia-type) fractal init: z = #pixel loop: complex summand = z int k=0 z = z+@c while (k<@degree) summand = sqr(summand) z = z + summand k = k+1 endwhile bailout: |z| < @bailout default: title = "Lacunary 1" int param degree caption="degree n" hint="If your coloring algorithm requires a 'power' or 'exponent', set it to 2^n. \ If you use large values here, you may have to set a small bailout or increase the precision." default=2 min=1 endparam complex param c default=0.0 endparam float param bailout caption = "Bailout value" hint = "Iteration stops when z becomes larger than this bailout." min=0 default=1e4 endparam switch: type="Lacunary1_M" degree=degree bailout=bailout } # - - - - - - - - - - - - - - - - - - - - - - - Lacunary1_M { ; iteration of f(z) = c + z + z^2 + z^4 + ... + z^(2^n) init: if (@degree==1) z = -0.5 elseif (@degree==2) z = -0.38545849852963 elseif (@degree==3) z = -0.3828986970212 else z = -0.38289643077689 endif loop: complex summand = z int k=0 z = z+#pixel while (k<@degree) summand = sqr(summand) z = z + summand k = k+1 endwhile bailout: |z| < @bailout default: title = "Lacunary 1 Mandel" int param degree caption="degree n" hint="If your coloring algorithm requires a 'power' or 'exponent', set it to 2^n. \ If you use large values here, you may have to set a small bailout or increase the precision." default=2 min=1 endparam float param bailout caption = "Bailout value" hint = "Iteration stops when z becomes larger than this bailout." min=0 default=1e4 endparam switch: type="Lacunary1" c=#pixel degree=degree bailout=bailout } # - - - - - - - - - - - - - - - - - - - - - - - Lacunary3_J { ; Lacunary power series, escape-time (Julia-type) fractal global: ; calculate coefficients int powers[@power+1] float coeff[@power+1] powers[0] = 0 ; Initial power and coefficient if (@formula=="exp") coeff[0] = 1 else coeff[0] = 0 endif float qfac = 1.0 ; factorial for 'exp' int n=1 ; current power int maxp=0 ; max. index into powers[] and coeff[] while (n <= @power) maxp = maxp+1 powers[maxp] = n if (@formula=="exp") coeff[maxp] = qfac elseif (@formula=="log") coeff[maxp] = recip(n) else coeff[maxp] = 1 endif int n2 = ceil( n*@gap ) ; next power if (n2<=n) n2 = n+1 endif if (@formula=="exp") while (n=2) && (@power<=5) && @p_perturb int param power caption="Power" hint="If your coloring algorithm requires a 'power' or 'exponent', \ set it to this value. \ If you use very large values here, you may have to set a small bailout \ or increase the precision." default=22 min=2 endparam complex param c caption="c" hint="Switch to the M-type formula to pick a value." default=(0.25,0) endparam complex param p_a caption = "a" default = 1 ; compatibility with older version endparam bool param p_monic caption = "Monic Polynomial" default = false hint = "Use '1' instead of '1/Exponent' as highest coefficient. Default is false for compatibility with older versions." endparam bool param p_perturb caption = "Use perturbation when possible" default = true endparam heading caption = "Bailout" endheading float param bailout caption="Bailout value" hint="Iteration stops when z becomes larger than this bailout. \ Values below 2 will 'clip' the fractal. \ If you set 'Power' to a large value, you may have to reduce the bailout or increase the precision." min=0 default=1e10 endparam switch: type="CircleLogo_M" power=power bailout=bailout p_a = p_a p_monic = p_monic p_perturb = p_perturb } # - - - - - - - - - - - - - - - - - - - - - - - CircleLogo_M { ; C-plane fractal for f(z) = (z^power)/power + a*z + c ; Parameter-plane (Mandel-type) fractal \$define DEBUG global: float rPower = 1 if ! @p_monic rPower = recip(@power) endif int iPowerM1 = @power-1 complex startpoint complex startpts[@power-1] if @crit_choice == "Chosen point" startpoint = @user_start else startPoint = calcOneCriticalPoint() int shift = 0 if @crit_choice == "One critical point" shift = (shift + @crit_index) % iPowerM1 endif if (shift != 0) startpoint = startpoint * unarg( shift/iPowerM1 ) endif \$ifdef DEBUG print("startpoint=",startpoint,", a=",@p_a) complex spwr = startpoint^(@power-1) if @p_monic print("dfdz(startpoint)=", @power * spwr + @p_a) else print("dfdz(startpoint)=", spwr + @p_a) endif \$endif if @crit_choice == "All critical points" startpts[0] = startpoint int stix = 1 while (stix1) ;; && (real(@power)==exponent) int dexp = exponent-1 if dexp<0 dexp = -dexp endif if (dexp > 0) complex startPwr = - @p_a if @p_monic startPwr = startPwr / @power endif if @p_mangle && (imag(startPwr) == 0) ; real parameter float flotStartPwr = real(startPwr) float flotRoot = abs(flotStartPwr) ^ (1/dexp) int which = dexp%4 if flotStartPwr>0 return flotRoot elseif (which==1) || (which==3) return - flotRoot elseif which==2 return flip(flotRoot) endif endif return startPwr^(1/dexp) endif endif print("No critical points!") return 1/0 endfunc complex func unarg( float value ) return exp( (0,2) * #pi * value ) endfunc init: #z = startpoint bool continue = true float rad0 = -1 if @crit_choice == "All critical points" bool alive[@power-1] int ix=0 while (ix0 else if @p_monic #z = #pixel + #z * (@p_a + (z^(@power-1))) else #z = #pixel + #z * (@p_a + rPower*(z^(@power-1))) endif endif bailout: (@crit_choice == "All critical points" && continue) || (@crit_choice != "All critical points" && |#z| < @bailout) perturbinit: #dz = 0 perturbloop: if @power==2 if @p_monic #dz = 2*#dz*#z + sqr(#dz) + @p_a*#dz + #dpixel else #dz = #dz*#z + sqr(#dz)/2 + @p_a*#dz + #dpixel endif elseif @power==3 complex sqdz = sqr(#dz) if @p_monic #dz = 3*#dz*#z*(#z+#dz) + sqdz*#dz + @p_a*#dz + #dpixel else #dz = #dz*#z*(#z+#dz) + sqdz*#dz/3 + @p_a*#dz + #dpixel endif elseif @power==4 complex sqdz = sqr(#dz) if @p_monic #dz = #dz*#z*(4*sqr(#z) + 6*#z*#dz + 4*sqdz) + sqr(sqdz) + @p_a*#dz + #dpixel else #dz = #dz*#z*(sqr(#z) + 1.5*#z*#dz + sqdz) + sqr(sqdz)/4 + @p_a*#dz + #dpixel endif elseif @power==5 complex sqdz = sqr(#dz) ; #dz^2 complex p3dz = #dz * sqdz ; #dz^3 complex temp1 = p3dz + #z * (2*sqdz + #z * (2*#dz + #z)) if @p_monic #dz = #dz*#z*5*temp1 + sqdz*p3dz + @p_a*#dz + #dpixel else #dz = #dz*#z*temp1 + sqdz*p3dz/5 + @p_a*#dz + #dpixel endif endif default: title = "Circle-Logo Mandel" rating = recommended periodicity = 0 perturb = (@power>=2) && (@power<=5) && (@crit_choice != "All critical points") && @p_perturb heading caption="Formula parameters" endheading int param power caption="Power" hint="If your coloring algorithm requires a 'power' or 'exponent', set it to this value. \ If you use very large values here, you may have to set a small bailout or increase the precision." default=22 min=2 endparam complex param p_a caption = "a" default = 1 ; compatibility with older version endparam bool param p_monic caption = "Monic Polynomial" default = false hint = "Use '1' instead of '1/Exponent' as highest coefficient. Default is false for compatibility with older versions." endparam bool param p_perturb caption = "Use perturbation when possible" default = true endparam heading caption = "Bailout" endheading float param bailout caption="Bailout value" hint="Iteration stops when z becomes larger than this bailout. \ If you set 'Power' to a large value, you may have to reduce the bailout or increase the precision." min=0 default=1e10 endparam param crit_choice caption="Point to iterate" hint="Deceides which point(s) to iterate" enum="One critical point" "All critical points" "Chosen point" endparam int param crit_index caption="Index of critical point" hint="Use a number between 0 and 'Power'-1" default=0 min=0 visible = @crit_choice=="One critical point" endparam bool param @p_mangle caption="'Manage' critical points" default=true visible = @crit_choice!="Chosen point" endparam complex param user_start caption = "Start of iteration" hint="Custom iteration start" default = (1,0) visible = @crit_choice=="Chosen point" endparam switch: type="CircleLogo_J" c=#pixel power=power p_a=p_a p_monic = p_monic bailout=bailout p_perturb = p_perturb } ############################################################################# ### non-Analytic functions ### ############################################################################# frothy_basin { ; The "frothy basin" formula discovered by James C. Alexander ; and used as 'type=frothybasin' in Fractint. ; Adapted for UltraFractal 3.0 by Andreas Hoppler ; ; Iteration of f(z) = z^@power - @c*conj(z) ; ; 'Periodicity checking' will mess up this formula - better turn it off. ; 'Guessing' doesn't work well with frothier values of c. global: ; convergence and divergence radii float rad_inner = @delta float rad_outer = @bailout if (@bailout <= 0) if (@delta <= 0) rad_outer = 1e16 ; pick a reasonable default rad_inner = 1e-16 else rad_outer = 1/@delta endif elseif (@delta <= 0) rad_inner = 1/@bailout endif if rad_inner > rad_outer float rad_tmp = rad_inner rad_inner = rad_outer rad_outer = rad_tmp endif float radi2 = sqr(rad_inner) ; rotation by 0, 120, 240 degrees complex rot120[3] ; rot120[i]^3 = 1 for any i rot120[0] = 1 ; not used in formulas rot120[1] = -0.5 + 0.5*flip(sqrt(3)) rot120[2] = -0.5 - 0.5*flip(sqrt(3)) float half_a = 0.5 * imag(@c) float xcrit = 1.0 - sqrt( 1 + 0.75*sqr(imag(@c)) ) int triage_step = 0 if imag(@c) > 0 triage_step = 1 else triage_step = 2 endif bool attr3 = abs(imag(@c)) >= sqrt( 2/3 ) && abs(imag(@c)) <= 1.028713768218725 bool stripes = @test_inner == "tri-step" || @test_inner == "triage" || @test_inner == "tribble" init: z=#pixel complex z_prev = z - (0,1) complex z_step[3] int step = 0; for stripes, which stripe are we currently stepping through? int found = -1; for stripes, wich stripe should bail? int iter = 0; iteration counter, to see if we've got enough steps ; if we want to treat UF the current point as 'inside'. bool force_inside = false loop: bool continue = true; bailout criterion ; once force_inside is set, this formula just loops empty if (! force_inside) if (!stripes) || (step == 0) if (@test_for != "convergence") if |z| >= rad_outer if (@test_for == 3) ; divergence inside ; #solid = false ; can't do that here. Curses! force_inside = true continue = true else continue = false endif endif endif if (@test_for != "divergence") && continue if (@test_inner == "in-a-line") || (@test_inner == "4-in-a-line") if (iter >= 3) continue = sqr(real(z_step[0]) * imag(z_step[1]) - imag(z_step[0]) * real(z_step[1])) >= radi2 * (|z_step[0]|*|z_step[1]|) if (@test_inner == "4-in-a-line") && continue && iter >= 4 continue = sqr(real(z_step[1]) * imag(z_step[2]) - imag(z_step[1]) * real(z_step[2])) >= radi2 * (|z_step[1]|*|z_step[2]|) endif endif elseif (@test_inner == "scalar product") || (@test_inner == "scalar product 4") if (iter >= 3) continue = abs(real(z_step[0]) * imag(z_step[1]) - imag(z_step[0]) * real(z_step[1])) >= rad_inner if (@test_inner == "scalar product 4") && continue && iter >= 4 continue = abs(real(z_step[1]) * imag(z_step[2]) - imag(z_step[1]) * real(z_step[2])) >= rad_inner endif endif elseif (@test_inner == "step") continue = iter < 1 || |z_step[0]| > rad_inner elseif (@test_inner == "point trap") continue = |z-@zero| > rad_inner else ; 'imag' or one of the 'triangle' variants float xfound = -10 if (abs(imag(z)-half_a) <= rad_inner) found = 0 xfound = real(z) endif if (@test_inner != "imag") int ix_f = 1 while (ix_f < 3) && (found < 0) complex zrot = z * rot120[ix_f] if abs(imag(zrot)-half_a) <= rad_inner found = ix_f xfound = real(zrot) endif ix_f = ix_f + 1 endwhile endif if (@test_inner=="triage") if (found > 0) || (found == 0 && attr3) if (xfound > xcrit) found = (found + triage_step)%3 endif endif elseif (@test_inner=="tribble") if ((iter % 2) == 1) && (found >= 0) found = (5-found) % 3 endif endif endif endif ; test for convergence if (continue) && (found<0) ; ok, we need to calculate the next z here z_prev = z int rept = 0 while rept < @iteration z = z^@power - @c*conj(z) rept = rept+1 endwhile ; and update the trail z_step[2] = z_step[1] z_step[1] = z_step[0] z_step[0] = z - z_prev iter = iter+1 endif endif ; step 0 if (continue) if (stripes) continue = found != step step = step+1 if step >= 3 step = 0 endif else continue = found < 0 endif endif endif ; force_inside bailout: ; all the tests have been done in the loop section continue default: title = "frothy basin" rating = recommended periodicity = 0 ; turn off periodicity checking by default as it messes up this formula heading caption = "Fractal shape" endheading complex param c caption="c" default=(1, 0.7) hint="The c in the formula z^k - c*conj(z). \ For 'standard' images, use real(c) == 1.0 \ Certain bailout tests (imag, tri...) require real(c)=1.0" endparam int param power caption = "power k" default = 2 hint="The Power k in the formula z^k - c*conj(z)" endparam int param iteration caption = "function iteration" default = 1 min = 1 endparam heading caption = "Bailout" endheading param test_for caption = "Test for" enum = "convergence" "divergence" "conv. and div." "conv. and inside div." default = 2 endparam param test_inner caption = "convergence" enum = "in-a-line" "4-in-a-line" "scalar product" "scalar product 4" \ "imag" "triangle" "tri-step" "triage" "tribble" \ "step" "point trap" hint = "Convergence test to be used. \ 'in-a-line': stop if 3 iterations are in a straight line. \ '4-in-a-line': stop if 4 iterations are in a straight line. \ 'scalar product', 'scalar product 4': simpler versions of the straight line tests . \ 'imag': tests for imag(z) close to 0.5*imag(c). Works only if power==2 and real(c)==1. \ 'triangle': version of 'imag' with 3 lines in 120° angles. \ 'tri-step','triage','tribble': triangle variants for use with 3 stripes. \ 'step': stop if z stops moving. \ 'point trap': stop if z gets close to a fixed point." default = 5 enabled = @test_for != 1 endparam float param delta caption = "small bailout" hint="If <= 0, use reciprocal of large bailout" default=1e-16 enabled = @test_for != 1 endparam float param bailout caption = "large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=-1 enabled = @test_for != 0 endparam complex param zero caption = "Trap point" default = (0.0, 0.0) enabled = @test_inner=="point trap" && @test_for != 1 endparam ; switch: ; no switch for this formula -- yet. } # --- noise ---------------------------------------------------------------- aho_Perlin { global: ; init randomizer int permut[@permsize+@permsize+2] float gradient[@permsize+@permsize+2,2] int permsize2 = @permsize+@permsize+2 int seed = @randseed int iperm = 0 while (iperm < permsize2) permut[iperm] = iperm iperm = iperm+1 endwhile iperm = 0 while (iperm < permsize2) seed = random(seed) int newperm = abs(seed) % permsize2 int tmpperm = permut[iperm] permut[iperm] = permut[newperm] permut[newperm] = tmpperm iperm = iperm+1 endwhile iperm = 0 float gscale while iperm < permsize2 float gx float gy repeat seed = random(seed) gx = seed / #randomrange seed = random(seed) gy = seed / #randomrange gscale = sqr(gx) + sqr(gy) until gscale <= 1.0 && gscale > 1e-16 gscale = gscale^(-0.5) gradient[iperm,0] = gx*gscale gradient[iperm,1] = gy*gscale iperm = iperm+1 endwhile complex oct_mult = @oct_power * exp( (0,2) * #pi * @angle / 180.0 ) complex oct_init = oct_mult^@first_oct init: z = #pixel int iterate = 0 bool continue = false loop: complex zold = z complex z1 = z z = z * @z_amplitude int octave = @first_oct float amp = @amplitude complex frequency = oct_init while octave <= @last_oct ; loop through octaves complex z2 = z1 * frequency ; scale point to grid int oshift = sqr(octave) % @permsize int modulus if (@tile) modulus = round(cabs(frequency)) else modulus = @permsize endif if (modulus > @permsize) ; cutoff high freqs if buffer too small octave = @last_oct ; force bail out else ; ok ; calculate corners (bx_, by_) and distance from corners (rx_, ry_) int bx0 = floor(real(z2)) int by0 = floor(imag(z2)) float rx0 = real(z2) - bx0 float ry0 = imag(z2) - by0 float rx1 = rx0-1 float ry1 = ry0-1 bx0 = bx0 % modulus by0 = by0 % modulus if (bx0 < 0) bx0 = bx0 + modulus endif if (by0 < 0) by0 = by0 + modulus endif int bx1 = (bx0 + 1) if (bx1 >= modulus) bx1 = 0 endif int by1 = (by0 + 1) if (by1 >= modulus) by1 = 0 endif ; interpolation float sx float sy if @smooth sx = sqr(rx0) * (3 - 2*rx0) ; sx = s_curve(rx0); sy = sqr(ry0) * (3 - 2*ry0) ; sy = s_curve(ry0); else sx = rx0 sy = ry0 endif bool is_imag = true repeat ; for real and imag part int ip0 = permut[bx0+oshift] int ip1 = permut[bx1+oshift] if (is_imag) ip0 = permut[ip0] ip1 = permut[ip1] endif ip0 = ip0 % @permsize ip1 = ip1 % @permsize ; gradient indexes for the four corners int b00 = permut[ ip0 + by0 ] int b01 = permut[ ip0 + by1 ] int b10 = permut[ ip1 + by0 ] int b11 = permut[ ip1 + by1 ] float u float v u = gradient[b00,0]*rx0 + gradient[b00,1]*ry0 v = gradient[b10,0]*rx1 + gradient[b10,1]*ry0 float a = u + sx*(v-u) u = gradient[b01,0]*rx0 + gradient[b01,1]*ry1 v = gradient[b11,0]*rx1 + gradient[b11,1]*ry1 float b = u + sx*(v-u) float res = a + sy*(b-a) res = res*amp ; print( "octave=", octave, "res=", res) if (is_imag) z = z + flip(res) is_imag = false else z = z + res is_imag = true endif until is_imag ; print( "zold=", zold, "z=", z) endif ; prepare for next octave amp = amp * @persistence frequency = frequency * oct_mult octave = octave + 1 endwhile iterate = iterate+1 if (@bailtype==0) continue = false else continue = |z-zold| > @delta endif \$ifdef DEBUG if (iterate >= #maxiter - 20) print(iterate, ": ", z, " d= ", |z-zold| ) endif \$endif bailout: continue default: title = "Perlin noise" ; [v1.2] rating = recommended center = (0.5, 0.5) magn = 2 periodicity = 0 float param formula_version caption = "formula version" visible = false enabled = false default = 1.2 endparam heading caption = "fractal" endheading complex param z_amplitude caption = "z amplitude" default = 1.0 hint = "z[n+1] = z_amplitude * z[n] + noise_amplitude * noise(z[n])" endparam bool param tile endparam param bailtype caption = "Bailout on" enum = "1 iteration" "z - z_prev" default = 1 endparam float param delta caption = "small bailout" default = 1e-6 enabled = @bailtype == 1 endparam heading caption = "Noise" endheading int param first_oct caption = "first octave" default = 0 endparam int param last_oct caption = "last octave" default = 6 endparam float param oct_power caption = "octave power" default = 2 endparam float param angle default = 0 endparam float param persistence default = 0.5 endparam float param amplitude caption = "noise amplitude" default = 1/3 hint = "z[n+1] = z_amplitude * z[n] + noise_amplitude * noise(z[n])" endparam bool param smooth default = true endparam heading caption = "Random" endheading int param randseed default = 12345 endparam int param permsize caption = "buffer size" default = 8192 hint = "set this larger than the octave_power^last_octave" endparam } # --- Hénon attractor ----------------------------------------------- Henon_escape_J { global: float a = real(@c) float b = imag(@c) if (@mode == "rotated") a = real(@c) / sqr(imag(@c) ) b = 1 / imag(@c) ; print( a, "/", b ) elseif (@mode == "backward") b = 1 / imag(@c) endif init: z=#pixel loop: float x1 = real(z) float y1 = imag(z) float xnew float ynew if (@mode == "forward" || @mode == "rotated") xnew = y1 + 1 - a*sqr(x1) ynew = x1 * b else xnew = y1 * b ynew = x1 - 1 + a*sqr(xnew) endif z = xnew + flip(ynew) bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue |z| < @bailout default: title = "Hénon escape fractal (J)" magn = 0.5 float param formula_version caption = "formula version" visible = false enabled = false default = 1.0 endparam heading caption = "Fractal shape" endheading param mode enum = "forward" "backward" "rotated" default = 1 endparam complex param c default = (1.4, 0.3) endparam heading caption = "Bailout" endheading float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=1e16 endparam switch: type = "henon_escape_m" z_init = #pixel bailout = bailout mode = mode } # - - - - - - - - - - - - - - - - - - - - - - - Henon_escape_M { ; global: init: z = @z_init float a = real(#pixel) float b = imag(#pixel) if (@mode == "rotated") a = real(#pixel) / sqr(imag(#pixel) ) b = 1 / imag(#pixel) ; print( a, "/", b ) elseif (@mode == "backward") b = 1 / imag(#pixel) endif loop: float x1 = real(z) float y1 = imag(z) float xnew float ynew if (@mode == "forward" || @mode == "rotated") xnew = y1 + 1 - a*sqr(x1) ynew = x1 * b else xnew = y1 * b ynew = x1 - 1 + a*sqr(xnew) endif z = xnew + flip(ynew) bailout: ; -- this isn't the bailout but the continuing condition ; false = bail; true = continue |z| < @bailout default: title = "Hénon escape fractal (M)" rating = recommended magn = 0.5 float param formula_version caption = "formula version" visible = false enabled = false default = 1.0 endparam heading caption = "Fractal shape" endheading param mode enum = "forward" "backward" "rotated" default = 1 endparam complex param z_init caption = "initial z" default = (0, 0) endparam heading caption = "Bailout" endheading float param bailout caption="large bailout" hint="Radius for escape to infinity. If <= 0, use reciprocal of small bailout" default=1e16 endparam switch: type = "henon_escape_j" c = #pixel bailout = bailout mode = mode } ##################################################################### ### Plug-in (UltraFractal 5+) ### ##################################################################### aho_GenericSwitchFormula { \$ifdef VER50 global: import "common.ulb" import "aho.ulb" aho_BaseSwitchFormula switchF = wrapFormula() aho_Orbiter orb = new @p_colorOrbiter(0) orb.setFormula(switchF) if @p_isMandel switchF.SetMandelMode() else switchF.SetJuliaMode(@p_c) endif aho_BaseSwitchFormula func wrapFormula() Formula one = new @formulaClass(0) aho_BaseSwitchFormula resultFormula = aho_BaseSwitchFormula(one) if resultFormula==0 aho_SwitchFormulaWrapper rapper = new aho_SwitchFormulaWrapper(0) rapper.wrap(one) return rapper elseif @p_isMandel GenericArray instances = 0 if @p_whichCriticalPoint == "All until done" instances = getInstanceArray(resultFormula) if instances != 0 aho_AllBailFormulaWrapper wrap1 = new aho_AllBailFormulaWrapper(0) wrap1.wrap(instances) return wrap1 endif elseif @p_whichCriticalPoint == "One until done" instances = getInstanceArray(resultFormula) if instances != 0 aho_AnyBailFormulaWrapper wrap2 = new aho_AnyBailFormulaWrapper(0) wrap2.wrap(instances) return wrap2 endif endif endif return resultFormula endfunc GenericArray func getInstanceArray(aho_BaseSwitchFormula one) int numCrits = one.getNumberOfCriticalPoints() if numCrits <= 1 return 0 else if numCrits>100 numCrits = 100 ; QQQ higher? endif GenericArray resultArray = new GenericArray(numCrits) resultArray.m_elements[0] = one int ix = 1 while ix0) return critArray.m_elements[(@p_criticalPointIndex+len) % len] endif endif return fallback endfunc loop: #z = switchF.Iterate(#z) bool didBail = switchF.IsBailedOut(#z) if didBail orb.Finish(#z, switchF.getTriggeredBailout()) #z = orb.ComplexResult() else orb.Iterate(#z) endif bailout: ! didBail default: title = "Plug-in (switchable)" rating = recommended bool param p_isJulia default = false visible = false endparam bool param p_isMandel default = true visible = false endparam Formula param formulaClass caption = "Fractal Formula" default = aho_HRingSwitchFormula endparam heading caption = "Mandel shape" visible = @p_isMandel endheading param p_whichCriticalPoint caption = "Iterate" enum = "Default" "Selected Critical Point" "All until done" "One until done" "Custom point" default = 0 visible = @p_isMandel endparam int param p_criticalPointIndex caption = "Critical point index" hint = "0 = first (default) critical point" min = 0 default = 0 visible = @p_isMandel && @p_whichCriticalPoint == "Selected Critical Point" endparam complex param p_iterationStart caption = "Start point" visible = @p_isMandel && @p_whichCriticalPoint == "Custom point" default = (0.4,0.6) endparam int param p_mandelVariant default = 1 visible = false endparam heading caption = "Julia shape" visible = @p_isJulia endheading complex param p_c caption = "c" hint = "Generic parameter" default = (0, 0) visible = @p_isJulia endparam heading caption = "Coloring" endheading aho_Orbiter param p_colorOrbiter caption = "Color orbiter" default = aho_NoneOrbiter endparam switch: type = "aho_GenericSwitchFormula" p_colorOrbiter = p_colorOrbiter formulaClass = formulaClass p_isMandel = p_isJulia p_isJulia = p_isMandel p_mandelVariant = p_mandelVariant p_power = p_power p_bailout = p_bailout p_c = #pixel \$endif } ##################################################################### ### THE END ### #####################################################################