//sfb 3 dec 21 adding comments

clear simulate
function [len1, len2, len3] =simulate(style)
	//r = 1;
	theta=0:0.01:2*%pi;
	scf(0); close();
	drawlater;
	f=scf(0);
	f.figure_position= [0,0];
	clf;
	drawnow;
	a=gca();
	a.isoview="on"; //sleep(1000)
	drawlater;
	plot([-1.2,6.2],[-1.4,-1.4],'w');
	plot(cos(theta),sin(theta),'k-','linewidth',3);      //circle 1
	xstring(-.5,-1.3,"Random Radii");
	plot(2.5+cos(theta),sin(theta),'k-','linewidth',3);  //circle 2
	xstring(1.8,-1.3,"Random Darts");
	plot(5+cos(theta),sin(theta),'k-','linewidth',3);    //circle 3
	xstring(4.3,-1.3,"Random Endpoints");
	a=gca();
	a.axes_visible=["off","off","off"];
	a.box="off";
	drawnow;
	n = 200;  //sample size
	n = 400;
	// data: len = length of chord, a and b endpoints
	len1=zeros(n,1); len2=zeros(n,1); len3=zeros(n,1);
	x1a=zeros(n,1); x2a=zeros(n,1); x3a=zeros(n,1);
	y1a=zeros(n,1); y2a=zeros(n,1); y3a=zeros(n,1);
	x1b=zeros(n,1); x2b=zeros(n,1); x3b=zeros(n,1);
	y1b=zeros(n,1); y2b=zeros(n,1); y3b=zeros(n,1);
    // random radii
	for i=1:n,
		t = rand()*%pi;   // between 0 <= t <= pi  -- direction
		h = rand()*2-1;   // between -1 <= h <= 1  -- distance from center
		len = 2*sqrt(1-h^2);  // sides h, 1 and sqrt(1-h^2) -- 2x
		len1(i) = len;
		if style then
			t=0;       // normalize with i-hat the random directions
		end;
		x0 = h*cos(t); y0 = h*sin(t);  // coordinates of point
		// chord coordinates
		x1a(i)= x0+len/2*sin(t); y1a(i)= y0-len/2*cos(t);
		x1b(i)= x0-len/2*sin(t); y1b(i)= y0+len/2*cos(t);
	end
    // random endpoints
	for i=1:n,
		t1 = rand()*2*%pi;   // point1
		t2 = rand()*2*%pi;   // point2
		if style then
			t2 = t2 - t1;
			t1 = 0;          // normalize with point1 at (1,0)
		end
        // what is h????
		// x0 = h*cos(t); y0 = h*sin(t);   
		// not used???
		x3a(i)= cos(t1); y3a(i) = sin(t1);
		x3b(i)= cos(t2); y3b(i) = sin(t2);
		len = sqrt( (x3a(i) - x3b(i)) .^2 + (y3a(i) - y3b(i)) .^2 );
		len3(i) = len;
	end
    // random darts
	i = 1;
	while  i <= n,   // because we throw away points outside the circl3
		x0 = rand()*2-1;    // -1 <= x0 <= 1
		y0 = rand()*2-1;    // -1 <= y0 <= 1
		r = sqrt(x0^2+y0^2);
		t = atan(y0,x0);      // polar coordinate theta 
		if style then
			x0 = sign(y0)*r; y0 = 0;
			t=0;     // normalize to lie on x-axis
		end
		// if (x0, y0) is inside the circle
		if r>0 & r < 1 then  
			len = 2*sqrt(1-r^2);  // rt triangle r, 1, sqrt(1-r^2)
			x2a(i)= x0+len/2*sin(t); y2a(i)= y0-len/2*cos(t);
			x2b(i)= x0-len/2*sin(t); y2b(i)= y0+len/2*cos(t);
			len2(i) = len;
			i = i + 1;
		end
	end
    // draw stuff
	for i=1:n,
		drawlater();
		//chord
		xpoly([x1a(i), x1b(i)],[y1a(i),y1b(i)]);
		e=gce(); set(e,"foreground",5);
		//chord
		xpoly(2.5+[x2a(i), x2b(i)],[y2a(i),y2b(i)]);
		e=gce(); set(e,"foreground",2);
		//chord
		xpoly(5+[x3a(i), x3b(i)],[y3a(i),y3b(i)]);
		e=gce(); set(e,"foreground",26);
		//a.axes_visible=["off","off","off"];
		drawnow();
		sleep(100);
	end
endfunction

[res1,res2,res3] = simulate(1);
	scf(1);close();f=scf(1);
	f.figure_position= [610,0];
	f.figure_size=[1320,648];
	clf;
	show_window();
	classes=0:0.1:2;
	subplot(1,3,1);
	histplot(classes,res1);
	med = median(res1);
	avg = sum(res1)/size(res1,1);
	title("$\mbox{Median Length = } \sqrt{3} = 1.73205$",'fontsize',3);
	xstring(0.1,2.0,"Random Radii");
	e=gce(); e.font_size=3;
	xstring(0.1,1.8,"Sample Average = " + string(avg));
	e=gce(); e.font_size=3;
	xstring(0.1,1.6,"Sample Median = " + string(med));1
	e=gce(); e.font_size=3;

	subplot(1,3,2);
	histplot(classes,res2);
	med = median(res2);
	avg = sum(res2)/size(res2,1);
	title("$\mbox{Median Length = } \sqrt{2} = 1.41421$",'fontsize',3);
	xstring(0.1,0.8,"Random Darts");
	e=gce(); e.font_size=3;
	xstring(0.1,0.7,"Sample Average = " + string(avg));
	e=gce(); e.font_size=3;
	xstring(0.1,0.6,"Sample Median = " + string(med));
	e=gce(); e.font_size=3;

	subplot(1,3,3);
	histplot(classes,res3);
	med = median(res3);
	avg = sum(res3)/size(res3,1);
	title("$\mbox{Median Length = } \sqrt{2} = 1.41421$",'fontsize',3);
	xstring(0.1,1.2,"Random Endpoints");
	e=gce(); e.font_size=3;
	xstring(0.1,1,"Sample Average = " + string(avg));
	e=gce(); e.font_size=3;
	xstring(0.1,0.8,"Sample Median = " + string(med));
	e=gce(); e.font_size=3;

