comment { Common Classes for Ultra Fractal 5

common.ulb 1.0
by Damien M. Jones, Frederik Slijkerman, Jim Blue, David Makin, Ron Barnett and others
May 19, 2017

This is a set of common classes that handle routine tasks within fractal formulas. They are designed to provide a common way to build these kinds of objects, to encourage authors to create inter-operable objects.

Although this file comes with Ultra Fractal 5, it is an active project that is updated from time to time. You may retrieve an updated library file (along with any other formulas published in the public database) by choosing Options > Update Public Formulas from inside Ultra Fractal, or by visiting the formula database and downloading the formula collection:

That site will also provide an online version of the library documentation, which may be easier to read than these comments.

Editing note: tabs are set to 4 for this code. } ; ----------------------------------------------------------------------------------------------------------------------------- ; Non-Generic classes class Math { ; Mathematical tools library. ; ;

; This class contains a number of useful functions for ; performing root-solving and other high-level math functions. ; They are implemented as static ; methods so that you do not have to maintain an instance ; of this class; you need only reference the class itself. ;

; For example, to solve 4x^2 - 4x + 1 = 0: ;

; myroots = new ComplexArray(2) ; Math.SolveRealQuadraticPolynomial(4,-4,1,0,myroots) ;

; myroots.m_Elements[0] will be 0.5 and myroots.m_Elements[1] ; will be -0.5. Note that you don't have to create a Math ; object first, you just use the function directly. ;

#### Root Solver General Notes

;

; If you are constructing polynomials inside a loop, you ; should create your roots array object just once, outside ; the loop, and re-use it each time. This is much faster ; that constructing a new roots object for each polynomial. ;

; Root solvers are provided for quadratic, cubic, and quartic ; functions, with either real or complex coefficients. Using the ; real version (if you know your coefficients are only real ; values) may be faster in some cases; note that real coefficients ; may still produce complex roots. ; ;

; In general, the pflags parameter controls some aspects of root-finding. ; It is the sum of these values: ;

; 1: Return real roots only (only available in the real-coefficient versions).
; 2: If a root occurs more than once, only return it once.
; 4: Use iterative improvement for roots (only available on cubic and quartic).
; 8: Emphasize the accuracy of smaller-magnitude roots, at the expense of accuracy in larger-magnitude roots (only available on cubic and quartic).
; 16: Force real roots by setting discriminant to zero if it's negative (quadratic only; used internall by the quartic solver). ;

; For example, pflags = 5 will return real roots only and use iterative improvement. ;

; Because of round-off errors, real roots (especially duplicate roots) may be ; found to have small imaginary parts. If real roots are specified, imaginary ; parts less than a threshold (in magnitude) are set to zero. This threshold is ; currently hard-coded inside these formulas and is generally acceptable, but ; may not be sufficient for your specific use. ;

; Similarly, if duplicate roots are to be removed and two roots differ by less ; than the threshold, they will be considered equal. ; ;

#### Specific Solver Notes

;

; Consider ;

; a*x^n + b*x^(n-1) + c*x^(n-2) + ... = 0 ;

##### Order n = 2, quadratic.
;

; This solution is well known. ;

; x = (-b +- sqrt(b*b-4*a*c))/(2*a) ;

; But there can be cancellation and loss of accuracy. Use this for the root that ; doesn't lose accuracy. Since x1*x2=c/a, use x2 = c/(a*x1) for the other root. ;

; For real coefficients, the two roots are either both real or a complex-conjugate ; pair. ; ;

##### Order n = 3, cubic.
;

; Divide by a, then use the substitution x = t-b/(3*a) to reduce to the "depressed" ; form, ;

; t^3 + aa*t + bb = 0 ;

; For real coefficients, the three roots are either all real or one real and a ; complex-conjugate pair. ; ;

##### Order n = 4, quartic.
;

; Divide by a, then use the substitution x = t-b/(4*a) to reduce to the "depressed" ; form, ;

; t^4 + aa*t^2 + bb*t + cc = 0 ; There are several ways (at least 5) to proceed. All of them require the solution ; of a cubic polynomial. ;

; For real coefficients, the four roots are in pairs. In each pair, there are ; either two reals or a complex-conjugate pair. ; ;

#### Accuracy

;

; If there are multiple roots, round-off may cause the roots to have a small ; imaginary part. Use the "real" flag setting if you want only real roots. ;

; Quadratic: no other problems. ;

; Cubic and Quartic: If there is a large variation in the magnitude of the roots, ; the smaller roots may be found less accurately. Some of this is inherent: the ; polynomial is ill-conditioned (small changes in the coefficients lead to larger ; changes in the roots). Sometimes it helps to solve for the roots of w = 1/x. ;

; Iterative refinement also sometimes helps. In solving p(x) = 0, suppose we have ; a "close enough" guess at a root, x0. Then Newton's method can be used to ; improve it, i.e., ;

; x1 = x0 - p(x0)/p'(x0) ;

; which converges quadratically if the root is a simple one. Quadratically means ; that the number of incorrect digits is halved at each iteration. ;

; However, Newton's method converges only linearly to a multiple root. (For a ; double root, the error is cut in half, approximately, at each iteration.) ;

; Instead we use a formula that includes the second derivative as well: ;

; x1 = x0 - p(x0)*p'(x0)/[ p'(x0)^2 - k*p(x0)*p"(x0)] with k = 1 ;

; This converges quadratically even with multiple roots. Note that setting k = 1/2 ; gives Halley's method, which converges cubically for simple roots, but only ; linearly with multiple roots. ;

; NOTE: If higher accuracy is needed, use the "Additional Precision" setting on ; the Formula tab in Layer Properties. ;

; Currently all thresholds are set to 1e-6, a "reasonable" value. ;

; ;

#### Bessel Function Notes

;

; These approximations to the Bessel functions J_0(x) and J_1(x) ; (for real x) are accurate to 15 or 16 digits. They come from
; Computer Approximations
; J. F. Hart, et al.
; Wiley, 1968
; Spot values have been checked from the text and as calculated by ; Mathematica. ;

; These "Bessel" functions J_0(z) and J_1(z) ; (for complex z) are not quite the actual Bessel functions. ; Let R = cabs(z). For R <= R1, these are the actual Bessel functions, ; computed by power series, to full precision except for round-off. ; For R > R1, these are ; similar to the actual Bessel functions, with appropriate asymptotic ; behavior for large R, and matched at R = R1. The actual Bessel functions ; cannot be calculated with full precision for large R without going to ; higher-precision calculations, and this is (a) not possible in UF and ; (b) not really necessary for making fractals. ;

; R1 = 8 is used below. ;

; Jim Blue, June 2008 public: ; Constructor ;

; You should never need to actually construct a Math object, as all of its ; methods are static and can be called without an object (just use the class ; name). func Math() endfunc ; Returns the sign of a number. ; ; @param px the number to test ; @return the sign; negative numbers give -1, positive numbers (including zero) return 1 static float func Sgn(const float px) if (px >= 0.0) return 1.0 else return -1.0 endif endfunc ; Determine the maximum of two numbers. ; ; @param px first number ; @param py second number ; @return the largest of the two numbers static float func MaxF(const float px, const float py) if (px >= py) return px else return py endif endfunc ; Determine the minimum of two numbers. ; ; @param px first number ; @param py second number ; @return the smallest of the two numbers static float func MinF(const float px, const float py) if (px <= py) return px else return py endif endfunc ; Solve (real) quadratic equation. ;

; See the class overview for details on the solvers. ; This solves a generic quadratic equation of the form ;

; (a1)(x^2) + (a2)(x) + (a3) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited) ; @return the number of roots found (will be 2 unless return values are limited) static int func SolveRealQuadraticPolynomial( \ const float pa1, const float pa2, const float pa3, \ const int pflags, ComplexArray &proots) bool forcerealroots = (floor(pflags/16)%2 == 1) ;bool smallrootsbest = (floor(pflags/8)%2 == 1) ;bool iterimprove = (floor(pflags/4)%2 == 1) bool removeduplicates = (floor(pflags/2)%2 == 1) bool realrootsonly = ( (pflags )%2 == 1) ; Degenerate case: not a quadratic if (pa1 == 0) if (pa2 == 0) return 0 else proots.m_Elements[0] = -pa3/pa2 return 1 endif endif ; Degenerate case: zero root if (pa3 == 0) proots.m_Elements[0] = 0 if (removeduplicates && pa2 == 0) return 1 endif proots.m_Elements[1] = -pa2/pa1 return 2 endif float RealThreshold = 1.e-6 float d = pa2*pa2 - 4*pa1*pa3 ; discriminant if (d < 0 && forcerealroots) d = 0 endif if (d < 0) int nroots = 2 float i2pa1 = 0.5 / pa1 float im = sqrt(-d) * i2pa1 if (realrootsonly || removeduplicates) if (abs(im) <= RealThreshold) im = 0 if (removeduplicates) nroots = 1 endif else return 0 endif endif float re = -pa2 * i2pa1 proots.m_Elements[0] = re + flip(im) proots.m_Elements[1] = re - flip(im) return nroots else ; d >= 0 so do carefully to avoid possible cancellation float t = -0.5*(pa2 + sgn(pa2)*sqrt(d)) proots.m_Elements[0] = t/pa1 ; this root is larger in magnitude proots.m_Elements[1] = pa3/t if (removeduplicates && \ abs(real(proots.m_Elements[0]-proots.m_Elements[1])) < RealThreshold) return 1 else return 2 endif endif endfunc ; end of SolveRealQuadraticPolynomial ; Solve (complex) quadratic equation. ;

; See the class overview for details on the solvers. ; This solves a generic quadratic equation of the form ;

; (a1)(x^2) + (a2)(x) + (a3) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited) ; @return the number of roots found (will be 2 unless return values are limited) static int func SolveComplexQuadraticPolynomial( \ const complex pa1, const complex pa2, const complex pa3, \ const int pflags, ComplexArray &proots) bool removeduplicates = (floor(pflags/2)%2 == 1) ; Degenerate case: not a quadratic if (pa1 == 0) if (pa2 == 0) return 0 else proots.m_Elements[0] = -pa3/pa2 return 1 endif endif ; Degenerate case: zero root if (pa3 == 0) proots.m_Elements[0] = 0 if (removeduplicates && pa2 == 0) return 1 endif proots.m_Elements[1] = -pa2/pa1 return 2 endif float Thresholdsq = sqr(1.e-6) complex sqrtd = sqrt(pa2*pa2 - 4*pa1*pa3) ; do carefully to avoid possible cancellation ; (this gives the same result as the real case when sqrtd is real) complex t1 = (pa2 + sqrtd) complex t2 = (pa2 - sqrtd) complex t if (|t1| >= |t2|) t = -0.5*t1 else t = -0.5*t2 endif proots.m_elements[0] = t/pa1 proots.m_elements[1] = pa3/t if (removeduplicates && \ |proots.m_elements[0] - proots.m_elements[1]| < Thresholdsq) return 1 else return 2 endif endfunc ; end of SolveComplexQuadraticPolynomial ; Solve (real) cubic equation. ;

; See the class overview for details on the solvers. ; This solves a generic cubic equation of the form ;

