%Code by Nathan Kutz, Kyle Mandli, Cleve Moler (polyinterp), and Eric
%Shea-Brown

clear all; close all;

% Load the data
load linefit.dat
x=linefit(:,1);
y=linefit(:,2);
figure(1)
set(gca,'FontSize',20)
plot(x,y,'O')
legend('Raw Data')

% These are the points we would like to interpolate to
xp=0:0.1:7;

% Polynomial fits -- LEAST SQUARES
clin=polyfit(x,y,1);
yp=polyval(clin,xp);
clin2=polyfit(x,y,2);
yp2=polyval(clin2,xp);

figure(2),
set(gca,'FontSize',20)
plot(x,y,'O', ...
    xp,yp,'m',xp,yp2,'g')
legend('Raw Data','Least Squares Line','Least Squares 2nd order polynomial')

% Polynomial fits -- EXACT INTERPOLATION FIT, ORDER N-1 POLYNOMIAL
% FOR N DATAPOINTS, VIA LAGRANGE POLYNOMIALS
ypn=polyinterp(x,y,xp);
figure(3),
set(gca,'FontSize',20)
plot(x,y,'O',xp,ypn,'m')
legend('Raw Data','Poly interpolation via Lagrange')

% Error
E_2_1 = sqrt( 1/length(x) * sum( abs(y-polyval(clin,x)).^2) )
E_2_2 = sqrt( 1/length(x) * sum( abs(y-polyval(clin2,x)).^2) )


% Piecewise interpolation
yint=interp1(x,y,xp); % Linear functions
yspline=spline(x,y,xp); % Cubic functions (splines)
figure(4),
set(gca,'FontSize',20)
plot(x,y,'O',...
    xp,yint,'m',xp,yspline,'k')
legend('Raw Data','Linear Interpolation','Spline interpolation')

load gaussfit.dat
x2=gaussfit(:,1) ;
y2=gaussfit(:,2) ;
figure(5),
set(gca,'FontSize',20)
plot(x2,y2,'O')
legend('Raw Data')

xga=-3:0.1:3;
coeff=fminsearch('gafit',[1 1]);
a=coeff(1); b=coeff(2) ;
yga=a*exp(-b*xga.^2);
figure(6), 
set(gca,'FontSize',20)
plot(x2,y2,'O',xga,yga,'m')
legend('Raw Data','Gaussian nonlinear least-squares fit')