Basic-tna { ; extension of basic coloring alrogithm in standard.ucl ; by Ted Nason astrolathe@comcast.net ; ; UNDER DEVELOPMENT SEPTEMBER 2004! USE AT YOUR OWN RISK! ; init: complex f[50000] ; array to hold the points of the orbit complex sum = (0,0) int iter = 0 int finaliter = 0 float rescape = 1000000000; if the stat being calced gets too big loop: iter=iter+1 f[iter]=#z final: finaliter=iter sum = (0,0) if(@type == 0); mod(lastz) #index=|f[finaliter]| elseif(@type == 1); phase(lastz) if imag(f[finaliter]) !=0 #index=atan(real(f[finaliter])/imag(f[finaliter])) endif elseif(@type == 2); zlast - zfirst #index = 10*(real(#z) * imag(#z)) elseif(@type == 3); mod(mean(last 3 pts)) #index = |f[finaliter-2] + f[finaliter-1] + f[finaliter]| elseif(@type == 4); phase(mean(last 3 pts)) sum = f[finaliter-2] + f[finaliter-1] + f[finaliter] if imag(sum) !=0; is this really necessary??? #index = real(sum)/imag(sum) endif elseif(@type == 5); angle(1st vertex) ; #index = xxx elseif(@type == 6); angle(last vertex) ; #index = xxx elseif(@type == 7); eValue(lastz) iter=0 sum=(1,1) \$define DEBUG while iter #maxiter nmaxx = #maxiter ; you can't calculate a fourier period which exceeds the max iterations :) endif while n <= nmaxx; no sense in computing ft elements for periods greater than plot range sum = 0 ; reset on every component k = 0 while k <= iterr sum = sum + f[k]*cos(2*#pi*n*k/iterr) ; print ("k = ", k, " ; sum = ", sum, " ; f[k] = ", f[k]) k = k + 1 endwhile ft[n] = (sum/(iterr+1))^2 ; print ("n = ", n, "; ft[n] = ", ft[n]) n = n + 1; this is the Fourier component of period n in index space ("time" domain) endwhile ; now calculate the peak in frequency domain (1/index, or 1/"time" domain) int mmin = floor(iterr/nmaxx) ; these are the time domain limits corresponding to nmin, nmax int mmax = floor(iterr/@nmin) ; within which to look for a peak in the spectrum int m = 0 ; frequency index int nn = 0 ; period corresponding to freq index m float fmax = 0 ; this is the largest element of ft[] in the range nmin < n < nmax int maxindex = 0 ; value of n at peak if mmin != mmax ; must be some range to plot m = mmin ; find the highest fourier component in range of mmin to mmax while m <= mmax nn = floor(iterr/m) ; calculate period from freq if ft[m] > fmax fmax = ft[m] ; continually find the greatest value of ft maxindex = nn ; and assign the period nn to the maxindex endif ; print ("m = ", m, "; ft[m] = ", ft[m]) m = m + 1 endwhile endif float indexx = 0 if(@stattype==0); only code by index at peak #index = maxindex elseif(@stattype==1) ; sum sqr(peaks) - sqr helps bring low values (<< 1) up m = mmin ; normalize the components in range of mmin to mmax while m <= mmax nn = floor(iterr/m) ; calculate period from freq ft[m] = ft[m]/fmax ; normalize ; print ("ft[n] = ", ft[m]) m = m + 1 endwhile sum = 0 m = mmin while m <= mmax sum = sum + sqrt(ft[m]) ; print ("sum = ", sum, "; ft[n] = ", ft[m]) m = m + 1 endwhile sum = sum - 1; remove influence of major peak (= 1) ; indexx = maxindex + floor(10*(indexx)) ; combine influence of major peak and minor peaks indexx = floor(maxindex*(10*(indexx))) ; print ("maxindex = ", maxindex, "; sum = ", sum, "; ft[n] = ", ft[m]) #index = indexx ; #index = maxindex/40 elseif(@stattype==2); sum sqr(sqr(peaks)) else; sum exp(peaks) endif default: title = "FourierFreq" helpfile = "xxx.html" int param nmin caption = "min Fourier cpt" default = 3 hint = "min Fourier cpt to plot - best if >2, else 2 dominates" endparam int param nmax caption = "max Fourier cpt" default = 50 hint = "max Fourier cpt to plot - MUST be less than max iterations or meaningless" endparam param vartype caption="variable" default=2 enum= "real" "imag" "magnitude" "imag/real" "angle to origin" "angle(k k-1)" "eigenvalue" "eigenphase" hint="variable for which statistics will be calculated" endparam param stattype caption="subordinate sum" default=1 enum="none" "sqr" "^1/4" "exp" hint="method of incorporating minor peaks" endparam } Autocorrelation { ; autocorrelation coloring algorithm ; by Ted Nason astrolathe@comcast.net ; ; UNDER DEVELOPMENT SEPTEMBER 2004! USE AT YOUR OWN RISK! ; init: int iterr = 0 float sum = 0 float f[500] ; array to hold the function of the iterate to analyze loop: iterr = iterr + 1 f[iterr] = |#z| ; put case select here - for now, analyze the FT spectrum of |#z| sum = sum + f[iterr] ; this is used to get the avg final: float avg = sum/iterr float auto[500] ; array to hold the autocorrelation components \$define DEBUG print ("x = ", real(#z), " ; y = ", imag(#z), "avg = ", avg) int k = 1 int n = 0 float sum1 = 0 float sum2 = 0 while n <= iterr-1 sum1 = 0 sum2 = 0 k = 1 while k <= iterr - n sum1 = sum1 + (f[k]-avg)*(f[k+n]-avg) sum2 = sum2 + (f[k]-avg)*(f[k]-avg) k = k + 1 endwhile auto[n] = sum1/sum2 n = n + 1 endwhile int nmaxx = @nmax; now calculate the peak in specified range if @nmax > iterr nmaxx = iterr endif int nn = @nmin float fmax = 0 ; largest element of auto[] in the range nmin < n < nmax int maxindex = 0 ; index at peak while nn < nmaxx if auto[nn] > fmax ; find the highest component in range of nmin to nmax fmax = auto[nn]; print ("nn = ", nn, "; auto[nn] = ", auto[nn]) maxindex = nn endif nn = nn + 1 endwhile print ("maxindex = ", maxindex) #index = maxindex default: title = "Autocorrelation" helpfile = "xxx.html" int param nmin caption = "min Autocor cpt" default = 3 hint = "min cpt to plot - best if >2, else 2 dominates" endparam int param nmax caption = "max Autocor cpt" default = 50 hint = "max cpt to plot - MUST be less than max iterations or meaningless" endparam } PowerSpectrum { ; power spectrum coloring algorithm ; by Ted Nason astrolathe@comcast.net ; ; UNDER DEVELOPMENT SEPTEMBER 2004! USE AT YOUR OWN RISK! ; init: int iterr = 0 float sum = 0 float f[500] ; array to hold the function of the iterate to analyze loop: iterr = iterr + 1 f[iterr] = |#z| ; put case select here - for now, analyze the FT spectrum of |#z| sum = sum + f[iterr] ; this is used to get the avg final: float avg = sum/iterr float auto[500] ; array to hold the autocorrelation components \$define DEBUG print ("x = ", real(#z), " ; y = ", imag(#z), "avg = ", avg) int k = 1 int n = 0 float sum1 = 0 float sum2 = 0 while n <= iterr-1 sum1 = 0 sum2 = 0 k = 1 while k <= iterr - n sum1 = sum1 + (f[k]-avg)*(f[k+n]-avg) sum2 = sum2 + (f[k]-avg)*(f[k]-avg) k = k + 1 endwhile auto[n] = sum1/sum2 n = n + 1 endwhile ; now we have the autocorrelation elements, auto[n] ; range is 0 #maxiter nmaxx = #maxiter ; you can't calculate a fourier period which exceeds the max iterations :) endif sum = 0 n = 0 ; this is period in index space ("time" domain) k = 0 ; index to count through while n <= iterr -1; should only go up to nmaxx? sum = 0 ; reset on every component k = 0 while k <= iterr -1 sum = sum + auto[k]*cos(2*#pi*n*k/iterr) ;print ("k = ", k, " ; sum = ", sum, " ; auto[k] = ", auto[k]) k = k + 1 endwhile f[n] = (sum/(iterr+1))^2; this is the Fourier component of period n in index space ("time" domain) print ("n = ", n, "; ft[n] = ", f[n]) n = n + 1 endwhile ; now calculate the peak in frequency domain (1/index, or 1/"time" domain) int mmin = floor(iterr/nmaxx); these are the time domain limits corresponding to nmin, nmax int mmax = floor(iterr/@nmin); within which to look for a peak in the spectrum int m = 0 ; frequency index int nn = 0 ; period corresponding to freq index m float fmax = 0 ; this is the largest element of ft[] in the range nmin < n < nmax int maxindex = 0 ; value of n at peak ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;if mmin <> mmax ; must be some range to plot m = mmin ; find the highest fourier component in range of mmin to mmax while m <= mmax nn = floor(iterr/m) ; calculate period from freq if f[m] > fmax fmax = f[m] ; continually find the greatest value of ft maxindex = nn ; and assign the period nn to the maxindex endif ; print ("m = ", m, "; ft[m] = ", f[m]) m = m + 1 endwhile print ("maxindex = ", maxindex) #index = maxindex default: title = "PowerSpectrum" helpfile = "xxx.html" int param nmin caption = "min Autocor cpt" default = 3 hint = "min cpt to plot - best if >2, else 2 dominates" endparam int param nmax caption = "max Autocor cpt" default = 50 hint = "max cpt to plot - MUST be less than max iterations or meaningless" endparam } Fourier { ; discrete Fourier transform coloring algorithm ; by Ted Nason astrolathe@comcast.net ; ; UNDER DEVELOPMENT SEPTEMBER 2004! USE AT YOUR OWN RISK! ; init: int iterr = 0 float f[500] ; array to hold the function of the iterate to analyze loop: f[iterr] = |#z| iterr = iterr + 1 ; put case select here - for now, analyze the FT spectrum of |#z| final: \$define DEBUG print ("x = ", real(#z), " ; y = ", imag(#z)) ; at this point #numiter (=iterr) is either = #maxiter (inside) ; or some smaller value where orbit escaped ; would like to dim ft dynamically & smaller if possible... float ft[500] ; array to hold the Fourier components float sum = 0 int n = 0 ; this is period in index space ("time" domain) int k = 0 ; index to count through float nmaxx = @nmax if @nmax > #maxiter nmaxx = #maxiter ; you can't calculate a fourier period which exceeds the max iterations :) endif while n <= iterr ; should only go up to nmaxx... ; (if you only want to find peaks below @nmax, ; no sense in computing ft's for periods greater than that...) sum = 0 ; reset on every component k = 0 while k <= iterr sum = sum + f[k]*cos(2*#pi*n*k/iterr) ; print ("k = ", k, " ; sum = ", sum, " ; f[k] = ", f[k]) k = k + 1 endwhile ft[n] = (sum/(iterr+1))^2 ; print ("n = ", n, "; ft[n] = ", ft[n]) n = n + 1 ; this is the Fourier component of period n in index space ("time" domain) endwhile ; now calculate the peak in period domain (index or "time" domain) int nn = @nmin ; find the peaks float fmax = 0 ; this is the largest element of ft[] in the range nmin < n < nmax int maxindex = 0 ; value of n at peak while nn < nmaxx if ft[nn] > fmax ; find the highest fourier component in range of nmin to nmax fmax = ft[nn] maxindex = nn endif ; print ("nn = ", nn, "; ft[nn] = ", ft[nn]) nn = nn + 1 endwhile #index = maxindex ; print ("maxindex = ", maxindex) default: title = "Fourier" helpfile = "xxx.html" int param nmin caption = "min Fourier cpt" default = 3 hint = "min Fourier cpt to plot - best if >2, else 2 dominates" endparam int param nmax caption = "max Fourier cpt" default = 50 hint = "max Fourier cpt to plot - MUST be less than max iterations or meaningless" endparam }