; (a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pa4 a4 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited) ; @return the number of roots found (will be 3 unless return values are limited) static int func SolveRealCubicPolynomial( \ const float pa1, const float pa2, const float pa3, const float pa4, \ const int pflags, ComplexArray &proots) if (pa1 == 0) ; degenerate-not a cubic return SolveRealQuadraticPolynomial(pa2, pa3, pa4, pflags, proots) endif if (pa4 == 0) ; one root of zero int n = SolveRealQuadraticPolynomial(pa1, pa2, pa3, pflags, proots) proots.m_Elements[n] = 0 return n+1 endif bool smallrootsbest = (floor(pflags/8)%2 == 1) bool iterimprove = (floor(pflags/4)%2 == 1) bool removeduplicates = (floor(pflags/2)%2 == 1) bool realrootsonly = ( (pflags )%2 == 1) ; scale out leading coefficient float a float b float c if (smallrootsbest) ; w = 1/x, solve for roots of w a = pa3/pa4 b = pa2/pa4 c = pa1/pa4 else a = pa2/pa1 b = pa3/pa1 c = pa4/pa1 endif ; Now have x^3 + aa*x + bb = 0 ; change of variable t = x - a/3 removes the quadratic term, giving ; t^3 - 3*Q*t + 2*R = 0 float a3 = a/3 float aa = b - 3*a3*a3 float bb = c - a3*(b - 2*a3*a3) float Q = -aa/3 ; (a*a-3*b)/9 float R = bb/2 ; (a*(2*a*a-9*b)+27*c)/54 float D = Q*Q*Q-R*R ; Could here do special cases Q=0, R=0, D=0 ; If two equal roots, D may be slightly negative, making them into a ; c.c. pair with small imaginary part. int nroots = 3 if (D >= 0) ; three real roots ; Note - theta usually done as acos, but atan2 is better float theta = atan2(R +flip(sqrt(D)))/3 ; float theta = acos(R/sqrt(Q*Q*Q))/3 ; is the usual version float cc = -2*sqrt(Q) ; D > 0, so Q > 0 proots.m_Elements[0] = cc*cos(theta ) -a3 proots.m_Elements[1] = cc*cos(theta + 2*#pi/3) -a3 proots.m_Elements[2] = cc*cos(theta - 2*#pi/3) -a3 else ; one real root, two complex-conjugate roots ; no duplicates possible float Ar = -Sgn(R) * ((abs(R) + sqrt(-D)))^(1/3) float Br = 0 if (Ar != 0) ; always Ar!=0 with real coefficients Br = Q/Ar endif proots.m_Elements[0] = Ar+Br - a3 float im = 0.5*sqrt(3)*(Ar-Br) proots.m_Elements[1] = -0.5*(Ar+Br) -a3 + flip(im) proots.m_Elements[2] = conj(proots.m_Elements[1]) endif if (smallrootsbest) int n = 0 while (n < nroots) proots.m_Elements[n] = 1/proots.m_Elements[n] n = n + 1 endwhile endif ; iterative refinement if (iterimprove) int n = 0 while (n < nroots) complex x = proots.m_Elements[n] complex v0 = pa4 + x*(pa3 + x*(pa2 + pa1*x)) complex v1 = pa3 + x*(2*pa2 + 3*pa1*x) if (v0*v1 != 0) complex v2 = 2*pa2 + 6*pa1*x ; complex xnew = x - v0*v1/(v1*v1-v0*v2) ; x = xnew x = x - v0*v1/(v1*v1-v0*v2) endif proots.m_Elements[n] = x n = n + 1 endwhile endif if (realrootsonly) nroots = RealOnly(nroots, proots) endif if (removeduplicates) nroots = SingleOnly(nroots, proots) endif return nroots endfunc ; end of SolveRealCubicPolynomial ; Solve (complex) cubic equation. ;

; See the class overview for details on the solvers. ; This solves a generic cubic equation of the form ;

; (a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pa4 a4 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited) ; @return the number of roots found (will be 3 unless return values are limited) static int func SolveComplexCubicPolynomial( \ const complex pa1, const complex pa2, const complex pa3, const complex pa4, \ const int pflags, ComplexArray &proots) if (pa1 == 0) ; degenerate-not a cubic return SolveComplexQuadraticPolynomial(pa2, pa3, pa4, pflags, proots) endif if (pa4 == 0) ; one root of zero int n = SolveComplexQuadraticPolynomial(pa1, pa2, pa3, pflags, proots) proots.m_Elements[n] = 0 return n+1 endif bool smallrootsbest = (floor(pflags/8)%2 == 1) bool iterimprove = (floor(pflags/4)%2 == 1) bool removeduplicates = (floor(pflags/2)%2 == 1) ;bool realrootsonly = ( (pflags )%2 == 1) ; scale out leading coefficient complex a complex b complex c if (smallrootsbest) ; w = 1/x, solve for roots of w a = pa3/pa4 b = pa2/pa4 c = pa1/pa4 else a = pa2/pa1 b = pa3/pa1 c = pa4/pa1 endif ; Now have x^3 + a*x^2 + b*x + c = 0 ; change of variable t = x - a/3 removes the quadratic term, giving ; t^3 - 3*Q*t + 2*R = 0 complex a3 = a/3 complex aa = b - 3*a3*a3 complex bb = c - a3*(b - 2*a3*a3) complex Q = -aa/3 ; (a*a-3*b)/9 complex R = bb/2 ; (a*(2*a*a-9*b)+27*c)/54 complex D = Q*Q*Q-R*R complex sm = sqrt(-D) complex t1 = r + sm complex t2 = r - sm complex Ar if (|t1| >= |t2|) Ar = -(t1^(1/3)) else Ar = -(t2^(1/3)) endif complex Br = 0 if (Ar != 0) Br = Q/Ar endif proots.m_Elements[0] = Ar + Br - a3 proots.m_Elements[1] = -0.5*(Ar+Br) -a3 + (0,0.5)*sqrt(3)*(Ar-Br) proots.m_Elements[2] = -0.5*(Ar+Br) -a3 - (0,0.5)*sqrt(3)*(Ar-Br) int nroots = 3 if (smallrootsbest) int n = 0 while (n < nroots) proots.m_Elements[n] = 1/proots.m_Elements[n] n = n + 1 endwhile endif ; iterative refinement if (iterimprove) int n = 0 while (n < nroots) complex x = proots.m_Elements[n] complex v0 = pa4 + x*(pa3 + x*(pa2 + pa1*x)) complex v1 = pa3 + x*(2*pa2 + 3*pa1*x) if (v0*v1 != 0) complex v2 = 2*pa2 + 6*pa1*x ; complex xnew = x - v0*v1/(v1*v1-v0*v2) ; x = xnew x = x - v0*v1/(v1*v1-v0*v2) endif proots.m_Elements[n] = x n = n+1 endwhile endif if (removeduplicates) nroots = SingleOnly(nroots, proots) endif return nroots endfunc ; end of SolveComplexCubicPolynomial ; Solve (real) quartic equation. ;

; See the class overview for details on the solvers. ; This solves a generic quartic equation of the form ;

; (a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pa4 a4 term ; @param pa5 a5 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited) ; @return the number of roots found (will be 4 unless return values are limited) static int func SolveRealQuarticPolynomial( \ const float pa1, const float pa2, const float pa3, const float pa4, \ const float pa5, const int pflags, ComplexArray &proots) if (pa1 == 0) ; degenerate: not a quartic return SolveRealCubicPolynomial(pa2, pa3, pa4, pa5, pflags, proots) endif if (pa5 == 0) ; one root of zero int n = SolveRealCubicPolynomial(pa1, pa2, pa3, pa4, pflags, proots) proots.m_Elements[n] = 0 return n+1 endif bool smallrootsbest = (floor(pflags/8)%2 == 1) bool iterimprove = (floor(pflags/4)%2 == 1) bool removeduplicates = (floor(pflags/2)%2 == 1) bool realrootsonly = ( (pflags )%2 == 1) float a float b float c float d ; scale out leading coefficient if (smallrootsbest) ; w = 1/x, solve for roots of w a = pa4/pa5 b = pa3/pa5 c = pa2/pa5 d = pa1/pa5 else a = pa2/pa1 b = pa3/pa1 c = pa4/pa1 d = pa5/pa1 endif ; Now have x^4 + a*x^3 + b*x^2 + c*x + d = 0 ; change of variable t = x - a/4 removes the cubic term, giving ; the "depressed quartic" ; t^4 + aa*t^2 + bb*t + cc = 0 float a4 = a/4 float aa = b - 6*a4*a4 float bb = c - 2*a4*(b - 4*a4*a4) float cc = d - a4*(c-a4*(b-3*a4*a4)) ; Could here do any special cases ; The next substitution is from ; http://en.wikipedia.org/wiki/Quartic_equation ; It gives the same result as Faucette, who is quoted on Wolfram's site, ; http://mathworld.wolfram.com/QuarticEquation.html ; ; Factor the depressed quartic into two quadratic terms ; (t^2 +p*t + q)*(t^2 + r*t + s) = 0 ; Equating terms gives ; 0 = p + r so r = -p ; aa = q - p^2 + s ; bb = p * (s - q) ; cc = q * s ; Some algebra gives a cubic in p^2 ; p^6 + 2*aa*p^4 + (aa^2-4*cc)*p^2 - b^2 = 0 ; Solve this and choose a real root for p^2. ; if p^2 > 0, one sign of the square root gives ; s and the other sign gives q, ; then 2s = aa +-b/p +p^2 ; then 2q = aa -+b/p +p^2 ; but we choose the better equation for s and then q = cc/s. ; if p^2 = 0, use the aa and cc equations to get a quadratic for s and q. float c1 = 2*aa float c2 = aa*aa-4*cc float c3 = -bb*bb int n1 = SolveRealCubicPolynomial(1, c1, c2, c3, 1, proots) ; Since -p is the sum of 2 roots of the quartic, at least one of the 3 roots ; of the cubic is real and non-negative (within round-off). ; If all roots of the quartic are real, all 3 roots of the cubic are real ; and non-negative. Also, n1 is 1 or 3. float p2 if (n1 == 1) p2 = real(proots.m_Elements[0]) else p2 = Maxf(real(proots.m_Elements[0]), real(proots.m_Elements[1])) p2 = Maxf(p2, real(proots.m_Elements[2])) ; minimum positive root maybe gives most accuracy ; float p2_0, float p2_1, float p2_2 ; p2_0 = real(proots.m_Elements[0]) ; p2_1 = real(proots.m_Elements[1]) ; p2_2 = real(proots.m_Elements[2]) ; Sort3f(p2_0, p2_1, p2_2) ; p2 = p2_0 ; if (p2 + RealThreshold < 0) ; p2 = p2_1 ; if (p2 + RealThreshold < 0) ; p2 = p2_2 ; endif ; endif endif if (p2 < 0) ; only because of roundoff p2 = 0 endif float s float q int n3 = SolveRealQuadraticPolynomial(1, -(aa+p2), cc, 16, proots) if (n3 == 0) ; can't happen -- means p2 = 0, but round-off... return 0 endif s = real(proots.m_Elements[0]) q = real(proots.m_Elements[1]) float p = sgn(bb)*sgn(s-q)*sqrt(p2) ; Now solve the two quadratics; since the leading term is 1, each has two ; roots, either both real or c.c. int n2 = SolveRealQuadraticPolynomial(1, p, q, 0, proots) proots.m_Elements[2] = proots.m_Elements[0] - a4 proots.m_Elements[3] = proots.m_Elements[1] - a4 int n0 = SolveRealQuadraticPolynomial(1, -p, s, 0, proots) proots.m_Elements[0] = proots.m_Elements[0] - a4 proots.m_Elements[1] = proots.m_Elements[1] - a4 int nroots = n0 + n2 if (nroots == 0) ; shouldn't happen return 0 endif if (smallrootsbest) int n = 0 while (n < nroots) proots.m_Elements[n] = 1/proots.m_Elements[n] n = n + 1 endwhile endif ; iterative refinement if (iterimprove) int n = 0 while (n < nroots) complex x = proots.m_Elements[n] complex v0 = pa5 + x*(pa4 + x*(pa3 + x*(pa2 + pa1*x))) complex v1 = pa4 + x*(2*pa3 + x*(3*pa2 + 4*pa1*x)) if (v0*v1 != 0) ; x = x - v0/v1 complex v2 = 2*pa3 + x*(6*pa2 + 12*pa1*x) x = x - v0*v1/(v1*v1-v0*v2) endif proots.m_Elements[n] = x n = n + 1 endwhile endif if (realrootsonly) nroots = RealOnly(nroots, proots) endif if (removeduplicates) nroots = SingleOnly(nroots, proots) endif return nroots endfunc ; end of SolveRealQuartic Polynomial ; Solve (complex) quartic equation. ;

; See the class overview for details on the solvers. ; This solves a generic quartic equation of the form ;

; (a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0 ; ; @param pa1 a1 term ; @param pa2 a2 term ; @param pa3 a3 term ; @param pa4 a4 term ; @param pa5 a5 term ; @param pflags options ; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited) ; @return the number of roots found (will be 4 unless return values are limited) static int func SolveComplexQuarticPolynomial( \ const complex pa1, const complex pa2, const complex pa3, const complex pa4, \ const complex pa5, const int pflags, ComplexArray &proots) if (pa1 == 0) ; degenerate: not a quartic return SolveComplexCubicPolynomial(pa2, pa3, pa4, pa5, pflags, proots) endif if (pa5 == 0) ; one root of zero int n = SolveComplexCubicPolynomial(pa1, pa2, pa3, pa4, pflags, proots) proots.m_Elements[n] = 0 return n+1 endif bool smallrootsbest = (floor(pflags/8)%2 == 1) bool iterimprove = (floor(pflags/4)%2 == 1) bool removeduplicates = (floor(pflags/2)%2 == 1) ;bool realrootsonly = ( (pflags )%2 == 1) complex a complex b complex c complex d ; scale out leading coefficient if (smallrootsbest) ; w = 1/x, solve for roots of w a = pa4/pa5 b = pa3/pa5 c = pa2/pa5 d = pa1/pa5 else a = pa2/pa1 b = pa3/pa1 c = pa4/pa1 d = pa5/pa1 endif ; Now have x^4 + a*x^3 + b*x^2 + c*x + d = 0 ; change of variable t = x - a/4 removes the cubic term, giving ; the "depressed quartic" ; t^4 + aa*t^2 + bb*t + cc = 0 complex a4 = a/4 complex aa = b - 6*a4*a4 complex bb = c - 2*a4*(b - 4*a4*a4) complex cc = d - a4*(c-a4*(b-3*a4*a4)) ; Could here do any special cases ; The next substitution is from ; http://en.wikipedia.org/wiki/Quartic_equation ; It gives the same result as Faucette, who is quoted on Wolfram's site, ; http://mathworld.wolfram.com/QuarticEquation.html ; ; Factor the depressed quartic into two quadratic terms ; (t^2 +p*t + q)*(t^2 + r*t + s) = 0 ; Equating terms gives ; 0 = p + r so r = -p ; aa = q - p^2 + s ; bb = p * (s - q) ; cc = q * s ; Some algebra gives a cubic in p^2 ; p^6 + 2*aa*p^4 + (aa^2-4*cc)*p^2 - b^2 = 0 ; Solve this and choose a real root for p^2. ; if p^2 > 0, one sign of the square root gives ; s and the other sign gives q, ; then 2s = aa +-b/p +p^2 ; then 2q = aa -+b/p +p^2 ; but we choose the better equation for s and then q = cc/s. ; if p^2 = 0, use the aa and cc equations to get a quadratic for s and q. complex c1 = 2*aa complex c2 = aa*aa-4*cc complex c3 = -bb*bb SolveComplexCubicPolynomial(1, c1, c2, c3, 0, proots) ; Since -p is the sum of 2 roots of the quartic, at least one of the 3 roots ; of the cubic is real and non-negative (within round-off). ; If all roots of the quartic are real, all 3 roots of the cubic are real ; and non-negative. Also, n1 is 1 or 3. complex p2 = proots.m_Elements[0] if (|p2| < |proots.m_Elements[1]|) p2 = proots.m_Elements[1] endif if (|p2| < |proots.m_Elements[2]|) p2 = proots.m_Elements[2] endif complex s complex q int n3 = SolveComplexQuadraticPolynomial(1, -(aa+p2), cc, 0, proots) if (n3 == 0) ; can't happen -- means p2 = 0, but round-off... return 0 endif s = proots.m_Elements[0] q = proots.m_Elements[1] complex p = sqrt(p2) complex t1 = bb - p * (s-q) complex t2 = bb = p * (s-q) if (|t1| > |t2|) p = -p endif ; Now solve the two quadratics int n2 = SolveComplexQuadraticPolynomial(1, p, q, 0, proots) proots.m_Elements[2] = proots.m_Elements[0] - a4 proots.m_Elements[3] = proots.m_Elements[1] - a4 int n0 = SolveComplexQuadraticPolynomial(1, -p, s, 0, proots) proots.m_Elements[0] = proots.m_Elements[0] - a4 proots.m_Elements[1] = proots.m_Elements[1] - a4 int nroots = n0 + n2 if (smallrootsbest) int n = 0 while (n < nroots) proots.m_Elements[n] = 1/proots.m_Elements[n] n = n + 1 endwhile endif ; iterative refinement if (iterimprove) int n = 0 while (n < nroots) complex x = proots.m_Elements[n] complex v0 = pa5 + x*(pa4 + x*(pa3 + x*(pa2 + pa1*x))) complex v1 = pa4 + x*(2*pa3 + x*(3*pa2 + 4*pa1*x)) if (v0*v1 != 0) ; x = x - v0/v1 complex v2 = 2*pa3 + x*(6*pa2 + 12*pa1*x) x = x - v0*v1/(v1*v1-v0*v2) endif proots.m_Elements[n] = x n = n + 1 endwhile endif if (removeduplicates) nroots = SingleOnly(nroots, proots) endif return nroots endfunc ; end of SolveRealQuartic Polynomial ; Complex Gamma() David Makin June 2008 ; @param pz complex argument ; @return Gamma function of argument static complex func ComplexGamma(complex pz) complex w = pz ; float gplus = 5.2421875 if pz==0.0 w = recip(0) else if real(pz)<0.0 w = -pz if imag(pz)==0.0 if (real(pz)%1)==0.0 w = recip(0) endif endif endif if real(w)<1e300 w = (sqrt(2.0*#pi)* \ (0.99999999999999709182 \ + 57.156235665862923517/(w+1.0) \ - 59.597960355475491248/(w+2.0) \ + 14.136097974741747174/(w+3.0) \ - 0.49191381609762019978/(w+4.0) \ + .33994649984811888699E-4/(w+5.0) \ + .46523628927048575665E-4/(w+6.0) \ - .98374475304879564677e-4/(w+7.0) \ + .15808870322491248884E-3/(w+8.0) \ - .21026444172410488319e-3/(w+9.0) \ + .21743961811521264320e-3/(w+10.0) \ - .16431810653676389022e-3/(w+11.0) \ + .84418223983852743293e-4/(w+12.0) \ - .26190838401581408670e-4/(w+13.0) \ + .36899182659531622704e-5/(w+14.0)) \ /w) * exp(-w-5.2421875)*(w+5.2421875)^(w+0.5) if real(pz)<0.0 w = pz*w*sin(#pi*pz) if (w==0.0) w = recip(0) else w = -#pi/w endif endif endif endif return w endfunc ; Real Gamma() David Makin June 2008 ; @param x float argument ; @return Gamma function of argument static float func RealGamma(float x) float w = x ; float gplus = 5.2421875 if x==0.0 w = recip(0) else if x<0.0 w = -x if (x%1)==0.0 w = recip(0) endif endif if w<1e300 w = (sqrt(2.0*#pi)* \ (0.99999999999999709182 \ + 57.156235665862923517/(w+1.0) \ - 59.597960355475491248/(w+2.0) \ + 14.136097974741747174/(w+3.0) \ - 0.49191381609762019978/(w+4.0) \ + .33994649984811888699E-4/(w+5.0) \ + .46523628927048575665E-4/(w+6.0) \ - .98374475304879564677e-4/(w+7.0) \ + .15808870322491248884E-3/(w+8.0) \ - .21026444172410488319e-3/(w+9.0) \ + .21743961811521264320e-3/(w+10.0) \ - .16431810653676389022e-3/(w+11.0) \ + .84418223983852743293e-4/(w+12.0) \ - .26190838401581408670e-4/(w+13.0) \ + .36899182659531622704e-5/(w+14.0)) \ /w) * exp(-w-5.2421875)*(w+5.2421875)^(w+0.5) if x<0.0 w = x*w*sin(#pi*x) if (w==0.0) w = recip(0) else w = -#pi/w endif endif endif endif return w endfunc ; Complex Beta() David Makin June 2008 ; @param a complex argument ; @param b complex argument ; @return Beta function of arguments static complex func ComplexBeta(complex a, complex b) return ComplexGamma(a)*ComplexGamma(b)/ComplexGamma(a+b) endfunc ; Real Beta() David Makin June 2008 ; @param a float argument ; @param b float argument ; @return Beta function of arguments static float func RealBeta(float a, float b) return RealGamma(a)*RealGamma(b)/RealGamma(a+b) endfunc ; Real Bessel J0. ; @param x float argument ; @return Bessel function of order 0 static float func RealBesselJ0(const float x) ; JZERO 5846 float ax = abs(x) if (ax <= 8.0) float xx = x*x return ( \ +.181083899214164308956057747e9 + xx*( \ -.4446475956150241214875200307e8 + xx*( \ +.2629297609723391264477929426e7 + xx*( \ -.6635120296335057670760789901e5 + xx*( \ +.9000010866385970012308555062e3 + xx*( \ -.7411779995119348858162812772e1 + xx*( \ +.39771954892980603699828143e-1 + xx*( \ -.1446351863454593575005971793e-3 + xx*( \ +.3612785949562383677336235392e-6 + xx*( \ -.6089358750808789415483968656e-9 + xx*( \ +.6424736784979100876734782112e-12 + xx*( \ -.3313360272836170910319585771e-15)))))))))))) / ( \ +.1810838992141643042977083347e9 + xx*( \ +.806215242038697916042541338e6 + xx*( \ +.1415495011707417263126447751e4 + xx))) else float z = sqr(8.0/ax) ; PZERO 6546 float P0 = ( \ +.8554822541506661710252074e4 + z*( \ +.8894437532960619440762804e4 + z*( \ +.2204501043965180428995069e4 + z*( \ +.128677585748714193298851e3 + z*( \ +.9004793474802880316384e0))))) / ( \ +.8554822541506662842462151e4 + z*( \ +.8903836141709595355210343e4 + z*( \ +.2214048851914710419825683e4 + z*( \ +.130884900499923882835109e3 + z )))) ; QZERO 6946 float Q0 = (8.0/x) * ( \ -.37510534954957111594836e2 + z*( \ -.46093826814625174976377e2 + z*( \ -.13990976865960680088016e2 + z*( \ -.1049732798234554833126e1 + z*( \ -.93525953294031893049e-2)))))/ ( \ +.2400674237117267479318819e4 + z*( \ +.2971983745208491990065486e4 + z*( \ +.921566975526530895082307e3 + z*( \ +.74428389741411178824152e2 + z )))) z = ax - #pi/4 return sqrt(2/#pi/ax)*(P0*cos(z)-Q0*sin(z)) endif endfunc ; Real Bessel J1. ; @param x float argument ; @return Bessel function of order 1 static float func RealBesselJ1(const float x) ; JONE 6045 float ax = abs(x) if (ax <= 8.0) float xx = x*x return x * ( \ +.695364226329838502166085207e8 + xx*( \ -.8356785487348914291918495672e7 + xx*( \ +.3209027468853947029888682298e6 + xx*( \ -.58787877666568200462726094e4 + xx*( \ +.6121876997356943874446879769e2 + xx*( \ -.3983107983952332023421699105e0 + xx*( \ +.1705769264349617107854016566e-2 + xx*( \ -.4910599276555129440130592573e-5 + xx*( \ +.9382193365140744507653268479e-8 + xx*( \ -.1107352224453730633782671362e-10 + xx*( \ +.63194310317443161294700346e-14 ))))))))))) / ( \ +.139072845265967685120764336e9 + xx*( \ +.6705346835482299302199750802e6 + xx*( \ +.1284593453966301898121332163e4 + xx))) else float z = sqr(8.0/ax) ; PONE 6747 float P1 = ( \ +.1290918471896188077350689e5 + z*( \ +.1309042051103506486292571e5 + z*( \ +.313275295635506951011069e4 + z*( \ +.17431379748379024599685e3 + z*( \ +.122850537643590432633e1 )))))/ ( \ +.1290918471896187879332737e5 + z*( \ +.1306678308784401036110575e5 + z*( \ +.310928141677002883350924e4 + z*( \ +.16904721775008609992033e3 + z )))) ; QONE 7147 float Q1 = (8.0/ax) * ( \ +.14465282874995208675225e3 + z*( \ +.1744291689092425885102e3 + z*( \ +.5173653281836591636536e2 + z*( \ +.379944537969806734901e1 + z*( \ +.36363466476034710809e-1 )))))/ ( \ +.308592701333231723110639e4 + z*( \ +.373434010601630179517765e4 + z*( \ +.11191098527047487025919e4 + z*( \ +.8522392064341340397334e2 + z )))) z = ax - #pi*3/4 z = sqrt(2/#pi/ax)*(P1*cos(z)-Q1*sin(z)) if (x < 0) return -z else return z endif endif endfunc ; Complex Bessel J0. ; @param z complex argument ; @return "Bessel" function of order 0 static complex func ComplexBesselJ0(const complex z) float R = cabs(z) if (R <= 8) return ComplexBesselJ0_Inner(z) else complex zz = z * 8/R return ComplexBesselJ0_Outer(z) * (ComplexBesselJ0_Inner(zz)/ComplexBesselJ0_Outer(zz)) endif endfunc ; Complex Bessel J1. ; @param z complex argument ; @return "Bessel" function of order 0 static complex func ComplexBesselJ1(const complex z) float R = cabs(z) if (R <= 8) return ComplexBesselJ1_Inner(z) else complex zz = z * 8/R return ComplexBesselJ1_Outer(z) * (ComplexBesselJ1_Inner(zz)/ComplexBesselJ1_Outer(zz)) endif endfunc private: ; Helper function to eliminate duplicate roots. static int func SingleOnly(int nroots, ComplexArray &proots) float DuplThreshold = 1.e-12 int n = nroots while (n > 0) n = n - 1 complex xx = proots.m_Elements[n] bool dupl = false int k = 0 while (k < n && !dupl) ; check earlier roots against nth if (|xx -proots.m_Elements[k]| < DuplThreshold) ; nth is duplicate nroots = nroots - 1 dupl = true if (n < nroots) ; need to move down int m = n while (m < nroots) proots.m_Elements[m] = proots.m_Elements[m+1] m = m + 1 endwhile endif endif k = k+1 endwhile endwhile return nroots endfunc ; Helper function to eliminate complex roots. static int func RealOnly(int nroots, ComplexArray &proots) float RealThreshold = 1.e-6 int n = nroots while (n > 0) n = n - 1 if (abs(imag(proots.m_Elements[n])) <= RealThreshold) ; real proots.m_Elements[n] = real(proots.m_Elements[n]) else ; not real nroots = nroots - 1 if (n < nroots) ; need to move down int m = n while (m < nroots) proots.m_Elements[m] = proots.m_Elements[m+1] m = m + 1 endwhile endif endif endwhile return nroots endfunc ; Power series for Bessel Function J0. ;

; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 360. ;

; Number of terms to use determined experimentally. ; ; @param z complex argument ; @return Bessel function of order 0 via power series for small argument static complex func ComplexBesselJ0_Inner(const complex z) int Nmax = 3 + round(6*sqrt(cabs(z))) complex zz = sqr(z/2) complex term = 1 complex sum = term int n = 1 while (n <= Nmax) term = -term * zz / sqr(n) sum = sum + term n = n + 1 endwhile return sum endfunc ; Function with appropriate asymptotic behavior for Bessel Function J0. ;

; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 364. ;

; Two terms of P0 and Q0 are used, rather arbitrarily. ; ; @param z complex argument ; @return "Bessel" function of order 0 via asymptotic series for large argument static complex func ComplexBesselJ0_Outer(const complex z) complex z1 = 1/(8*z) complex z2 = sqr(z1) complex arg = z - #pi/4 return sqrt(2/#pi/z)*((1-z2*9/2)*cos(arg)+z1*(1-z2*9*25/6)*sin(arg)) endfunc ; Power series for Bessel Function J1. ;

; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 360. ;

; Number of terms to use determined experimentally. ; ; @param z complex argument ; @return Bessel function of order 1 via power series for small argument static complex func ComplexBesselJ1_Inner(const complex z) int Nmax = 3 + round(6*sqrt(cabs(z))) complex zz = sqr(z/2) complex term = z/2 complex sum = term int n = 1 while (n <= Nmax) term = -term * zz / (n*(n+1)) sum = sum + term n = n + 1 endwhile return sum endfunc ; Function with appropriate asymptotic behavior for Bessel Function J1. ;

; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 364. ;

; Two terms of P1 and Q1 are used, rather arbitrarily. ; ; @param z complex argument ; @return "Bessel" function of order 1 via asymptotic series for large argument static complex func ComplexBesselJ1_Outer(const complex z) complex z1 = 1/(8*z) complex z2 = sqr(z1) complex arg = z - #pi*3/4 return sqrt(2/#pi/z)*((1-z2*3*5/2)*cos(arg)+3*z1*(1-z2*5*21/6)*sin(arg)) endfunc } class Vector { ; Vector class. ; ;

; The Vector class allows you to easily define a four-dimensional vector. ; These have a variety of uses in fractal generation, including use as a ; container for a 3D vector. ;

; For performance reasons, none of the Vector methods return a vector. ; The last vector argument is the return value. For example, ;

; Vector V1 = new Vector(....)
; V1.Cross(V2,V3) ;

; is equivalent to ;

; V3 = V1 x V2 ;

; Also for performance reasons, you may access the vector elements ; directly, as m_x, m_y, m_z, and m_w. ;

; Ron Barnett public: ; Constructor ; ; @param x x term of vector ; @param y y term of vector ; @param z z term of vector ; @param w w term of vector (use 1.0 if you are storing a 3D vector) func Vector(float x, float y, float z, float w) Init(x,y,z,w) endfunc ; Re-initialize a Vector. ;

; If you need to completely change the value of a Vector, you can use this function rather than deleting the old Vector and creating a new one. (Re-use is faster.) ; ; @param x x term of vector ; @param y y term of vector ; @param z z term of vector ; @param w w term of vector (use 1.0 if you are storing a 3D vector) func Init(float x, float y, float z, float w) m_x = x m_y = y m_z = z m_w = w endfunc ; Copy a Vector to another Vector. ;

; If you simply assign the value of one Vector to another, you actually create a reference in both variables to the same Vector object, and changes in one variable's vector will appear in the other. See the UF help file topic on "reference" ("Objects") for details. To make a true copy from one Vector to another Vector, use this function. ; ; @param T target vector that will receive a copy of the Vector's data func Copy(Vector T) T.m_x = m_x T.m_y = m_y T.m_z = m_z T.m_w = m_w endfunc ; Get a Vector element by index. ;

; Sometimes you need to access a Vector's elements by index number (0-3) rather than by name. If you do, you may use this function. Note that this is slower than accessing the member variables directly. ; ; @param n element number (0-3) to retrieve ; @return the value of the element float func Get(int n) if (n == 0) return m_x elseif (n == 1) return m_y elseif (n == 2) return m_z elseif (n == 3) return m_w else return 0 endif endfunc ; Set a Vector element by index. ;

; Sometimes you need to set a Vector's elements by index number (0-3) rather than by name. If you do, you may use this function. Note that this is slower than accessing the member variables directly. ; ; @param n element number (0-3) to set ; @param v the new value for the element func Set(int n, float v) if (n == 0) m_x = v elseif (n == 1) m_y = v elseif (n == 2) m_z = v elseif (n == 3) m_w = v endif endfunc ; Normalize a Vector. ;

; Normalization gives a unit vector. The x, y, z and w values are ; direction cosines. ; ; @param V target vector that will receive the normalized vector func Normalize(Vector V) float vd = 1/sqrt(m_x^2+m_y^2+m_z^2+m_w^2) V.m_x = m_x*vd V.m_y = m_y*vd V.m_z = m_z*vd V.m_w = m_w*vd endfunc ; Determine a Vector's size. ; ; @return the size (or length) of the vector float func Size() return sqrt(m_x^2+m_y^2+m_z^2+m_w^2) endfunc ; Add two Vectors. ; ; @param V1 the second vector to add ; @param V2 the target vector for the added vectors to be stored in func Add(Vector V1, Vector V2) V2.m_x = m_x+V1.m_x V2.m_y = m_y+V1.m_y V2.m_z = m_z+V1.m_z V2.m_w = m_w+V1.m_w endfunc ; Subtract two Vectors. ; ; @param V1 the second vector to subtract ; @param V2 the target vector for the differenced vectors to be stored in func Sub(Vector V1, Vector V2) V2.m_x = m_x-V1.m_x V2.m_y = m_y-V1.m_y V2.m_z = m_z-V1.m_z V2.m_w = m_w-V1.m_w endfunc ; Multiply a vector by a constant. ; ; @param c the constant to multiply by ; @param V the target vector for the multiplied values to be stored in func MConst(float c, Vector V) V.m_x = m_x*c V.m_y = m_y*c V.m_z = m_z*c V.m_w = m_w*c endfunc ; Shift Vector elements. ;

; This function shifts the vector elements so that x becomes y, y becomes z, z becomes w, and w becomes x. You may shift by any number of positions, although any multiple of four will return the elements to their starting positions. Negative values will be interpreted as shifts in reverse. ; ; @param shift number of places to shift (negative values shift in reverse) ; @param V the target vector to store the shifted vector in func Shift(int shift, Vector V) float x = m_x float y = m_y float z = m_z float w = m_w shift = shift % 4 if shift < 0 shift = shift + 4 endif if shift == 1 V.m_x = w V.m_y = x V.m_z = y V.m_w = z elseif shift == 2 V.m_x = z V.m_y = w V.m_z = x V.m_w = y elseif shift == 3 V.m_x = y V.m_y = z V.m_z = w V.m_w = x else V.m_x = x V.m_y = y V.m_z = z V.m_w = w endif endfunc ; Vector Dot Product ;

; This gives the cosine of the angle between the vectors. ; ; @param V the second vector in the dot product ; @return the dot product float func Dot(Vector V) return V.m_x*m_x+V.m_y*m_y+V.m_z*m_z+V.m_w*m_w endfunc ; Vector Cross Product (3D) ;

; This gives the vector cross product, which given two 3D vectors will produce a third vector which is at right angles to the plane containing the two vectors. Note that changing the order of the two vectors will generally result in the cross product pointing in the opposite direction: ;

; V x V1 = - V1 x V ;

; This version of the cross product assumes the Vector objects hold a 3D vector, and returns a 3D vector. ; ; @param V1 the second vector in the cross product ; @param V2 the target vector to store the cross product vector in func Cross(Vector V1, Vector V2) float x = m_y*V1.m_z-m_z*V1.m_y float y = m_z*V1.m_x-m_x*V1.m_z float z = m_x*V1.m_y-m_y*V1.m_x V2.m_x = x V2.m_y = y V2.m_z = z V2.m_w = 1 endfunc ; Vector Cross Product (4D) ;

; This gives the vector cross product, which given three 4D vectors will produce a fourth vector which is at right angles to the hyperplane containing the three vectors. Note that changing the order of the three vectors will generally result in a different orientation of the result. ; ; @param V1 the second vector in the cross product ; @param V2 the third vector in the cross product ; @param V3 the target vector to store the cross product vector in func Cross4D(Vector V1, Vector V2, Vector V3) float a = (V1.m_x * V2.m_y) - (V1.m_y * V2.m_x); float b = (V1.m_x * V2.m_z) - (V1.m_z * V2.m_x); float c = (V1.m_x * V2.m_w) - (V1.m_w * V2.m_x); float d = (V1.m_y * V2.m_z) - (V1.m_z * V2.m_y); float e = (V1.m_y * V2.m_w) - (V1.m_w * V2.m_y); float f = (V1.m_z * V2.m_w) - (V1.m_w * V2.m_z); float x = (m_y * f) - (m_z * e) + (m_w * d); float y = - (m_x * f) + (m_z * c) - (m_w * b); float z = (m_x * e) - (m_y * c) + (m_w * a); float w = - (m_x * d) + (m_y * b) - (m_z * a); V3.m_x = x V3.m_y = y V3.m_z = z V3.m_w = w endfunc ; Rotate Vector in the YZ plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotYZ(float theta, Vector V) theta = theta*#pi/180 float y = m_y*cos(theta) - m_z*sin(theta) float z = m_y*sin(theta) + m_z*cos(theta) V.m_y = y V.m_z = z V.m_x = m_x V.m_w = m_w endfunc ; Rotate Vector in the ZX plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotZX(float theta, Vector V) theta = theta*#pi/180 float z = m_z*cos(theta) - m_x*sin(theta) float x = m_z*sin(theta) + m_x*cos(theta) V.m_z = z V.m_x = x V.m_y = m_y V.m_w = m_w endfunc ; Rotate Vector in the XY plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotXY(float theta, Vector V) theta = theta*#pi/180 float x = m_x*cos(theta) - m_y*sin(theta) float y = m_x*sin(theta) + m_y*cos(theta) V.m_x = x V.m_y = y V.m_z = m_z V.m_w = m_w endfunc ; Rotate Vector in the XW plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotXW(float theta, Vector V) theta = theta*#pi/180 float x = m_x*cos(theta) - m_w*sin(theta) float w = m_x*sin(theta) + m_w*cos(theta) V.m_x = x V.m_w = w V.m_z = m_z V.m_y = m_y endfunc ; Rotate Vector in the YW plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotYW(float theta, Vector V) theta = theta*#pi/180 float w = m_w*cos(theta) - m_y*sin(theta) float y = m_w*sin(theta) + m_y*cos(theta) V.m_y = y V.m_w = w V.m_x = m_x V.m_z = m_z endfunc ; Rotate Vector in the ZW plane. ; ; @param theta angle to rotate, in radians ; @param V target vector to store resulting rotated vector in func RotZW(float theta, Vector V) theta = theta*#pi/180 float w = m_w*cos(theta) - m_z*sin(theta) float z = m_w*sin(theta) + m_z*cos(theta) V.m_w = w V.m_z = z V.m_x = m_x V.m_y = m_y endfunc float m_x float m_y float m_z float m_w default: } class Array { ; Generic dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class. Normally ; you cannot pass arrays (in their entirety) as parameters ; to functions; if you instead create an Array object of ; the appropriate type, you can store all the values in ; the array and still pass it to a function. ;

; Also, with dynamic arrays, you cannot assign the contents ; of one array to another in a single assignment. If you ; do this with an Array object, you will be making both ; variables refer to the same array. Instead, use the ; Copy() function. ;

; Because accessing elements in an Array object require a ; reference to the object, there is an ever-so-slight ; performance cost for using an Array instead of a built-in ; dynamic array. This cost should be negligible. ;

; To use an Array class, create an object of the required ; type: ;

; myobject = new ComplexArray(3) ;

; You set the initial size when you create the object. You ; can change the size using SetArrayLength() just as you ; would use setLength() on a built-in dynamic array type: ;

; myobject.SetArrayLength(12) ;

; To access the array elements, use the m_Elements member: ;

; complex p = myobject.m_Elements[5] public: ; Constructor ; ; @param plength initial length of the array func Array(int plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return 0 endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) endfunc ; Copy (must be implemented in each type) ; func Copy(Array &dest) default: } class BooleanArray(Array) { ; Generic Boolean dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Boolean data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func BooleanArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest BooleanArray object to copy all values into; previous contents of dest will be discarded. func Copy(BooleanArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new BooleanArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) bool m_Elements[] default: } class IntegerArray(Array) { ; Generic Integer dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Integer data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func IntegerArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest IntegerArray object to copy all values into; previous contents of dest will be discarded. func Copy(IntegerArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new IntegerArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) int m_Elements[] default: } class FloatArray(Array) { ; Generic Float dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Float data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func FloatArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest FloatArray object to copy all values into; previous contents of dest will be discarded. func Copy(FloatArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new FloatArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) float m_Elements[] default: } class ComplexArray(Array) { ; Generic Complex dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Complex data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func ComplexArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest ComplexArray object to copy all values into; previous contents of dest will be discarded. func Copy(ComplexArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new ComplexArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) complex m_Elements[] default: } class ColorArray(Array) { ; Generic Color dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Color data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func ColorArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest ColorArray object to copy all values into; previous contents of dest will be discarded. func Copy(ColorArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new ColorArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) color m_Elements[] default: } class VectorArray(Array) { ; Generic Vector dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Vector data. See the Array base class for more ; information. public: ; Constructor ; ; @param plength initial length of the array func VectorArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest VectorArray object to copy all values into; previous contents of dest will be discarded. func Copy(VectorArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new VectorArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) if (m_Elements[j] == 0) dest.m_Elements[j] = 0 else if (dest.m_Elements[j] == 0) dest.m_Elements[j] = new Vector(m_Elements[j].m_x,m_Elements[j].m_y,m_Elements[j].m_z,m_Elements[j].m_w) else m_Elements[j].Copy(dest.m_Elements[j]) endif endif j = j + 1 endwhile endfunc ; Array elements (for direct access) Vector m_Elements[] default: } class ArrayArray(Array) { ; Generic Array dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Array data. Note that with this, you can create ; arrays of arrays, with each nested array having ; a possibly different number of elements. You can ; also create arrays of arrays of arrays. ;

; See the Array base class for more information. public: ; Constructor ; ; @param plength initial length of the array func ArrayArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy ;

; This performs a "deep" copy, where each element Array ; is fully copied to the target ArrayArray object. For ; large arrays, this will not be fast. ; ; @param dest ArrayArray object to copy all values into; previous contents of dest will be discarded. func Copy(ArrayArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new ArrayArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) if (BooleanArray(m_Elements[j])) BooleanArray ba = BooleanArray(dest.m_Elements[j]) if (!ba) dest.m_Elements[j] = 0 endif BooleanArray(m_Elements[j]).Copy(ba) dest.m_Elements[j] = ba elseif (IntegerArray(m_Elements[j])) IntegerArray ia = IntegerArray(dest.m_Elements[j]) if (!ia) dest.m_Elements[j] = 0 endif IntegerArray(m_Elements[j]).Copy(ia) dest.m_Elements[j] = ia elseif (FloatArray(m_Elements[j])) FloatArray fa = FloatArray(dest.m_Elements[j]) if (!fa) dest.m_Elements[j] = 0 endif FloatArray(m_Elements[j]).Copy(fa) dest.m_Elements[j] = fa elseif (ComplexArray(m_Elements[j])) ComplexArray ca = ComplexArray(dest.m_Elements[j]) if (!ca) dest.m_Elements[j] = 0 endif ComplexArray(m_Elements[j]).Copy(ca) dest.m_Elements[j] = ca elseif (ColorArray(m_Elements[j])) ColorArray pa = ColorArray(dest.m_Elements[j]) if (!pa) dest.m_Elements[j] = 0 endif ColorArray(m_Elements[j]).Copy(pa) dest.m_Elements[j] = pa elseif (VectorArray(m_Elements[j])) VectorArray va = VectorArray(dest.m_Elements[j]) if (!va) dest.m_Elements[j] = 0 endif VectorArray(m_Elements[j]).Copy(va) dest.m_Elements[j] = va elseif (ArrayArray(m_Elements[j])) ArrayArray aa = ArrayArray(dest.m_Elements[j]) if (!aa) dest.m_Elements[j] = 0 endif ArrayArray(m_Elements[j]).Copy(aa) dest.m_Elements[j] = aa elseif (GenericArray(m_Elements[j])) GenericArray ga = GenericArray(dest.m_Elements[j]) if (!ga) dest.m_Elements[j] = 0 endif GenericArray(m_Elements[j]).Copy(ga) dest.m_Elements[j] = ga else dest.m_Elements[j] = m_Elements[j] endif j = j + 1 endwhile endfunc ; Array elements (for direct access) Array m_Elements[] default: } class GenericArray(Array) { ; Generic object dynamic array object wrapper. ;

; This is a generic dynamic array wrapper class for ; Generic data. That is, you can create an array ; to hold any kind of object that is derived from ; the Generic base class. Unlike the other Array ; classes, elements in a GenericArray do not have ; to be of the same type. This gives you a lof of ; flexibility. However, when you access your objects ; from the m_Elements array, they will appear to be ; Generic objects rather than whatever class they ; really are. It is up to you to cast them to the ; correct type in your code; if you try to make a ; cast that is not valid (because the object is not ; of that type or derived from that type) UF will ; give a null (0) result for the cast. ;

; See the Array base class for more information. public: ; Constructor ; ; @param plength initial length of the array func GenericArray(int plength) setLength(m_Elements, plength) endfunc ; Length wrapper ; ; @return the number of elements in the array int func GetArrayLength() return length(m_Elements) endfunc ; setLength wrapper ;

; This is analagous to UF's native setLength() function ; for dynamic arrays. ; ; @param plength new length of the array func SetArrayLength(int plength) setLength(m_Elements, plength) endfunc ; Copy (must be implemented in each type) ; ; @param dest GenericArray object to copy all values into; previous contents of dest will be discarded. func Copy(GenericArray &dest) int l = this.GetArrayLength() if (dest == 0) dest = new GenericArray(l) else dest.SetArrayLength(l) endif int j = 0 while (j < l) dest.m_Elements[j] = m_Elements[j] j = j + 1 endwhile endfunc ; Array elements (for direct access) Generic m_Elements[] default: } ; ----------------------------------------------------------------------------------------------------------------------------- ; Base classes class Generic { ; Generic base class. ;

; This is a generic object class. Most of the other classes ; derive from this class. Its primary purpose is to act as ; glue between all of the objects that are used within a ; formula. ;

; Each class in UF derives from the built-in Object class, ; so this class may seem a bit redundant. But eventually ; this class will be extended to include some parameter-passing ; functions, which will require all major classes to derive ; from a common extensible parent. If your class is intended ; to be user-selectable and has any user-selectable parameters ; at all, you should probably derive it from this class. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Generic(Generic pparent) m_Parent = pparent endfunc ; Fetch parent object ; ; @return a reference to the object that created this object; ; you can repeatedly apply this function to the resulting object ; reference until a null (0) is returned, and the object that ; gave the null return value is the "top" object. Generic func GetParent() return m_Parent endfunc protected: ; A reference to the parent object. Generic m_Parent default: int param v_generic caption = "Version (Generic)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_generic < 100 endparam } class Transform(Generic) { ; Transform base class. ;

; This is a generic transform class. Its purpose is to take ; an arbitrary point and transform it into some other point. ; This roughly corresponds to UF's native transformation ; formula type. The advantage to writing a transform as a ; class is that then it can be used in many other places. ; Transform classes also support the notion that a sequence ; may yield a "solid" point, just as UF's transformation ; formulas do. ;

; If you are migrating an existing transformation formula ; to a Transform-derived class, the process is fairly ; straightforward. Look at your global: and transform: ; sections. Any variables that are shared between the two ; should be declared in your protected: class section. ; The code in your global: section goes into the constructor. ; The code in your transform: section goes into the Iterate() ; function. Use m_Solid in place of #solid and use pz in ; place of #pixel. Whenever possible, try to avoid using ; other built-in values like #width or #screenpixel as this ; will make your class behave oddly when not being used to ; transform an entire fractal. ;

; As an extra optional step, you may look at what you've ; placed in your Iterate() function and decide to refactor it. ; With a transformation formula, the transform: section is ; used only once per pixel. However, a Transform object can ; be used for this AND used to transform an entire sequence ; of values (for example, the entire sequence of iterated ; values from a fractal formula). If there are some computations ; you perform for your transformation that will be exactly the ; same for each of the values in that sequence, you may move ; them to the Init() function; this will improve the performance ; of your formula. Any variables holding the results of these ; computations will need to be declared in your protected: ; section so that they will be available to the Iterate() ; function. ;

; IMPORTANT NOTE: If your transform is intended to be a ; user-selected class, you should probably derive from the ; UserTransform class rather than directly from Transform. ; The Transform family also includes UtilityTransform and ; ClipShape, so objects looking to allow Transforms will ; normally require them to be in the UserTransform subclass. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Transform(Generic pparent) Generic.Generic(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). It will be called even if only one value is being ; transformed (e.g. a normal transformation formula). Use ; this to perform any setup that is exactly the same for ; each value in the sequence being transformed. ; ; @param pz the value representing the sequence; for a normal transformation formula use, this will be #pixel. In some cases this may differ from the first value passed to Iterate() if the calling code applies some other transformations. func Init(complex pz) ; As part of the base class functionality, we set up ; a simple iteration counter and a flag to indicate ; whether the sequence yields a solid color or not. ; If you would like to re-use this, you will need to ; call Transform.Init() in your derived class's Init() ; function if you provide an Init() function. (If you ; don't have an Init() function, UF will call the base ; class Transform.Init() automatically.) m_Iterations = 0 m_Solid = false endfunc ; Transform a single point within a sequence ;

; After a sequence has been set up with Init(), this function ; will be called once for each value in the sequence. Note ; that all values in the sequence must be processed in order ; (they cannot be processed out of order). If the sequence ; contains only one value, Init() will still be called and ; then Iterate() will be called just once. ; ; @param pz the complex value to be transformed ; @return the transformed value complex func Iterate(complex pz) ; Update the iteration count. As with Init(), if you ; provide your own Iterate() function, as you almost ; certainly will, you would need to call Transform.Iterate() ; to get this functionality. However, since you would not ; need the default transformation, it is better to call ; Transform.IterateSilent() instead to update this value. ; Even higher performance could be obtained by just ; updating the iteration count manually; however, if you ; do this and the base Transform class is updated, your ; code might not inherit the functionality because you ; are not calling the base class. m_Iterations = m_Iterations + 1 return pz endfunc ; Update internal counters without transforming a point ;

; For some Transform classes, the actual transformation being ; performed changes for each value in the sequence (e.g. ; TrapTransform, which can rotate each iteration by a different ; amount). In some cases the calling code may determine in ; advance that a particular point does not need to be transformed ; (perhaps because it is not being used) but it still needs to be ; accounted for. This function is used in place of Iterate() for ; those situations; the Transform should update any internal ; state that changes between iterations within a sequence. Since ; no value is being transformed, no parameters are passed and no ; return value is provided. func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc ; Test whether a sequence is solid-colored or not. ;

; This function takes the place of the #solid flag in a native ; UF transformation formula. Since this is an object which may be ; used in a context other than a transformation formula, #solid ; may not be available to us. Instead, Iterate() should set an ; internal flag m_Solid and the calling code should use IsSolid() ; to test the value of this flag. For some Transform classes this ; test may be very involved, which is why it is available as a ; function rather than simply giving the calling code access to ; the m_Solid variable directly. ; ; @return a Boolean flag indicating whether the sequence is solid-colored or not. bool func IsSolid() return m_Solid endfunc protected: int m_Iterations ; count the number of iterations bool m_Solid ; flag indicating whether last pixel was solid or not default: int param v_transform caption = "Version (Transform)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_transform < 100 endparam } class UserTransform(Transform) { ; User-selectable Transform base class. ;

; This class adds no new functionality to the Transform class. ; Its purpose is to provide a grouping for user-selectable ; Transform objects that actually transform complex points. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func UserTransform(Generic pparent) Transform.Transform(pparent) endfunc default: int param v_usertransform caption = "Version (UserTransform)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_usertransform < 100 endparam } class Formula(Generic) { ; Formula base class. ;

; This is a generic formula class. Its purpose is to take ; an arbitrary initial point and produce a sequence of ; points of indeterminate length. This roughly corresponds ; to UF's native fractal formula type. The advantage to ; writing a fractal formula as a class is that it can then ; be used in many other places. Formula classes also ; support the notion of "bailed out", as a way of indicating ; when the sequence of values is finished. ;

; If you are migrating an existing fractal formula to a ; Formula-derived class, the process is fairly straightforward. ; Any variables you set in your global: section that are used ; elsewhere in your formula should be declared in your protected: ; section. Move the code from your global: section into the ; constructor. Move the code from your init: section into the ; Init() function. Move the code from your loop: section into the ; Iterate() function. Move the code from your bailout: test into ; the IsBailedOut() function (but be sure it returns the logically ; opposite value; see the explanation on IsBailedOut()). ;

; If you have complicated bailout testing, you may instead put ; that code into your Iterate() function and save the result in ; the m_BailedOut variable and omit implementing IsBailedOut(). ; If you do this, remember to set m_BailedOut to false in your ; Init() function, as the base class implementation sets it to ; true to end iterating right away. ;

; Note that if you are creating a formula whose outside points ; escape towards infinity, you should derive from DivergentFormula ; as it will automatically include a bailout parameter and make ; that accessible to other objects. Similarly, if your formula ; is convergent, you should derive from ConvergentFormula. If ; your formula has both kinds of outside points, derive from ; ConvergentDivergentFormula. Otherwise, you will have to write ; your bailout test yourself. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Formula(Generic pparent) Generic.Generic(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence; for a normal fractal formula, this will be #pixel ; @return first value in the sequence; this corresponds to #z in a fractal formula complex func Init(complex pz) ; Base class implementation flags the sequence to end ; immediately. As long as you provide your own implementation ; of IsBailedOut(), this is irrelevant. We also clear the ; iteration counter. m_Iterations = 0 m_BailedOut = true return pz endfunc ; Produce the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. ; ; @param pz previous value in the sequence; corresponds to #z in a fractal formula. Note that you should always use this value for computing the next iteration, rather than a saved value, as the calling code may modify the returned value before passing it back to the next Iterate() call. ; @return the next value in the sequence complex func Iterate(complex pz) ; Base class implentation automatically updates the ; iteration counter and returns an unmodified pz ; (that is, the default sequence is the same value ; repeated infinitely). m_Iterations = m_Iterations + 1 return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; This is typically called once per iteration to test whether ; the sequence has completed. If producing this result is ; difficult to do separately from the Iterate() processing, ; you should produce the result in Iterate() and set the ; m_BailedOut flag appropriately, and leave IsBailedOut() ; unimplemented in your class to inherit this behavior. ;

; Note that this test is the OPPOSITE sense of the bailout: ; section in UF's fractal formulas. A bailout: section ; returns TRUE to continue iterating. IsBailedOut() must ; return FALSE to continue iterating. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Formula classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(complex pz) return m_BailedOut endfunc ; Determine the primary exponent. ;

; Many fractals can be characterized by an exponent value that ; is useful to other formulas, so we provide that here. If ; your formula does not need or use this value, override the ; p_power parameter and make it hidden. ; ; @return the primary exponent parameter complex func GetPrimaryExponent() return @p_power endfunc ; Determine the upper bailout boundary. ;

; By default, we return a simple fixed value. Divergent ; formulas will override this function and return the ; real parameter. ; ; @return the upper bailout parameter float func GetUpperBailout() return 256 endfunc ; Determine the lower bailout boundary. ;

; By default, we return a simple fixed value. Convergent ; formulas will override this function and return the ; real parameter. ; ; @return the lower bailout parameter float func GetLowerBailout() return 0 endfunc protected: int m_Iterations ; count the number of iterations bool m_BailedOut ; flag indicating whether sequence has bailed out or not default: int param v_formula caption = "Version (Formula)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_formula < 100 endparam complex param p_power caption = "Exponent" default = (2,0) hint = "Defines the primary exponent for the fractal." endparam } class DivergentFormula(Formula) { ; Divergent Formula base class. ;

; This class extends Formula to provide automatic, simple ; bailout tests. Formulas for which the outside points ; become arbitrarily large should derive from this class ; so as to avoid having to write the bailout tests and ; the GetUpperBailout() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func DivergentFormula(Generic pparent) Formula.Formula(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence; for a normal fractal formula, this will be #pixel ; @return first value in the sequence; this corresponds to #z in a fractal formula complex func Init(complex pz) Formula.Init(pz) ; Base class implementation flags the sequence to end ; immediately. We clear the flag here. m_BailedOut = false return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; Since this is a divergent fractal, the test is easy: if it's ; bigger than the bailout, the sequence is done. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Formula classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(complex pz) return |pz| > @p_bailout endfunc ; Determine the upper bailout boundary. ; ; @return the upper bailout parameter float func GetUpperBailout() return @p_bailout endfunc default: int param v_divergentformula caption = "Version (DivergentFormula)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_divergentformula < 100 endparam float param p_bailout caption = "Bailout" default = 1.0e20 hint = "Defines how soon an orbit bails out, i.e. doesn't belong to the fractal set anymore. Note: larger values result in more iterations being calculated." endparam } class ConvergentFormula(Formula) { ; Convergent Formula base class. ;

; This class extends Formula to provide automatic, simple ; bailout tests. Formulas for which the outside points ; converge on points should derive from this class ; so as to avoid having to write the bailout tests and ; the GetLowerBailout() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func ConvergentFormula(Generic pparent) Formula.Formula(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence; for a normal fractal formula, this will be #pixel ; @return first value in the sequence; this corresponds to #z in a fractal formula complex func Init(complex pz) Formula.Init(pz) ; Base class implementation flags the sequence to end ; immediately. We clear the flag here. m_BailedOut = false m_ZOld = pz return pz endfunc ; Produce the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. ; ; @param pz previous value in the sequence; corresponds to #z in a fractal formula. Note that you should always use this value for computing the next iteration, rather than a saved value, as the calling code may modify the returned value before passing it back to the next Iterate() call. ; @return the next value in the sequence complex func Iterate(complex pz) Formula.Iterate(pz) ; We need to update our previously-saved Z value. m_ZOld = pz return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; Since this is a divergent fractal, the test is easy: if it's ; bigger than the bailout, the sequence is done. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Formula classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(complex pz) return |pz-m_ZOld| < @p_bailout endfunc ; Determine the lower bailout boundary. ; ; @return the lower bailout parameter float func GetLowerBailout() return @p_bailout endfunc protected: complex m_ZOld default: int param v_convergentformula caption = "Version (ConvergentFormula)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_convergentformula < 100 endparam float param p_bailout caption = "Bailout" default = 1.0e-15 hint = "Defines how soon an orbit bails out, i.e. doesn't belong to the fractal set anymore. Note: smaller values result in more iterations being calculated." endparam } class ConvergentDivergentFormula(Formula) { ; Convergent/Divergent Formula base class. ;

; This class extends Formula to provide automatic, simple ; bailout tests. Formulas for which the outside points ; converge on points OR grow without bound should derive ; from this class so as to avoid having to write the ; bailout tests and the GetUpperBailout() and ; GetLowerBailout() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func ConvergentDivergentFormula(Generic pparent) Formula.Formula(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence; for a normal fractal formula, this will be #pixel ; @return first value in the sequence; this corresponds to #z in a fractal formula complex func Init(complex pz) Formula.Init(pz) ; Base class implementation flags the sequence to end ; immediately. We clear the flag here. m_BailedOut = false m_ZOld = pz return pz endfunc ; Produce the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. ; ; @param pz previous value in the sequence; corresponds to #z in a fractal formula. Note that you should always use this value for computing the next iteration, rather than a saved value, as the calling code may modify the returned value before passing it back to the next Iterate() call. ; @return the next value in the sequence complex func Iterate(complex pz) Formula.Iterate(pz) ; We need to update our previously-saved Z value. m_ZOld = pz return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; Since this is a divergent fractal, the test is easy: if it's ; bigger than the bailout, the sequence is done. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Formula classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(complex pz) return |pz| > @p_upperbailout || |pz-m_ZOld| < @p_lowerbailout endfunc ; Determine the upper bailout boundary. ; ; @return the upper bailout parameter float func GetUpperBailout() return @p_upperbailout endfunc ; Determine the lower bailout boundary. ; ; @return the lower bailout parameter float func GetLowerBailout() return @p_lowerbailout endfunc protected: complex m_ZOld default: int param v_convergentdivergentformula caption = "Version (ConvergentDivergentFormula)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_convergentdivergentformula < 100 endparam float param p_upperbailout caption = "Bailout" default = 1.0e20 hint = "Defines how soon an orbit bails out, i.e. doesn't belong to the fractal set anymore. Note: larger values result in more iterations being calculated." endparam float param p_lowerbailout caption = "Bailout" default = 1.0e-15 hint = "Defines how soon an orbit bails out, i.e. doesn't belong to the fractal set anymore. Note: smaller values result in more iterations being calculated." endparam } class Coloring(Generic) { ; Coloring base class. ;

; This is a generic coloring class. Its purpose is to take ; a sequence of points and reduce it to either a gradient ; index or an exact color (precisely which is left up to ; the derived class). ;

; For the most part, you will not use this formula directly. ; If you are migrating an existing coloring formula to a ; Coloring-derived class, you would normally choose to derive ; from either GradientColoring or DirectColoring. These sub- ; classes of Coloring automatically indicate whether they ; return an index or a color. In fact, other classes which ; use Coloring objects will likely specify which type of ; coloring they can be used with, so if at all possible, a ; class should be derived from one of these sub-classes ; rather than directly from Coloring. ;

; The one exception is for those very, very few colorings ; that can behave as either a gradient coloring or a direct ; coloring. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Coloring(Generic pparent) Generic.Generic(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz first value for the sequence; for a normal coloring formula, this will be #z ; @param ppixel seed value for the sequence; for a normal coloring formula, this will be #pixel func Init(complex pz, complex ppixel) ; Base class implementation just resets the iteration counter ; and clears the solid-color flag. m_Iterations = 0 m_Solid = false m_Pixel = ppixel endfunc ; Process the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. Note ; that such processing generally will not know in advance ; precisely how long the sequence is, and should be prepared ; to deal with sequences of arbitrary length. ;

; Your coloring may determine at some point that a solid color ; should be used rather than an index value. ; ; @param pz next value in the sequence; corresponds to #z in a coloring formula func Iterate(complex pz) ; Base class implementation only updates the iteration counter. ; Some colorings may not need to do per-iteration updates, but ; if their coloring depends on the final sequence value, they ; will need to save pz in a member variable so it can be used ; in the Result() function. m_Iterations = m_Iterations + 1 endfunc ; Produce a resulting color index after a sequence is finished ;

; This corresponds to the final: section in a coloring formula. ; Once it is called, no further calls to Iterate() should be ; made without calling Init() first. ; ; @return the color (corresponding to #color in a coloring formula) color func Result(complex pz) ; Base class implementation just returns black. ; Every coloring class should override this function. return rgb(0,0,0) endfunc ; Test whether the sequence produced a solid-color value rather than an index. ;

; This test is usually fairly complex and should be done as part of ; Iterate() or Result(), with the result saved in m_Solid. This function ; can then be left as the base class implementation, which just returns ; the value of that flag. ; ; @return whether the sequence produced a solid-color value bool func IsSolid() return m_Solid endfunc ; Test whether the formula produces gradient colors or direct colors. ;

; Note: this MUST NOT change over the lifetime of the object; other ; objects will call this and alter their behavior because of it, and ; it must give the same result regardless of the sequence. The only ; factors that should influence the return value of this function are ; parameters that the user has selected. ; ; @return whether the class provides gradient colors (true) or direct colors (false) bool func IsGradient() ; Base class implementation always declares itself as returning direct colors. return false endfunc ; Get saved value of pixel. ;

; Sometimes it's useful for another object to be able to retrieve the ; value of ppixel that's been passed in during Init(). This provides ; that access. ; ; @return the value of ppixel from the Init() for this sequence complex func GetPixel() return m_Pixel endfunc protected: int m_Iterations ; count the number of iterations bool m_Solid ; flag indicating whether sequence was solid or not complex m_Pixel ; saved value of pixel from Init() default: int param v_coloring caption = "Version (Coloring)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_coloring < 100 endparam } class GradientColoring(Coloring) { ; Gradient-based Coloring base class. ;

; This is a generic coloring class. Its purpose is to take ; a sequence of points and reduce it to a gradient index. ; This gradient index can then be assigned to #index or ; used with the gradient() function or in any other fashion. ;

; If you are migrating an existing coloring formula to a ; Coloring-derived class, the process is fairly straightforward. ; Any variables you set in your global: section that are used ; elsewhere in your formula should be declared in your protected: ; section. Move the code from your global: section into the ; constructor. Move the code from your init: section into the ; Init() function. Move the code from your loop: section into the ; Iterate() function. Move the code from your final: section ; into the Result() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func GradientColoring(Generic pparent) Coloring.Coloring(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz first value for the sequence; for a normal coloring formula, this will be #z ; @param ppixel seed value for the sequence; for a normal coloring formula, this will be #pixel func Init(complex pz, complex ppixel) Coloring.Init(pz, ppixel) endfunc ; Process the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. Note ; that such processing generally will not know in advance ; precisely how long the sequence is, and should be prepared ; to deal with sequences of arbitrary length. ;

; Your coloring may determine at some point that a solid color ; should be used rather than an index value. ; ; @param pz next value in the sequence; corresponds to #z in a coloring formula func Iterate(complex pz) ; Base class implementation only updates the iteration counter. ; Some colorings may not need to do per-iteration updates, but ; if their coloring depends on the final sequence value, they ; will need to save pz in a member variable so it can be used ; in the Result() function. m_Iterations = m_Iterations + 1 endfunc ; Produce a resulting color index after a sequence is finished ;

; This corresponds to the final: section in a coloring formula. ; Once it is called, no further calls to Iterate() should be ; made without calling Init() first. ; ; @return the color (corresponding to #color in a coloring formula) color func Result(complex pz) ; Base class implementation looks the color up in the gradient. ; Gradient-based coloring does not need to override this ; function unless different behavior when used as a direct ; coloring is required. return gradient(ResultIndex(pz)) endfunc ; Produce a resulting color index after a sequence is finished ;

; This corresponds to the final: section in a coloring formula. ; Once it is called, no further calls to Iterate() should be ; made without calling Init() first. ; ; @return the gradient index (corresponding to #index in a coloring formula) float func ResultIndex(complex pz) ; Base class implementation just returns the iteration counter. ; Every coloring class should override this function. return m_Iterations endfunc ; Test whether the sequence produced a solid-color value rather than an index. ;

; This test is usually fairly complex and should be done as part of ; Iterate() or Result(), with the result saved in m_Solid. This function ; can then be left as the base class implementation, which just returns ; the value of that flag. ; ; @return whether the sequence produced a solid-color value bool func IsSolid() return m_Solid endfunc ; Test whether the formula produces gradient colors or direct colors. ;

; Note: this MUST NOT change over the lifetime of the object; other ; objects will call this and alter their behavior because of it, and ; it must give the same result regardless of the sequence. The only ; factors that should influence the return value of this function are ; parameters that the user has selected. ; ; @return whether the class provides gradient colors (true) or direct colors (false) bool func IsGradient() ; Base class implementation always declares itself as returning gradient colors. return true endfunc protected: default: int param v_gradientcoloring caption = "Version (GradientColoring)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_gradientcoloring < 100 endparam } class DirectColoring(Coloring) { ; Direct Coloring base class. ;

; This is a generic coloring class. Its purpose is to take ; a sequence of points and reduce it to an explicit color. ; This gradient index can then be assigned to #color or ; or used in any other fashion. ;

; If you are migrating an existing coloring formula to a ; Coloring-derived class, the process is fairly straightforward. ; Any variables you set in your global: section that are used ; elsewhere in your formula should be declared in your protected: ; section. Move the code from your global: section into the ; constructor. Move the code from your init: section into the ; Init() function. Move the code from your loop: section into the ; Iterate() function. Move the code from your final: section ; into the Result() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func DirectColoring(Generic pparent) Coloring.Coloring(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz first value for the sequence; for a normal coloring formula, this will be #z ; @param ppixel seed value for the sequence; for a normal coloring formula, this will be #pixel func Init(complex pz, complex ppixel) Coloring.Init(pz, ppixel) endfunc ; Process the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. Note ; that such processing generally will not know in advance ; precisely how long the sequence is, and should be prepared ; to deal with sequences of arbitrary length. ;

; Your coloring may determine at some point that a solid color ; should be used rather than an index value. ; ; @param pz next value in the sequence; corresponds to #z in a coloring formula func Iterate(complex pz) m_Iterations = m_Iterations + 1 endfunc ; Produce a resulting color index after a sequence is finished ;

; This corresponds to the final: section in a coloring formula. ; Once it is called, no further calls to Iterate() should be ; made without calling Init() first. ; ; @return the gradient index (corresponding to #index in a coloring formula) color func Result(complex pz) return rgb(0,0,0) endfunc ; Test whether the sequence produced a solid-color value rather than an index. ;

; This test is usually fairly complex and should be done as part of ; Iterate() or Result(), with the result saved in m_Solid. This function ; can then be left as the base class implementation, which just returns ; the value of that flag. ; ; @return whether the sequence produced a solid-color value bool func IsSolid() return m_Solid endfunc protected: default: int param v_directcoloring caption = "Version (DirectColoring)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_directcoloring < 100 endparam } class IntegerTransfer(Generic) { ; This is a generic integer transfer class. It accepts one or more ; integer values in sequence and returns an integer value for each. ; Compare this to the Transform and Transfer base classes. public: ; constructor func IntegerTransfer(Generic pparent) Generic.Generic(pparent) endfunc ; call this to begin processing a sequence ; NOTE: although IntegerTransfer computes integers, it is still initialized ; with a complex value, usually corresponding to the point this ; is being applied to func Init(complex pz) m_Iterations = 0 endfunc ; call this to process another value in the sequence int func Iterate(int pn) m_Iterations = m_Iterations + 1 return pn endfunc ; Update internal counters without transforming a value ; func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc protected: int m_Iterations ; count the number of iterations default: int param v_integertransfer caption = "Version (IntegerTransfer)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_integertransfer < 100 endparam } class Transfer(Generic) { ; This is a generic transfer class. It accepts one or more ; real values in sequence and returns a real value for each. ; Compare this to the Transform base class. public: ; constructor func Transfer(Generic pparent) Generic.Generic(pparent) endfunc ; call this to begin processing a sequence ; NOTE: although Transfer computes reals, it is still initialized ; with a complex value, usually corresponding to the point this ; is being applied to func Init(complex pz) m_Iterations = 0 endfunc ; call this to process another value in the sequence float func Iterate(float pr) m_Iterations = m_Iterations + 1 return pr endfunc ; Update internal counters without transforming a value ; func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc protected: int m_Iterations ; count the number of iterations default: int param v_transfer caption = "Version (Transfer)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_transfer < 100 endparam } class ColorTransfer(Generic) { ; This is a generic color transfer class. It accepts one or more ; color values in sequence and returns a color value for each. ; Compare this to the Transform and Transfer base classes. public: ; constructor func ColorTransfer(Generic pparent) Generic.Generic(pparent) endfunc ; call this to begin processing a sequence ; NOTE: although ColorTransfer computes colors, it is still initialized ; with a complex value, usually corresponding to the point this ; is being applied to func Init(complex pz) m_Iterations = 0 endfunc ; call this to process another value in the sequence color func Iterate(color pcolor) m_Iterations = m_Iterations + 1 return pcolor endfunc ; Update internal counters without transforming a value ; func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc protected: int m_Iterations ; count the number of iterations default: int param v_colortransfer caption = "Version (ColorTransfer)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_colortransfer < 100 endparam } class ColorMerge(Generic) { ; This is a generic color blending class. It accepts two ; colors and produces a merged result. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func ColorMerge(Generic pparent) Generic.Generic(pparent) endfunc ; Color merging function ; ; @param pbottom color on the bottom ; @param ptop color on the top ; @return merged color; alpha value should be from the top color or composing won't work right color func Merge(color pbottom, color ptop) return ptop endfunc ; Opaque test function ;

; In some cases, formulas want to know if, given a particular ; color, merging it onto another color will leave any of the ; other color visible (i.e. does all the color come from the ; top color). This function should return true in that case. ; The default implementation always returns false; if you ; write a merging class which can identify opaque colors, you ; should override this function. (Note this isn't as simple ; as just looking at the top color's opacity; in non-Normal ; merge modes the bottom color still influences the output.) ; ; @param ptop top color to test ; @return true if top color is fully opaque bool func IsOpaque(color ptop) return false endfunc ; Static composing/blending helper function ;

; The explicit formula provided in the UF help for doing a ; color blend is: ;

; compose(b, blend(t, mergeX(b, t), alpha(b)), o) ;

; This at first seems overly complicated, but in fact is ; necessary to account for alpha in both the bottom and ; top colors. This function wraps up the compose and blend ; steps so you don't have to remember the correct formula. ; To use this, you will need to have the top color both ; before and after the merge function has been applied. ;

; This is a static function, so you can call it without ; creating a ColorMerge object. Just call ColorMerge.Stack() ; directly. ; ; @param pbottom color on the bottom ; @param ptop color on the top (pre-merge) ; @param pmerged color on the top (post-merge) ; @param popacity opacity; 0 indicates all bottom, 1 indicates all top static color func Stack(color pbottom, color ptop, color pmerged, float popacity) return compose(pbottom, blend(ptop, pmerged, alpha(pbottom)), popacity) endfunc ; Composing/blending helper function ;

; Like Stack(), this is a helper function to make it easier ; to merge colors layer-style. However, Stack() still requires ; you to produce the merged color. FullMerge() will do the ; merging for you, but you have to use a ColorMerge object ; (not just call a static function). ; ; @param pbottom color on the bottom ; @param ptop color on the top (pre-merge) ; @param popacity opacity; 0 indicates all bottom, 1 indicates all top color func FullMerge(color pbottom, color ptop, float popacity) return Stack(pbottom, ptop, Merge(pbottom, ptop), popacity) endfunc default: int param v_colormerge caption = "Version (ColorMerge)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_colormerge < 100 endparam } class Generator(Generic) { ; Generator base class. ;

; This is a generic generator class. Its purpose is to take ; an arbitrary initial value and produce a sequence of ; values of indeterminate length. This roughly corresponds ; to UF's native fractal formula type, except that it starts ; with and produces real values, not complex values. Compare ; this to the Formula base class. ;

; Note that by convention, values produced from a Generator ; class should be in the range [0,1]. If you adhere to this ; convention, other classes that use Generator objects will ; know what range to expect. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Generator(Generic pparent) Generic.Generic(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence ; @return first value in the sequence float func Init(float pz) ; Base class implementation flags the sequence to end ; immediately. As long as you provide your own implementation ; of IsBailedOut(), this is irrelevant. We also clear the ; iteration counter. m_Iterations = 0 m_BailedOut = true return pz endfunc ; Set up for a sequence of values ;

; This alternate form of Init will cause the Generator to ; start with its "default" sequence. This can be used if a ; generator allows a user to select a "seed" value itself. ; Such Generator classes should override this function. ; ; @return first value in the sequence float func InitDefault() return Init(0) endfunc ; Produce the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. ; ; @param pz previous value in the sequence. Note that you should always use this value for computing the next iteration, rather than a saved value, as the calling code may modify the returned value before passing it back to the next Iterate() call. ; @return the next value in the sequence float func Iterate(float pz) ; Base class implentation automatically updates the ; iteration counter and returns an unmodified pz ; (that is, the default sequence is the same value ; repeated infinitely). m_Iterations = m_Iterations + 1 return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; This is typically called once per iteration to test whether ; the sequence has completed. If producing this result is ; difficult to do separately from the Iterate() processing, ; you should produce the result in Iterate() and set the ; m_BailedOut flag appropriately, and leave IsBailedOut() ; unimplemented in your class to inherit this behavior. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Generator classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(float pz) return m_BailedOut endfunc protected: int m_Iterations ; count the number of iterations bool m_BailedOut ; flag indicating whether sequence has bailed out or not default: int param v_generator caption = "Version (Generator)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_generator < 100 endparam } class IntegerGenerator(Generic) { ; Integer Generator base class. ;

; This is a generic generator class. Its purpose is to take ; an arbitrary initial value and produce a sequence of ; values of indeterminate length. This roughly corresponds ; to UF's native fractal formula type, except that it starts ; with and produces integer values, not complex values. Compare ; this to the Formula and Generator base classes. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func IntegerGenerator(Generic pparent) Generic.Generic(pparent) endfunc ; Set up for a sequence of values ;

; This function will be called at the beginning of each ; sequence of values (e.g. at the beginning of each fractal ; orbit). ; ; @param pz seed value for the sequence ; @return first value in the sequence int func Init(int pz) ; Base class implementation flags the sequence to end ; immediately. As long as you provide your own implementation ; of IsBailedOut(), this is irrelevant. We also clear the ; iteration counter. m_Iterations = 0 m_BailedOut = true return pz endfunc ; Set up for a sequence of values ;

; This alternate form of Init will cause the Generator to ; start with its "default" sequence. This can be used if a ; generator allows a user to select a "seed" value itself. ; Such Generator classes should override this function. ; ; @return first value in the sequence int func InitDefault() return Init(0) endfunc ; Produce the next value in the sequence ;

; As long as the sequence has not bailed out, this function ; will be continually called to produce sequence values. ; ; @param pz previous value in the sequence. Note that you should always use this value for computing the next iteration, rather than a saved value, as the calling code may modify the returned value before passing it back to the next Iterate() call. ; @return the next value in the sequence int func Iterate(int pz) ; Base class implentation automatically updates the ; iteration counter and returns an unmodified pz ; (that is, the default sequence is the same value ; repeated infinitely). m_Iterations = m_Iterations + 1 return pz endfunc ; Test whether the formula has bailed out (i.e. the sequence is complete) ;

; This is typically called once per iteration to test whether ; the sequence has completed. If producing this result is ; difficult to do separately from the Iterate() processing, ; you should produce the result in Iterate() and set the ; m_BailedOut flag appropriately, and leave IsBailedOut() ; unimplemented in your class to inherit this behavior. ; ; @param pz last sequence value to test; this should be the value returned from the previous Iterate() call. Note that it is acceptable to ignore pz and use m_BailedOut, but any code calling IsBailedOut() should pass in the correct pz for Generator classes which do not use m_BailedOut. ; @return true if the sequence has bailed out (i.e. should be terminated) bool func IsBailedOut(int pz) return m_BailedOut endfunc protected: int m_Iterations ; count the number of iterations bool m_BailedOut ; flag indicating whether sequence has bailed out or not default: int param v_integergenerator caption = "Version (IntegerGenerator)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_integergenerator < 100 endparam } class GradientWrapper(Generic) { ; GradientWrapper base class. ;

; This class has a simple job: to convert a real index ; value to a color. The base implementation does this ; using UF's native gradient() function. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func GradientWrapper(Generic pparent) Generic.Generic(pparent) endfunc ; Lookup function for code that does not support ; multiple channels; by default, use channel 0 color func getColor(float pindex) return getColorChannel(pindex, 0) endfunc ; Lookup function supporting multiple channels ; By default, we use the UF gradient color func getColorChannel(float pindex, int pchannel) return gradient(pindex) endfunc default: int param v_gradientwrapper caption = "Version (GradientWrapper)" default = 101 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_gradientwrapper < 100 endparam } class ImageWrapper(Generic) { ; ImageWrapper base class. ;

; This class provides access to UF's native Image object ; in an extensible class. If you use an Image object ; directly, you allow no opportunity for other Image-like ; classes to be used with your code. If instead you use ; an ImageWrapper object, you will get the Image features ; in the default case and still allow Image-like classes ; to be used with your code. public: ; constructor func ImageWrapper(Generic pparent) Generic.Generic(pparent) endfunc ; member functions ; these are pass-throughs to the Image class color func getColor(complex pz) return rgba(0,0,0,0) endfunc color func getPixel(int px, int py) return rgba(0,0,0,0) endfunc bool func getEmpty() return true endfunc int func getWidth() return 0 endfunc int func getHeight() return 0 endfunc ; this is not a pass-through; it aspect-corrects a pixel to the ; correct range. UF's default Image object maps the entire image ; into the (-1,-1) .. (1,1) range regardless of size; this will ; shrink the real or imaginary component of pz as appropriate to ; restore the original image proportions. ; Note: UF's default Image coloring algorithm does this as well. complex func NormalizePixel(complex pz) if (this.getEmpty()) return pz elseif (this.getWidth() > this.getHeight()) return real(pz) + flip(imag(pz) * this.getWidth() / this.getHeight()) else return real(pz) * this.getHeight() / this.getWidth() + flip(imag(pz)) endif endfunc protected: default: int param v_imagewrapper caption = "Version (ImageWrapper)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_imagewrapper < 100 endparam } ; ----------------------------------------------------------------------------------------------------------------------------- ; Specialized Base classes class TrapShape(Generic) { ; Trap Shape base class. ;

; This is a generic trap shape class. The purpose of a trap shape ; is very simple: to take one or more complex values and determine ; a "distance" to a trap shape for those values. There are some ; extra features provided to facilitate many of the things that ; have been done in the past with trap shapes and trapping ; algorithms. public: ; constructor func TrapShape(Generic pparent) Generic.Generic(pparent) endfunc ; call this before each sequence of values to be trapped func Init(complex pz) m_Iterations = 0 m_LastZ = pz m_LastChannel = 0 endfunc ; call this for each iteration being trapped float func Iterate(complex pz) m_Iterations = m_Iterations + 1 m_LastZ = pz return real(pz) endfunc ; Update internal counters without transforming a value ; func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc ; call this to set the threshold selected by the user ; it's up to the calling code to apply the threshold, ; but some trap shapes or merges need to know the ; threshold to work properly func SetThreshold(float pthreshold) m_Threshold = pthreshold endfunc ; Get transformed point. ;

; Some trap modes and trap coloring modes require ; access to the "transformed" point. This is the ; final point that is actually used by the trap shape ; class to test against the trap. For most trap shape ; classes this will be exactly the same as the pz ; passed into Iterate(), and those classes can use ; this base class functionality. However, some trap ; shape classes (e.g. TrapShapeMerge, TrapShapeBlock) ; perform some internal transformations on pz before ; passing it to the trap shape; this gives them the ; opportunity to call their trap shape's function to ; get the final transformed point. ; ; @return the last transformed point used (i.e. call GetTransformedPoint() after calling Iterate(), not before) complex func GetTransformedPoint() return m_LastZ endfunc ; Get texture value. ;

; Ordinarily, a trap shape does not have a native texture ; (it is flat). Some formulas may pair trap shapes with ; trap textures, and may need more information about trap ; textures when trap shapes are nested via TrapShapeMerge. ; This function provides support for accessing texture ; information. ; ; @return the texture value for the last point used (i.e. call GetTextureValue() after calling Iterate(), not before) float func GetTextureValue() return 0 endfunc ; Get color channel. ;

; Some trap shapes may be inherently multi-colored. Such ; trap shapes may either override this function to return ; a different color index, or may simply store the index ; in m_LastChannel and this function will return it. The ; color index is just an indication of which color channel ; (in a GradientWrapper object) should be used. ; ; @return the color channel value for the last point used (i.e. call GetColorChannel() after calling Iterate(), not before) int func GetColorChannel() return m_LastChannel endfunc protected: int m_Iterations float m_Threshold complex m_LastZ int m_LastChannel default: int param v_trapshape caption = "Version (TrapShape)" default = 101 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapshape < 100 endparam } class ColorTrap(Generic) { ; This is a generic color trap formula. It takes a complex ; point and returns a color. public: ; constructor func ColorTrap(Generic pparent) Generic.Generic(pparent) endfunc ; call this before each sequence of values to be trapped func Init(complex pz) m_Iterations = 0 endfunc ; call this to produce a result ; you should always override this function color func Iterate(complex pz) m_Iterations = m_Iterations + 1 return gradient(cabs(pz)) endfunc ; Update internal counters without transforming a value ; func IterateSilent() ; Perform the same iteration-count update as Iterate(), but ; don't return a value. See the comments in Iterate() for ; more information on when to call this. m_Iterations = m_Iterations + 1 endfunc protected: int m_Iterations default: int param v_colortrap caption = "Version (ColorTrap)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_colortrap < 100 endparam } class TrapMode(Generic) { ; This is a generic trap mode formula. It takes a sequence of ; complex point pairs (untransformed and transformed) and ; distances and produces a final result public: ; constructor func TrapMode(Generic pparent) Generic.Generic(pparent) endfunc ; call this at the beginning of each sequence func Init(complex pz) m_Iterations = 0 m_Solid = false int j = 0 while (j < 4) m_UntransformedPoints[j] = (0,0) m_TransformedPoints[j] = (0,0) m_IterationPoints[j] = 0 m_Distances[j] = 0.0 m_Textures[j] = 0.0 j = j + 1 endwhile endfunc ; call this for each point func Iterate(complex pz, complex pzt, float pdistance, float ptexture) m_Iterations = m_Iterations + 1 endfunc ; call this for each point that is ignored func IterateSilent() m_Iterations = m_Iterations + 1 endfunc ; call this to compute final results func Result() endfunc ; call this to determine if the last sequence produced a solid point bool func IsSolid() return m_Solid endfunc ; get a final untransformed point complex func GetUntransformedPoint(int pindex) return m_UntransformedPoints[pindex] endfunc ; get a final transformed point complex func GetTransformedPoint(int pindex) return m_TransformedPoints[pindex] endfunc ; get a final distance float func GetDistance(int pindex) return m_Distances[pindex] endfunc ; get a final texture float func GetTexture(int pindex) return m_Textures[pindex] endfunc ; get a final iteration float func GetIteration(int pindex) return m_IterationPoints[pindex] endfunc ; get threshold value float func GetThreshold() return 1.0 endfunc func SetThreshold(float pthreshold) endfunc ; get whether threshold is even used bool func UsesThreshold() return false endfunc protected: int m_Iterations ; count the number of iterations bool m_Solid ; flag indicating whether sequence was solid or not complex m_UntransformedPoints[4] complex m_TransformedPoints[4] float m_IterationPoints[4] float m_Distances[4] float m_Textures[4] default: int param v_trapmode caption = "Version (TrapMode)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapmode < 100 endparam } class TrapModeWithThreshold(TrapMode) { ; Variant of TrapMode that is thresholdable public: func TrapModeWithThreshold(Generic pparent) TrapMode.TrapMode(pparent) m_Threshold = 1.0 endfunc float func GetThreshold() return m_Threshold endfunc func SetThreshold(float pthreshold) m_Threshold = pthreshold endfunc bool func UsesThreshold() return true endfunc protected: float m_Threshold default: int param v_trapmodewiththreshold caption = "Version (TrapModeWithThreshold)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapmodewiththreshold < 100 endparam } class TrapColoring(Generic) { ; This is a generic trap coloring mode formula. It takes ; the results of a trap mode and returns a color index public: ; constructor func TrapColoring(Generic pparent) Generic.Generic(pparent) endfunc ; get a final result float func Result(TrapMode ptrapmode) return ptrapmode.GetDistance(0) endfunc default: int param v_trapcoloring caption = "Version (TrapColoring)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapcoloring < 100 endparam } class UtilityTransform(Transform) { ; Utility Transform base class. ;

; This class adds no new functionality to the Transform class. ; Its purpose is to provide a grouping for utility transforms ; that other classes may use internally to provide extra ; features. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func UtilityTransform(Generic pparent) Transform.Transform(pparent) endfunc default: int param v_utilitytransform caption = "Version (UtilityTransform)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_utilitytransform < 100 endparam } class ClipShape(Transform) { ; This is a generic clip shape class. It takes an arbitrary ; point and determines whether it is inside (solid) or ; outside (not solid) (and returns the original point). ;

; This differs from a more generic Transform in that it has ; additional functions to allow it to manipulate a Handles ; object, and allows an aspect ratio to be passed in from ; a ScreenRelative transform and a generalized rotation. public: func ClipShape(Generic pparent) Transform.Transform(pparent) m_Aspect = 1.0 endfunc func SetHandles(Handles phandles) m_Handles = phandles endfunc func SetAspect(float paspect) m_Aspect = paspect endfunc protected: Handles m_Handles float m_Aspect default: int param v_clipshape caption = "Version (Clipshape)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_clipshape < 100 endparam } class MappingShape(Generic) { ; MappingShape base class. ;

; MappingShape represents a shape in 3D space with a texture (usually a ; fractal) mapped on it. It provides the GetTextureCoordinates function which ; can be called to obtain the intersection point of the shape and a given line, ; and the texture coordinates at that point. ;

; MappingShape uses the left-handed coordinate system, where X points to the ; right, Y points up, and Z points into the screen. To determine the direction ; of a rotation, hold your left hand such that your thumb points into the ; direction of the axis around which the rotation takes place, and then your ; curled fingers point into the direction of positive rotation. public: ; GetTextureCoordinates tests a ray against the mapping shape. If the ray ; intersects the shape, it returns true, and returns the coordinates of the ; point on the shape that intersects the ray in the result argument. If the ; shape does not intersect the ray, returns false. ;

; A concrete shape descendant should override this function to define the ; mapping shape. ;

; The ray is defined as:
; [ sx ] [ dx ]
; [ sy ] + lambda [ dy ]

; This is a utility transform used to display "handles". ; Handles are visual aids that show where important points ; that influence a formula are located within an image. ; Provided you have defined the handle data, it will return a ; solid flag for each transformed point. It is up to your ; code to decide what to do with that information. As an ; example, in a transform you might toggle the state of the ; #solid flag based on the handle rendering result; in a ; coloring formula, you might use #solid or you might use ; a specific color. ;

; This class provides three types of visual aids. First are ; the handles themselves; these are located at a specific ; point, and can be a box, a circle, a diamond, an X, a ; point with an arrow, or a point with a double arrow. Each ; handle may also be labeled with a number (integer), a ; combination of letters and numbers (up to five), or an ; 8x8 pixel icon. Handles may also be rotated to arbitrary ; angles (especially useful for arrows). See SetHandlePoint() ; for more information. ;

; The second type of visual aid are lines. Often these lines ; are used to connect handle points together, but there is ; absolutely no requirement that this be done; handle lines ; can be constructed anywhere a formula needs. Lines may ; be solid, dotted, short dashed, or long dashed. See ; SetHandleLine() for more information. ;

; The third type of visual aid are circles. Circles may be ; placed anywhere and be of any size; as with lines, they ; may be solid, dotted, short dashed, or long dashed. Note ; that circles will always be circles, even if the aspect ; ratio of the fractal image is not 1:1. ;

; Handles are intended as visual aids for the creation of ; fractal images, not as part of the final created image. ; For this reason, the default settings for the Handles ; class will automatically disable handle rendering for ; any rendered image. Also, handle sizes are always measured ; in pixels, which means when the image size shrinks (as it ; does for previews) the handles do not shrink. This makes ; them too large for most uses of the preview window, so ; the default settings for Handles disables them on preview ; images as well. The user can override these settings if ; they choose to. ;

; To use the Handles class, create a Handles parameter for ; your formula. (You probably do not want it to be selectable.) ; Then, call SetHandleCount() to indicate the number of points, ; lines, and circles you will use. Call SetHandlePoint(), ; SetHandleLine(), and SetHandleCircle() for each element. ; Then, while processing, use Init() to set up each pixel and ; Iterate() to compute the handle solid flag value. Use ; IsSolid() to test the result of the handle rendering. public: ; Constructor ; ; @param pparent a reference to the object creating the new object; typically, 'this' func Handles(Generic pparent) UtilityTransform.UtilityTransform(pparent) ; pixel coordinate transform m_PixelCoordinate = new PixelCoordinate(pparent) ; digit/letter bitmaps m_DigitsDetail[0] = (409387979.0,3552798232.0) ; 0 m_DigitsDetail[1] = (406329368.0,404232252.0) ; 1 m_DigitsDetail[2] = (1011221251.0,102261247.0) ; 2 m_DigitsDetail[3] = (1011221262.0,50546236.0) ; 3 m_DigitsDetail[4] = (271605860.0,3439266828.0) ; 4 m_DigitsDetail[5] = (2126576646.0,50562684.0) ; 5 m_DigitsDetail[6] = (476110054.0,3284362812.0) ; 6 m_DigitsDetail[7] = (4270197772.0,202905624.0) ; 7 m_DigitsDetail[8] = (1013335612.0,1656972860.0) ; 8 m_DigitsDetail[9] = (1013367747.0,1731921464.0) ; 9 m_DigitsDetail[10] = (0.0,0.0) ; blank m_DigitsDetail[11] = (404237356.0,1182696323.0) ; A m_DigitsDetail[12] = (4240885500.0,3334719228.0) ; B m_DigitsDetail[13] = (1013170368.0,3233833788.0) ; C m_DigitsDetail[14] = (4173775811.0,3284387576.0) ; D m_DigitsDetail[15] = (4274045176.0,3233857791.0) ; E m_DigitsDetail[16] = (4290822392.0,3233857728.0) ; F m_DigitsDetail[17] = (1013170383.0,3284362044.0) ; G m_DigitsDetail[18] = (3284386815.0,3284386755.0) ; H m_DigitsDetail[19] = (1008211992.0,404232252.0) ; I m_DigitsDetail[20] = (252052998.0,101092472.0) ; J m_DigitsDetail[21] = (3335313648.0,3637298883.0) ; K m_DigitsDetail[22] = (3233857728.0,3233858047.0) ; L m_DigitsDetail[23] = (2177099775.0,3687039939.0) ; M m_DigitsDetail[24] = (2210653171.0,3687827395.0) ; N m_DigitsDetail[25] = (1013367747.0,3284362812.0) ; O m_DigitsDetail[26] = (4240884678.0,4240490688.0) ; P m_DigitsDetail[27] = (1013367747.0,3553322555.0) ; Q m_DigitsDetail[28] = (4240884678.0,4241278659.0) ; R m_DigitsDetail[29] = (1013080124.0,100894332.0) ; S m_DigitsDetail[30] = (4279769112.0,404232252.0) ; T m_DigitsDetail[31] = (3284386755.0,3284362812.0) ; U m_DigitsDetail[32] = (2206418502.0,741087256.0) ; V m_DigitsDetail[33] = (2206434118.0,1448487980.0) ; W m_DigitsDetail[34] = (2210819128.0,946652867.0) ; X m_DigitsDetail[35] = (3284362812.0,404232216.0) ; Y m_DigitsDetail[36] = (4286975000.0,405823999.0) ; Z m_DigitsDetail[37] = (417038043.0,402653184.0) ; * m_DigitsDetail[38] = (101059596.0,404238384.0) ; / m_DigitsDetail[39] = (0.0,6168.0) ; . endfunc ; Compute handle state for a single point within a sequence ;

; Handles is a Transform, so to test the handle state for ; a particular point, the standard approach is to Init(), ; Iterate() (just once), and then test IsSolid(). ; ; @param pz the complex value to test ; @return the same complex value (Handles is a transform but the value is in the solid test; the complex value is not transformed) complex func Iterate(complex pz) UtilityTransform.Iterate(pz) float sx float sy float vx float vy float dx float dy bool nolines = false ; note that we completely ignore pz; handles are NOT transformed along with ; the underlying image ; if handles are turned off, stop now if (!@p_showhandles || \ (#calculationPurpose == 1 && @p_handlesnopreview) || \ (#calculationPurpose == 3 && @p_handlesnopreview) || \ (#calculationPurpose == 2 && @p_handlesnorender)) return pz endif \$ifdef debug int testx = 400 int testy = 600 \$endif ; see if any handle points apply int j = 0 int k = 0 int l = length(m_HandleTypes) int label = 0 int base = 0 int digit = 0 float n = 0 while (j < l) if (m_HandleTypes[j] > 0 || m_HandleLabels[j] > 0) sx = real(m_HandleCenters[j]) sy = imag(m_HandleCenters[j]) vx = #x-sx vy = #y-sy if (m_HandleRotations[j] != 0) complex v = (vx + flip(vy)) * m_HandleRotations[j] vx = real(v) vy = imag(v) endif dx = abs(vx) dy = abs(vy) ; if we're inside any handle area, flag that no lines be drawn if ((m_HandleTypes[j] < 5 && dx < @p_handlesize+@p_handlewidth*2 && dy < @p_handlesize+@p_handlewidth*2) || \ (m_HandleTypes[j] >= 5 && dx < @p_handlewidth*2.5 && dy < @p_handlewidth*2.5)) nolines = true endif ; handle if ((m_HandleTypes[j] == 1 && \ (dx >= @p_handlesize || dy >= @p_handlesize) && \ (dx < @p_handlesize+@p_handlewidth && dy < @p_handlesize+@p_handlewidth)) || \ (m_HandleTypes[j] == 2 && \ dx*dx+dy*dy >= @p_handlesize*@p_handlesize && dx*dx+dy*dy < sqr(@p_handlesize+@p_handlewidth)) || \ (m_HandleTypes[j] == 3 && \ dx+dy >= @p_handlesize && dx+dy < @p_handlesize+@p_handlewidth) || \ (m_HandleTypes[j] == 4 && \ 2*abs(dx+0.5-dy) < @p_handlewidth && (dx < @p_handlesize+@p_handlewidth && dy < @p_handlesize+@p_handlewidth)) || \ (m_HandleTypes[j] == 5 && \ ((dy < @p_handlewidth*0.5 && vx >= 0 && vx < @p_handlesize+@p_handlewidth) || \ (vx >= @p_handlesize*0.75 && dx+dy >= @p_handlesize && dx+dy < @p_handlesize+@p_handlewidth) || \ (dx*dx+dy*dy <= 1.5*@p_handlewidth*@p_handlewidth))) || \ (m_HandleTypes[j] == 6 && \ ((dy < @p_handlewidth*0.5 && dx < @p_handlesize+@p_handlewidth) || \ (dx >= @p_handlesize*0.75 && dx+dy >= @p_handlesize && dx+dy < @p_handlesize+@p_handlewidth) || \ (dx*dx+dy*dy <= 1.5*@p_handlewidth*@p_handlewidth))) \ ) m_Solid = !m_Solid endif ; number if (@p_handleswithnumbers && \ m_HandleLabelSizes[j] > 0 && \ #x >= sx+@p_handlesize+@p_handlewidth && #x < sx+@p_handlesize+@p_handlewidth*(1+5*m_HandleLabelSizes[j]) && \ #y >= sy+@p_handlesize+@p_handlewidth && #y < sy+@p_handlesize+@p_handlewidth*5) nolines = true ; no lines through number dx = floor((#x-sx-@p_handlesize-@p_handlewidth)/@p_handlewidth*2) ; pixel coordinate within label area dy = floor((#y-sy-@p_handlesize-@p_handlewidth)/@p_handlewidth*2) label = m_HandleLabels[j] if (label > 0) base = 10 else base = 40 endif digit = m_HandleLabelSizes[j]-floor(dx/10)-1 ; digit position within label dx = dx % 10 ; pixel coordinate within digit k = floor(abs(label) / base^digit + 0.00001) % base ; actual glyph number (includes a bit of round-off fudge) if (dx > 0 && dx < 9) ; pixel lies within glyph boundary if (dy < 4) ; top half of glyph n = real(m_DigitsDetail[k]) ; bit values dx = (8-dx) + 8*(3-dy) ; bit number else ; bottom half of glyph n = imag(m_DigitsDetail[k]) ; bit values dx = (8-dx) + 8*(7-dy) ; bit number endif dy = 2^dx ; mask bit if (floor(n / dy) % 2 > 0) ; test bit m_Solid = !m_Solid endif endif endif ; icon ; note the same flag that controls number rendering controls icons if (@p_handleswithnumbers && \ (real(m_HandleIcons[j]) > 0 || imag(m_HandleIcons[j]) > 0) && \ #x >= sx+@p_handlesize+@p_handlewidth && #x < sx+@p_handlesize+@p_handlewidth*5 && \ #y >= sy+@p_handlesize+@p_handlewidth && #y < sy+@p_handlesize+@p_handlewidth*5) nolines = true ; no lines through number dx = floor((#x-sx-@p_handlesize-@p_handlewidth)/@p_handlewidth*2) ; pixel coordinate within icon area dy = floor((#y-sy-@p_handlesize-@p_handlewidth)/@p_handlewidth*2) if (dx >= 0 && dx < 8) ; pixel lies within glyph boundary if (dy < 4) ; top half of glyph n = real(m_HandleIcons[j]) ; bit values dx = (7-dx) + 8*(3-dy) ; bit number else ; bottom half of glyph n = imag(m_HandleIcons[j]) ; bit values dx = (7-dx) + 8*(7-dy) ; bit number endif dy = 2^dx ; mask bit if (floor(n / dy) % 2 > 0) ; test bit m_Solid = !m_Solid endif endif endif endif j = j + 1 endwhile ; see if any handle lines apply float linelength = 0 j = 0 l = length(m_LineTypes) while (j < l && !nolines && @p_handleswithlines) if (!m_LineConnects[j]) ; this line does not connect to the previous line linelength = 0 ; reset the accumulated line length (helps line up dashes/dots) endif if (m_LineTypes[j] > 0) dx = ((real(m_LineEnds[j]) - real(m_LineStarts[j])) * (#y - imag(m_LineStarts[j])) - (#x - real(m_LineStarts[j])) * (imag(m_LineEnds[j]) - imag(m_LineStarts[j]))) / m_LineDistances[j] dy = ((imag(m_LineEnds[j]) - imag(m_LineStarts[j])) * (#y - imag(m_LineStarts[j])) - (#x - real(m_LineStarts[j])) * (real(m_LineStarts[j]) - real(m_LineEnds[j]))) / m_LineDistances[j] if (abs(dx) < @p_handlewidth*0.5 && dy > 0 && dy < m_LineDistances[j]) if (m_LineTypes[j] == 1) m_Solid = !m_Solid elseif (m_LineTypes[j] == 2) if (dy % (@p_handlewidth*2) > @p_handlewidth) m_Solid = !m_Solid endif elseif (m_LineTypes[j] == 3) if (dy % (@p_handlesize) > @p_handlesize*0.5) m_Solid = !m_Solid endif elseif (m_LineTypes[j] == 4) if (dy % (@p_handlesize*2) > @p_handlesize) m_Solid = !m_Solid endif endif endif endif j = j + 1 endwhile linelength = linelength + 1 ; **** silences compiler warning ; see if any handle circles apply j = 0 l = length(m_CircleTypes) while (j < l && !nolines && @p_handleswithlines) if (m_CircleTypes[j] > 0) dx = cabs((#x + flip(#y)) - m_CircleCenters[j]) - m_CircleRadii[j] dy = atan2((#x + flip(#y)) - m_CircleCenters[j]) if (dy < 0) dy = dy + 2*#pi endif dy = dy * m_CircleRadii[j] \$ifdef debug if (#x == testx && #y == testy) print("test: ", #x, ",", #y) print("circle ", j, " center: ", m_CircleCenters[j], " radius: ", m_CircleRadii[j]) print("dx: ", dx, " dy: ", dy) endif \$endif if (abs(dx) < @p_handlewidth*0.5) if (m_CircleTypes[j] == 1) m_Solid = !m_Solid elseif (m_CircleTypes[j] == 2) if (dy % (@p_handlewidth*2) > @p_handlewidth) m_Solid = !m_Solid endif elseif (m_CircleTypes[j] == 3) if (dy % (@p_handlesize) > @p_handlesize*0.5) m_Solid = !m_Solid endif elseif (m_CircleTypes[j] == 4) if (dy % (@p_handlesize*2) > @p_handlesize) m_Solid = !m_Solid endif endif endif endif j = j + 1 endwhile return pz endfunc ; Set the number of handles to render ;

; You must call this function before you can define any handles. ; Specify the maximum number of each type of visual aid you ; will use. Any handles not defined will not be rendered, but ; defining more than you need will slow down rendering slightly. ; ; @param phandles number of handle points ; @param plines number of handle lines ; @param pcircles number of handle circles ; @param pcurves number of special handles; this is being reserved for future use, so you should always use 0 func SetHandleCount(int phandles, int plines, int pcircles, int pcurves) if (phandles >= 0) setLength(m_HandleTypes, phandles) setLength(m_HandleLabels, phandles) setLength(m_HandleLabelSizes, phandles) setLength(m_HandleIcons, phandles) setLength(m_HandleCenters, phandles) setLength(m_HandleRotations, phandles) endif if (plines >= 0) setLength(m_LineTypes, plines) setLength(m_LineConnects, plines) setLength(m_LineDistances, plines) setLength(m_LineStarts, plines) setLength(m_LineEnds, plines) endif if (pcircles >= 0) setLength(m_CircleTypes, pcircles) setLength(m_CircleCenters, pcircles) setLength(m_CircleRadii, pcircles) endif endfunc ; Define a handle point ;

; This function defines all the information about a handle ; point. You will call this once for each handle point. ;

; Handles can be one of several pre-defined shapes. You may ; specify a blank handle, which is useful if you simply want ; to place a label within your image. Squares are, by ; convention, used to indicate points that are on an edge ; that the user can directly select. Circles are used to ; indicate points that are NOT on an edge that the user may ; select. Diamonds are used to indicate interior or center ; points; the user may or may not be allowed to directly ; select this point. X should be used for an implied or ; computed point that the user CANNOT select. Arrows should ; be used for points where a particular direction needs to ; be indicated (specify the angle as something other than ; zero). Double-arrows should be used for points where the ; location along the line perpendicular to the arrows is ; irrelevant, and only the location in the direction the ; arrows show is relevant; they may also be used to mark ; points on a handle circle. A single triangle is used to ; mark a point on a line that the user may set indirectly. ; A double triangle should be used to indicate a point on ; a line where the line is set automatically but the ; position along the line is under the user's control ; (different from the double-arrow, where the user sets the ; line position; here they are setting the position of a ; point along a line). ;

; Any handle shape may be rotated, but rotation is primarily ; intended for arrows. Rotation does not affect labels or ; icons, only the handle shape. ;

; Labels can be numbers or short letter/number combinations. ; For positive integers, simply provide the number. For ; words, the process is a bit more complicated. Words may ; be five characters or less (including spaces). Each character ; is encoded in a base-40 encoding scheme: ;

; 0 to 9: digits 0 to 9
; 10: space
; 11 to 36: letters A to Z
; 37: asterisk (*)
; 38: slash (/)
; 39: period (.)
;

; As an example, to encode the word "TEST", convert each letter ; to its number: 30 15 29 30 ;

; Next, treat these as base 40 values: ; (((30*40+15)*40+29)*40+30) = 1945190 ;

; Then, pass in the value as a negative integer to indicate it ; should be interpreted as characters rather than just a number: ; -1945190 ;

; Icons may be used in place of letters. (If you place an icon ; and characters on the same handle point, they will overlap.) ; Icons are 8x8 pixel bitmaps encoded as a single complex value: ;

; row 1: real(icon) bits 31 to 24
; row 2: real(icon) bits 23 to 16
; row 3: real(icon) bits 15 to 8
; row 4: real(icon) bits 7 to 0
; row 5: imag(icon) bits 31 to 24
; row 6: imag(icon) bits 23 to 16
; row 7: imag(icon) bits 15 to 8
; row 8: imag(icon) bits 7 to 0 ;

; Because 32 bits are used for each icon, you must specify the ; values as explicitly real (include .0 on the end). Note that ; internally, the letters, numbers, and symbols are specified ; with this same format. ; ; @param phandle slot number, starting at 0; if you indicate -1, the next slot will be used automatically ; @param ptype handle type: 0 = blank (just number), 1 = square, 2 = circle, 3 = diamond, 4 = X, 5 = single arrow, 6 = double arrow, 7 = single triangle, 8 = double triangle ; @param plabel label number: 0 = blank, >0 = number, <0 = number/letter ; @param picon icon (use 0 for blank) ; @param pz location of handle ; @param pr rotation of handle in degrees func SetHandlePoint(int phandle, int ptype, int plabel, complex picon, complex pz, float pr) float logbase = 10.0 ; work-around for log(10.0) optimization introducing round-off errors if (phandle == -1) phandle = m_HandleNext m_HandleNext = m_HandleNext + 1 endif ; save point metadata m_HandleTypes[phandle] = ptype m_HandleLabels[phandle] = plabel m_HandleIcons[phandle] = picon m_HandleRotations[phandle] = (0,1)^(pr/90.0) ; convert point to pixel coordinates now m_PixelCoordinate.Init(pz) m_HandleCenters[phandle] = m_PixelCoordinate.Iterate(pz) ; determine label length now if (plabel > 0) m_HandleLabelSizes[phandle] = floor(log(plabel)/log(logbase))+1 elseif (plabel < 0) logbase = 40.0 m_HandleLabelSizes[phandle] = floor(log(-plabel)/log(logbase))+1 else m_HandleLabelSizes[phandle] = 0 endif endfunc ; Define a handle line segment ;

; Handle lines can be used to complement handle points and ; show relationships between them. They can also be used by ; themselves (there is no requirement that they be used to ; connect handle points). ;

; Lines are automatically hidden in the area around handle ; points and labels/icons. This is so that the lines can ; never obscure the more important parts of the handle ; system. ;

; Lines may be drawn in one of four styles. Solid lines ; should be used only for shapes or elements that are ; otherwise completely invisible. Dotted lines should be ; used to indicate constraining or boundary lines. ; Short dashed lines should be used for a shape boundary ; that is drawn some distance from a visible edge. ;

; When lines are being drawn end-to-end without handles ; at the connecting points, you may wish to indicate that ; the lines are connected. This allows the handle system ; to automatically adjust the pattern on the lines to be ; continuous. If handles are rendered at the connecting ; points they will obscure the join, making this setting ; irrelevant. ; ; @param phandle slot number, starting at 0; if you indicate -1, the next slot will be used automatically ; @param pconnect whether segment connects to previous line segment ; @param ptype line type: 0 = blank, 1 = solid, 2 = dotted, 3 = short dashed, 4 = long dashed ; @param pstart starting point of line segment ; @param pend ending point of line segment func SetHandleLine(int phandle, int ptype, bool pconnect, complex pstart, complex pend) if (phandle == -1) phandle = m_LineNext m_LineNext = m_LineNext + 1 endif ; save line metadata m_LineTypes[phandle] = ptype m_LineConnects[phandle] = pconnect ; convert points to pixel coordinates now ; we include a bit of fudge so that lines don't precisely ; fall on pixel boundaries; this guarantees that lines ; will not be an extra pixel thick due to each side being ; exactly the right number of pixels from the center m_PixelCoordinate.Init(pstart) m_LineStarts[phandle] = m_PixelCoordinate.Iterate(pstart) - (0.001, 0.001) m_PixelCoordinate.Init(pend) m_LineEnds[phandle] = m_PixelCoordinate.Iterate(pend) - (0.001, 0.001) ; compute line length in pixel coordinates now m_LineDistances[phandle] = cabs(m_LineEnds[phandle] - m_LineStarts[phandle]) endfunc ; Define a handle circle ;

; Circles are similar to lines in that they have different ; styles. In the current implementation, circles will always ; appear circular regardless of the image aspect ratio. ; ; @param phandle slot number, starting at 0; if you indicate -1, the next slot will be used automatically ; @param ptype circle type: 0 = blank, 1 = solid, 2 = dotted, 3 = short dashed, 4 = long dashed ; @param pcenter center point of the circle ; @param pradius radius of the circle func SetHandleCircle(int phandle, int ptype, complex pcenter, float pradius) if (phandle == -1) phandle = m_CircleNext m_CircleNext = m_CircleNext + 1 endif ; save circle metadata m_CircleTypes[phandle] = ptype ; convert points to pixel coordinates now m_PixelCoordinate.Init(pcenter) m_CircleCenters[phandle] = m_PixelCoordinate.Iterate(pcenter) m_PixelCoordinate.Init(pcenter+pradius) m_CircleRadii[phandle] = cabs(m_CircleCenters[phandle]-m_PixelCoordinate.Iterate(pcenter+pradius)) endfunc ; use this function to set the screen-relative coordinate style ; if you don't do this, and you allow screen-relative coordinates, ; your handles will appear in the wrong place func SetScreenRelative(int ptype) m_PixelCoordinate.SetScreenRelative(ptype) endfunc protected: PixelCoordinate m_PixelCoordinate complex m_DigitsDetail[40] int m_HandleNext int m_HandleTypes[] int m_HandleLabels[] int m_HandleLabelSizes[] complex m_HandleIcons[] complex m_HandleCenters[] complex m_HandleRotations[] int m_LineNext int m_LineTypes[] bool m_LineConnects[] float m_LineDistances[] complex m_LineStarts[] complex m_LineEnds[] int m_CircleNext int m_CircleTypes[] complex m_CircleCenters[] float m_CircleRadii[] default: title = "Handles utility class" int param v_handles caption = "Version (Handles)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_handles < 100 endparam bool param p_showhandles caption = "Show Handles" default = true hint = "If enabled, shows 'handles' in the image to mark where the control points of the shape are located. This is very useful while you are creating your image, but you will probably not want handles in your final image render. By default, handles should not render in your final image if the render is large, but if you do a small final render you may need to use this master switch to disable handles." endparam bool param p_handlesnopreview caption = "Not on previews" default = true visible = (@p_showhandles) hint = "If checked, handles will not be shown on previews. Handles don't scale with window size, so they appear disproportionately large in the preview window." endparam bool param p_handlesnorender caption = "Not on renders" default = true visible = (@p_showhandles) hint = "If checked, handles will be disabled for any image above a certain size. This is useful for automatically disabling handles on disk renders." endparam bool param p_handleswithnumbers caption = "Show Numbers" default = true visible = (@p_showhandles) hint = "If checked, handles will be numbered." endparam bool param p_handleswithlines caption = "Show Lines" default = true visible = (@p_showhandles) hint = "If checked, lines may also appear." endparam int param p_handlesize caption = "Handle Size" default = 8 visible = (@p_showhandles) hint = "Sets the size of the inside area of the handle. You may need to adjust this to make the handles more visible on some fractals." endparam int param p_handlewidth caption = "Handle Thickness" default = 2 visible = (@p_showhandles) hint = "Sets the thickness of the handles. You may need to adjust this to make the handles more visible on some fractals." endparam } ; ----------------------------------------------------------------------------------------------------------------------------- ; Formula classes ; ----------------------------------------------------------------------------------------------------------------------------- ; Coloring classes ; ----------------------------------------------------------------------------------------------------------------------------- ; Trap Shape classes class TrapShapeBlock(TrapShape) { ; This is a generic trap "block". This is essentially a small ; wrapper around a transformation, a trap shape, and a transfer ; function. You can use it anywhere a TrapShape is needed. public: func TrapShapeBlock(Generic pparent) TrapShape.TrapShape(pparent) m_TrapTransform = new @f_traptransform(this) m_TrapShape = new @f_trapshape(this) m_TrapTransfer = new @f_traptransfer(this) endfunc func Init(complex pz) TrapShape.Init(pz) m_TrapTransform.Init(pz) m_TrapShape.Init(pz) m_TrapTransfer.Init(pz) endfunc float func Iterate(complex pz) complex zt = m_TrapTransform.Iterate(pz) float distance = m_TrapShape.Iterate(zt) distance = m_TrapTransfer.Iterate(distance) return distance endfunc func SetThreshold(float pthreshold) m_TrapShape.SetThreshold(pthreshold) endfunc complex func GetTransformedPoint() return m_TrapShape.GetTransformedPoint() endfunc float func GetTextureValue() return m_TrapShape.GetTextureValue() endfunc int func GetColorChannel() return m_TrapShape.GetColorChannel() endfunc protected: UserTransform m_TrapTransform TrapShape m_TrapShape Transfer m_TrapTransfer default: title = "TrapShape Block" int param v_trapshapeblock caption = "Version (TrapShapeBlock)" default = 101 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapshapeblock < 100 endparam UserTransform param f_traptransform caption = "Trap Position" default = TrapTransform expanded = false endparam TrapShape param f_trapshape caption = "Trap Shape" default = TrapShapePoint endparam Transfer param f_traptransfer caption = "Trap Transfer" default = NullTransfer endparam } class TrapShapePoint(TrapShape) { ; Simple point trap shape. public: func TrapShapePoint(Generic pparent) TrapShape.TrapShape(pparent) endfunc float func Iterate(complex pz) TrapShape.Iterate(pz) return cabs(pz) endfunc default: title = "Point" int param v_trapshapepoint caption = "Version (TrapShapePoint)" default = 100 hint = "This version parameter is used to detect when a change has been made to the formula that is incompatible with the previous version. When that happens, this field will reflect the old version number to alert you to the fact that an alternate rendering is being used." visible = @v_trapshapepoint < 100 endparam } class TrapShapeNoise(TrapShape) { ; Simple noise trap shape. ;

; Most trap shapes have a clearly-defined shape. However, in some ; cases where texturing is desired, a noise basis function may be ; required. This class provides a simple, modified Perlin noise ; basis function which can be used as a stand-alone trap shape or ; as the foundation on which to build more sophisticated textures. ;

; If you are writing a function that needs a noise basis function ; you should use this class (or any of its derivatives) to provide ; it. Then, users may select any other TrapShapeNoise-derived ; class. If you use TrapShape as your parameter, the user may be ; confused by the appearance of other trap shapes which will not ; work well as noise basis functions. ;

; If you want to write your own noise basis function, you should ; derive from this class so that your noise basis can be selected ; to replace any other TrapShapeNoise objects. These should give ; a result over the entire complex plane ranging from 0 to 1 (but ; occasional deviations outside this range should be accepted). ;