lines(0);

// do not change this function
function dwdt = lvm(t,w)
// w(1) contains C, w(2) contains R
// global alpha Beta lambda mu;
dwdt=[alpha*w(1,1) - lambda*w(1,1)*w(2,1); - Beta*w(2,1) + mu*w(1,1)*w(2,1)];
endfunction;

// old stuff
function z = datazero(w)
// usage datazero(w)
// assumes domain is the column one w(:,1) and the range is column two w(:,2)
// solves for domain values x where the range is 0
	x = w(:,1); y = w(:,2); ay = abs(y);
	jy = ay(1:$-1)+ay(2:$);
	oscilate = y(1:$-1) .* y(2:$);
	i = find(oscilate < 0);
	g0 = (ay(i+1).*x(i)+ay(i).*x(i+1))./jy(i);
	d0 = x(find(y==0));
	z = gsort([d0; g0],'g','i');
endfunction;

// From datasolve.m
function x0 = datasolve(y0, w)
// usage datasolve(y0, w)
// assumes domain is column one w(:,1) range is column two w(:,2) 
// domain is often the x value, range is often the y value.
// evaluates solves for x when the range y = x0
	w(:,2) = w(:,2) - y0;
	x0 = datazero(w);
endfunction;

// From dataeval.m
function y0 = dataeval(x0, w)
// usage dataeval(x0, w)
// assumes domain is column one w(:,1) and the range is column two w(:,2)
// evaluates data at x = x0
	x = w(:,1); y = w(:,2); xx = x - x0;
	eval0 = y(find(xx==0));
	oscilate = xx(1:$-1) .* xx(2:$);
	i = find(oscilate < 0);
	eval2 = (xx(i+1).*y(i)-xx(i).*y(i+1))./(xx(i+1)-xx(i));
	y0 = (gsort([eval0 eval2],'g','i'))';
endfunction;

function datausage()
	printf("\nGiven x0 want y0:\n");
	printf("Use dataeval(x0, [x-data-column, y-data-column])\n");
	printf("\nGiven y0 want x0:\n");
	printf("Use datasolve(y0, [x-data-column, y-data-column])\n");
endfunction

// matrix m is an array of frames
// n is the number of times to play the movie
// if n < 0, go forward and backward |n| times
// fps is frames per second default 12 fps
// delay of 40 milliseconds gives 12 fps on one computer
// pixmap command means the graphs are drawn on a hidden
// campus and switch to the front with the show_pixmap 
// command
function movie(m,n,delay)
	[x y d] = size(m);
	f = gcf();
	f.pixmap="on";
	square(0,0,x,y);
	for k = 1:abs(n),
		for i = 1:d,
			clf;
			Matplot(m(:,:,i));
			if (delay > 0), sleep(delay); end;
			show_pixmap;
		end
		if n < 0,
			for i = d:-1:1,
				clf;
				Matplot(m(:,:,i));
				if (delay > 0), sleep(delay); end;
				show_pixmap;
			end
		end
	end
	f.pixmap="off";
endfunction;

MaxIt=60;
n = 256;
t = 2*n + 1;

function A = manarray2(x0, y0, scale)
MaxIt = 60;
n = 256;
xmin = x0 - scale; xmax = x0 + scale;
ymin = y0 - scale; ymax = y0 + scale;
//printf("xmin = %f; xmax=%f; ymin = %f; ymax = %f\n", xmin, xmax, ymin, ymax);
//[X, Y] = meshgrid(xmin:scale/n:xmax, ymin:scale/n:ymax);
[X, Y] = meshgrid(xmin:scale/n:xmax, ymax:-scale/n:ymin);
Z = X + %i*Y;
A = MaxIt*ones(Z);
value = zeros(Z);
//printf("Computing .");
for i = 1:MaxIt;
	value = value .* value + Z;	
	current_mask = (abs(value) > 2);
	B = i * current_mask + MaxIt * (~current_mask);
	A = min(A,B);
	//value = value .* (~current_mask);
	//if modulo(i,5)  == 0, printf("."); end;
end;
endfunction;

function setcolormap()
MaxIt=60;
cm = ones(MaxIt+1,3); // everything white
cm(MaxIt,:) = [ 0 0 0 ]; // black
cm(:,1:2) = zeros(MaxIt+1,2); //blue
cm(1,:) = [ 1 1 1 ]; //white
cm(2,:) = [ 1 1 1 ]; //white
cm(3,:) = [ 1 1 1 ]; //white
cm(4,:) = [ 0 1 0 ]; //green
cm(5,:) = [ 0 1 0 ]; //green
for i = 6:10, cm(i,:) = [0 0.5 0.5]; end; 
for i = 11:20, cm(i,:) = [0 0.75 0.75]; end;
for i = 21:40, cm(i,:) = [0 1 1]; end; //cyan
for i = 40:(MaxIt-1), cm(i,:) = [0 0 1]; end; //blue
cm(MaxIt+1,:) = [ 1 0 0 ]; // red
scf(0); clf; square(0,0,t,t);
f=gcf();f.color_map = cm;
endfunction;


function m=mandelbrot()
MaxIt=60;
n = 256
t = 2*n + 1;
// Create the color map
setcolormap();
//
scale = 2;
stacksize('max');
m=zeros(t,t,10);
//x0 = -1/2;
x0 = 0;
y0 = 0;
for i = 1:10, 
    a = manarray2(x0, y0, scale);
	//a(n+1,:) = MaxIt+1;
	//a(:,n+1) = MaxIt+1;
	drawlater;
	clf; Matplot(a); square(0,0,t,t);
    s=gca(); s.tight_limits="on"; s.auto_ticks=["off","off","off"];
	//s.axes_reverse=["off","on","off"];
	s.axes_reverse=["off","off","off"];
	drawnow;
    m(:,:,i) = a;
	if i < 10,
	disp('Click at Zoom-in Location'); here = locate(1) // announce
	//printf("%f, %f\n", here(1), here(2));
    column = int(0.5+here(1)); row = int(0.5 + here(2));
    x0 = x0-scale + (column-1)*scale/n;
    y0 = y0-scale + (row-1)*scale/n;
    scale = scale/2;
	end;
	//printf("Center at %f, %f\n", x0, y0);
end;
endfunction;
