%% Polynomial Data Fitting & Trigonometric Functions 
%
% Given a set of points sampled form a trig function in R^2, find the polynomial that contains them.  
 
%% Define the trig function you want to integrate
	Given_f = @(x) sin(x);
	%Given_f = @(x) exp(x.^2);
	%Given_f = @(x) sin(pi*x.^2)/2;
	
%% We choose random points on the graph of $$\sin{\pi x^2 / 2}$.
	k = 10; % number of points
    x_given = -2+4*rand(1,k);
	y_given = Given_f(x_given);

	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 = max(abs(x_given))+6; % 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
	disp('The approximating polynomial is: ')
	disp(s)
	PolyApprox = @(x) eval(s);
	
%% Plot the points and the curve in a psudo-animation

	close all; % close all other figures  
	figure; % create a new figure
	ry = max(abs(Given_f(x_given)))+20;
	ylim([-ry ry]); % 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('A trig function and its approximation'); % puts the string s into the title

	% Plot the given points as large blue stars,
	plot(x_given,y_given,'*b','MarkerSize',12);
	disp('The given points are in blue')
	legend('Given Points') % create a legend
	pause

	xx = [-r:.001:r]; % sample x-values, the domain of the function
	yy = Given_f(xx);
	plot(xx,yy,'sk','MarkerSize',1) % plot the exact function in black
	disp('The given function is in black')
	legend('Given Points', 'Given Function')
	pause % pause the printing for greater appreciation

	% Plot the graph of the polynomial along with the given points
	plot(x,PolyApprox(x),'-r','MarkerSize',2);
	plot(x_given,Given_f(x_given),'*b','MarkerSize',12);
	disp('The polynomial approximation is in red')
	legend('Given Points','Given Function', 'Polynomial Approx')

%%  One could evaluate the exact and approximate integrals by calling the script CompareIntegrals afterward.
