%       This is a sample matlab program.
%       This file should have name surface.m
%       It runs with the data file surface1

%       To run it from matlab, just type the file
%       name, i.e., surface_fit

%       Prof. Strikwerda April 1995, modified Dec. 1999.

%       If any matlab statements are unclear to you,
%       either use the help command in matlab, e.g.,
%        help load
%       or e-mail to Prof. Strikwerda (strik) or one
%       of the TAs.

%              Clear old variables.
clear
%              Release current plots.
hold off

%              Load the data file.
load surface1

%              How many rows of data?
N = length( surface1 ) 

%              Rename the data columns.
xx = surface1(:,1)  ;
yy = surface1(:,2)  ;
zz = surface1(:,3)  ;

%              Clear unneeded storage.
clear surface1

%              Set up the least squares matrix.
for i = 1:N
%                       constant term
   A(i,1 ) = 1 ;
%                       linear terms
   A(i,2 ) = xx(i) ;
   A(i,3 ) = yy(i) ;
%                       quadratic terms
   A(i,4  ) = xx(i)^2       ;
   A(i,5  ) = yy(i)^2       ;
   A(i,6  ) = xx(i) * yy(i) ;
end

%              Print out the matrix. (Remove % before A.)
%  A

%              Number of columns of least squares matrix.
M = length(A(1,:)) 

%              Get the QR decomposition of A.
%              Q is an orthogonal matrix.
%              R is upper triangular.

[Q,R] = qr(A) ;

%              Invert the Q matrix.
%              Q' is the transpose of Q, which is also
%              the inverse since Q is orthogonal.

bb = Q'*zz   ;

%              Get the top M by M matrix from R
RR = R( 1:M, : )   ;

%              Solve RR * c = bb
%              Take only the top M entries of bb.

coef = RR\ bb(1:M)  

%              Compute the residual, method 1
%              Compute the norm of the last N - M entries
%              of bb.  These are the equations that can not
%              be solved.

r1 = norm( bb(M + 1: N ) ) 

%              Compute the residual, method 2
%              Find the norm of the vector of the difference
%              of zz and the fit.

r2 = norm( zz - A* coef ) 

%              r1 should be equal to r2

%              Prepare to plot the surface.
%              (xs, ys, zs) are mesh points on the surface.

for i = 0:10 
    x = i/10 ;
    for j = 0:10
        y = j/10 ;
        xs(i+1,j+1) = x ;
	ys(i+1,j+1) = y ;

%                                 ...  is for continuation to next line.
        zs(i+1,j+1) =  coef(1)     + coef(2)*x   + coef(3)*y  ...
                    +  coef(4)*x*x + coef(5)*y*y + coef(6)*x*y   ;
    end
end

%              Plot the surface.
mesh( xs, ys, zs)

%              Hold this plot to have the next plots on as well.
hold on

%              Plot the original data points.
plot3( xx, yy, zz, 'o')

%              Put a grid in for easier viewing.
grid on

%              Release the plot.
hold off

%              Use 
%  plot(xx, yy, 'x')
%              to see the placement of the data points in the (x,y) plane.

%              To send the plot to a file (as a postscript file)
%              use the print command.  See spline_fit.m
