%
% fractal1d.t
%
% Recursive algorithm to generate "smoothed" 1-D fractional Brownian motion
%
const	maxsize := 1024
const	norands := 5				% no of random variables

type output: array 0..maxsize of real 
var Fh: output					% global var for efficiency


function gauss(s: real) : real
	var r : real
	var sum : real := 0

	for i: 1..norands
		rand(r)
		sum := sum + (r-0.5)*2
	end for
	r := s*sum/norands
	result r
end gauss


proc subdivide( f1,f2: int, std, ratio: real)
	var fmid:  int
	var stdmid: real

	fmid := (f1 + f2) div 2
	if fmid not= f1 and fmid not= f2 then
		Fh(fmid) := (Fh(f1)+Fh(f2))/2.0 + gauss(std);
		stdmid := std*ratio
		subdivide(f1,fmid,stdmid,ratio)
		subdivide(fmid,f2,stdmid,ratio)
	end if
end subdivide

	
		
procedure fractal(maxlevel: int, h, scale: real)
	var first, last: int
	var ratio, std:  real

	first := 0
	last  := 2**maxlevel
	Fh(first) := gauss(scale)
	Fh(last)  := gauss(scale)
	ratio     := 2**(-h)
	std       := scale*ratio
	subdivide(first,last, std,ratio)
end fractal



var scale, h: real

const level := 10
setscreen("graphics:1024;800")
put "Fractal demonstration."
put "Enter smoothness value [0,1]: " ..
get h
setscreen("nocursor")
fractal(level,h,256)
randomize
for i : 0..2**level-2
	drawline(i*1,round(400+Fh(i)), (i+1)*1,round(400+Fh(i+1)),1);	
end for