#hw11.txt.  February 28, 2014.  Katie McKeon.

read(`ExpMath/C11.txt`):

with(plots):

Help11:=proc()
print(`NumerBifOrbF(f,x,x0,K1,K2,d1,n,eps) FirstBif(f,x) PlotMS(reso,K)`):
end:


#NumerBitOrbF inputs a function f(x), initial point x0, integer list sizes
#K1 and K2, digit place d1, and outputs the intervals surrounding n bifurcation
#points g of x -> gf(x) with interval size at most eps
NumerBifOrbF:=proc(f,x,x0,K1,K2,d1,n,eps) local g, B, i, orbSize,period:

g:=2.9:
B:=[]:
period:=1:

while nops(B)<n do
	g:=g+eps:
	orbSize:=nops({op(OrbF(f,x,g,x0,K1,K2,d1))}):

	if orbSize>2*period then 
		period:=2*period:
		B:=[op(B), [g-eps/2,g+eps/2]]:
	elif orbSize = 2*period then
		period:=2*period:
		B:=[op(B), [g-eps, g]]:	
		
	fi:
od:

evalf(B,5):

end:



#DATA
#NumerBifOrbF returns the following intervals for
#x^2*(1-x)   # [[5.32, 5.34], [5.76, 5.78], [5.86, 5.88], [5.8900, 5.9100]]
#x*(1-x)^2   # [[4.00, 4.02], [4.9900, 5.0100], [5.22, 5.24], [5.2900, 5.3100]]
#x^2*(1-x)^2 # [[10.40, 10.42], [11.810, 11.830], [12.14, 12.16], [12.22, 12.24]]



#FirstBif(f,x) inputs a function f(x) and determines the first bifurcation point
#of the generalized logistic map x-> g*f(x)
FirstBif:=proc(f,x) local h,g,Sol,i:

h:= unapply(g*f, x):
Sol:=[solve(discrim(normal(((h@h)(x)-x)/(h(x)-x)), x), g)] :
for i from 1 to nops(Sol) do
	if Sol[i] > 0 then
		RETURN(Sol[i]):
	fi:
od:

end:



#DATA
#FirstBif returns the following first bifurcation points
# x^2*(1-x)	# 16/3
# x*(1-x)^2	# 4
# x^2*(1-x)^2	# 125/12
# x^3*(1-x)^3	# -returns nothing because the discriminant appears to only
#		   have 0 as a root




#PlotMS(reso,K) inputs a resolution reso and positive integer K and 
#plots an approximation of the Mandelbrot set by iterating z-z^2+c K times
#for every point c in the disc of radius 2 with grid size reso
PlotMS:=proc(reso,K) local i, j, k, z, c, S:

S:=[]:
for i from 0 to 4/reso do
	for j from 0 to 4/reso do
		c:=Complex(-2+i*reso,-2+j*reso):
		if abs(c)<=2 then
			z:=0:
			for k from 1 to K do
				z:=z^2+c:
			od:
			if abs(z)<=2 then
				S:=[op(S),c]:
			fi:
		fi:
	od:
od:

complexplot(S,x=-2..2,y=-2..2,style=point):

end: