%% Polynomial Data Fitting 
% Given a set of points sampled form a trig function in R^2, find the polynomial that contains them.  
% One could evaluate the exact and approximate integrals by calling the script CompareIntegrals afterward.

% Possibly create an alternative to polyfit, and ask students to model climate data.


%% We choose random points on the graph of $$\sin{\pi x^2 / 2}$$.
	%x_given = -2+4*rand(1,20);
	%y_given = sin ((pi*x_given.^2 )/2);

	x_given = [0:2:24];
	y_given = 1/100*[34 50 60 64 63 59 53 46 40 36 34 37 45];

	n = length(x_given);

%% The "Coefficient" Matrix
% The polynomial is a_0 + a_1*x + ... + a_m*x^m where m = n-1;
% This general form determines one equation for every given pair
	M = ones(n,n+1);
	for point = 1:n
		for i = 1:n;
			M(point,i) = x_given(point)^(i-1);
		end
		M(point,n+1) = y_given(point);	
	end

%% Solve the system
	A = rref(M);
	a = A(:,end);

%% The polynomial
% For example, y = a(1) + a(2)*x + a(3)*x.^2;
	r = 70; % How wide is the range of the plotting window
	x = -r:.0001:r; % a sample of x values in the chosen range 
	s = [num2str(a(1))]; % make a string out of the number held in a(1)

	% term by term, make a polynomial whose coeffs are taken from a
	for i = 2:length(a); 
	 	% s is a character string that looks like the poly we want
		s = [s,' + ',num2str(a(i)),'*x.^',num2str(i-1)];
	end
	% now we turn the character string into a set of y-values of the polynomial
	y = eval(s);
	
%% Plot the points and the curve in a psudo-animation

	close all; % close all other figures  
	figure; % create a new figure
	ylim([-r r]); % the max/min of the y-axis
	xlim([-r r]); % the max/min of the x-axis
	hold on % put all the plots onto the same figure
	title(s); % puts the string s into the title

	%xx = [-r:.001:r]; % sample x-values, the domain of the function
	%yy = sin((pi*xx.^2)/2); % the exact function
	%plot(xx,yy,'sk','MarkerSize',1) % plot the exact function in black
	%legend('Exact') % create a legend
	pause % pause the printing for greater appreciation

	% Plot the given points as large blue stars,
	plot(x_given,y_given,'*b','MarkerSize',12);
	legend('Exact', 'Interpolating Points')
	pause

	% Plot the graph of the polynomial along with the given points
	plot(x,y,'or','MarkerSize',2);
	plot(x_given,y_given,'*b','MarkerSize',12);
	%legend('Exact','Interpolating Points', 'Polynomial Approx')
