征解:
着色代码:
dmj-Lyapunov {
;
; This algorithm computes the Lyapunov exponent
; for Mandelbrot types. This exponent is usually
; negative for divergent orbits (outside points),
; and positive for convergent orbits (inside
; points).
;
; This might have interesting results when used
; with other fractal types, although the results
; would not be mathematically accurate.
;
; Optimizations suggested by Charles Vassallo.
;
init:
float oldsum = 0
float sum = 1
float v = 0
float il = 1/log(real(@power))
float lp = log(log(@bailout)/2.0)
float f = 0.0
loop:
IF (@trackvariable == 0) ; |z|
v = cabs(#z)
ELSEIF (@trackvariable == 1) ; real(z)
v = real(#z)
ELSEIF (@trackvariable == 2) ; imag(z)
v = imag(#z)
ELSEIF (@trackvariable == 3) ; real(z)/imag(z)
v = real(#z)/imag(#z)
ELSEIF (@trackvariable == 4) ; imag(z)/real(z)
v = imag(#z)/real(#z)
ELSEIF (@trackvariable == 5) ; arg(z)
v = atan2(#z)
ELSEIF (@trackvariable == 6) ; 1/real(z)
v = 1.0/real(#z)
ELSEIF (@trackvariable == 7) ; 1/imag(z)
v = 1.0/imag(#z)
ENDIF
oldsum = sum
; sum = sum + log(abs(2*v)) ; sum the Lyapunov exponent (slow method)
sum = sum * (abs(2*v)) ; sum the Lyapunov exponent
final:
oldsum = log(oldsum)
sum = log(sum)
IF (@negative == 1)
sum = -sum/#numiter
oldsum = -oldsum/(#numiter-1)
ELSEIF (@negative == 2)
sum = abs(sum/#numiter)
oldsum = abs(oldsum/(#numiter-1))
ELSE
sum = sum/#numiter
oldsum = oldsum/(#numiter-1)
ENDIF
IF (@smooth)
f = il*lp - il*log(log(cabs(#z)))
#index = oldsum + (sum-oldsum) * (f+1)
ELSE
#index = sum
ENDIF
default:
title = "Lyapunov"
helpfile = "dmj-pub\dmj-pub-uf-lyapunov.htm"
param trackvariable
caption = "Variable to Track"
default = 0
enum = "magnitude of z" "real part of z" "imaginary part of z" \
"real / imag" "imag / real" "angle of z" "1 / real(z)" "1 / imag(z)"
hint = "Indicates which variable to measure the Lyapunov exponent for."
endparam
param negative
caption = "Sign"
default = 2
enum = "positive" "negative" "absolute value"
hint = "Affects the sign of the exponent. 'Negative' and 'absolute \
value' are useful for inside coloring."
endparam
param power
caption = "Exponent"
default = 2.0
hint = "This should be set to match the exponent of the \
formula you are using. For Mandelbrot, this is 2."
endparam
param bailout
caption = "Bailout"
default = 1e20
min = 1
hint = "This should be set to match the bailout value in \
the Formula tab. Use a very high bailout!"
endparam
param smooth
caption = "Smooth Coloring"
default = false
hint = "If set, results will be 'smoothed' to hide iteration bands."
endparam
} |