#Header Begins--Do not remove##################################################
#
#	Maple code to create 1-D, deterministic fractal curves starting from
#	initial sequence of line segments.
#
#	(c) Eugene Fiume, University of Toronto, 1995.
#
#	You are permitted to use this code for noncommercial, educational
#	purposes only.  This header must not be removed from this file.
#
#Header Ends--Do not remove####################################################

#
#	Replace the sequence of line segments in "current" using the
#	replacement rule called "rule".  This rule is the name of
#	a Maple procedure that is repeatedly called to perform the
#	replacement.
#

Replace := proc(current,rule)
	local n,i,new;

	n := nops(current);
	new := NULL;
	for i from 1 to n-1 do
		new := new, rule(current[i],current[i+1]);
	od;
	[new];
end:


#
#	Replacement rule for von Koch curve.
#
RepVonKoch := proc(P,Q)
	local P1, Q1, C;

	P1 := Lirp(P,Q,1/3);
	Q1 := Lirp(P,Q,2/3);
	C  := Rotate(Q1,P1,Pi/3);
	P,P1,C,Q1,Q;
end:


#
#	Rotate P=[px,py] about Q=[qx,qy] through angle ang.
#
Rotate := proc( P, Q, ang)
	local s, c, xx,x, y;

	s := evalhf(sin(ang));
	c := evalhf(cos(ang));
	x := P[1]-Q[1];
	y := P[2]-Q[2];
	xx:=  (c*x - s*y + Q[1]);
	y :=  (s*x + c*y + Q[2]);
	[xx,y];
end:

#
#	Return point that is the fraction t between P and Q.
#	"Lirp" stands for linear interpolation.
#
Lirp := proc(P,Q, t)
	[P[1] + t*(Q[1]-P[1]), P[2] + t*(Q[2]-P[2])]:
end:

#
#	Apply rule n times starting with initial geometry.
#
plotfractal := proc(geom,n,rule)
	local i,g;

	g := geom;
	for i from 1 to n do
		g := Replace(g,rule):
	od:
	print(plot(g,axes=NONE,scaling=CONSTRAINED));
end:

#
#	Here are some initial sequences that are fun to start with.
#
InitVonKoch2 := [[0,0], [1,0], Rotate([1,0],[0,0], evalf(Pi/3)),[0,0]]:
InitVonKoch1 := [[0,0], [1,0], Rotate([1,0],[0,0], evalf(-Pi/3)),[0,0]]:
InitVonKoch  := [[0,0], [1,0], Rotate([2,0],[1,0], evalf(Pi/3)),[2,0], [3,0]]:
Line         := [[0,0], [1,0]]:

#
#	Sample script.
#
#	g := Line:
#	for i from 1 to 3 do
#		g := Replace(g,RepVonKoch);
#	od;
#	plot(g,axes=none,scaling=constrained);
#
#
#	set g to other sequences for variety.
#	put "print(plot(g,axes=none,scaling=constrained));" within the
#	for loop if you want to see intermediate results.
#