function [fitresult,gof] = KPgSurfaceFit(xK, yK, zK)
%CREATESURFACEFIT(XK,YK,ZK)
%  Fit surface to data.
%
%  Data for 'untitled fit 1' fit:
%      X Input : xK
%      Y Input : yK
%      Z output: zK
%      Weights : (none)
%
%  Output:
%      fitresult : an sfit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, SFIT.

%  Auto-generated by MATLAB on 31-Jul-2013 11:26:46
%  Edited by Tom Tobin shortly afterward to make a standalone m-file
%  Added raw data for constructing K-Pg surface
%  Requires sftool.m
%
%  HOW TO USE:
%  1. Place this file in your working matlab directory (or possibly set
%  path to it)
%  2. Type [model,gof]=KPgSurfaceFit()
%  3. You will see a plot of the surface fit that you can play with 
%  4. 
%  5. use command: predint(model,[X,Y], 0.95, 'Observation') to find vertical 95% confidence
%  interval bounds for the predicted surface height at any given location
%  described by X,Y as easting and northing coordinates in UTM 13N - check
%  that site is reasonably within areal extent of controlling points
%  6. use model(X,Y) to get calculated surface value for that location
%
%  N.B. Technically some trig should be used to project your sampling
%  location perpendicular to the stratigraphic plane, but due to the
%  shallowness of the dip (0.2 deg) this actually ends up being irrelevant
%
%  methods(model)
%

format long e 

%Overrides any other x,y,z data to be the K-Pg surface data
east_x=[    3.550258530000000e+005
    3.239850780000000e+005
    3.366234050000000e+005
    3.367136130000000e+005
    3.378606570000000e+005
    3.253545950000000e+005
    3.252666400000000e+005
    3.865441540000000e+005
    3.865783770000000e+005
    3.482130770000000e+005
    3.540239730000000e+005
    3.828867180000000e+005
    3.489263380000000e+005
    3.933583190000000e+005
    3.595201380000000e+005];
north_y=[    5.266644810000000e+006
    5.273757185000000e+006
    5.267092774000000e+006
    5.266620670000000e+006
    5.266777047000000e+006
    5.276110774000000e+006
    5.276781925000000e+006
    5.283054277000000e+006
    5.283046457000000e+006
    5.266129535000000e+006
    5.265291846000000e+006
    5.278011040000000e+006
    5.265778457000000e+006
    5.264212543000000e+006
    5.268269448000000e+006];
elev_z=[    8.312411521000000e+002
    9.185127498000001e+002
    8.770822933000000e+002
    8.772815448000000e+002
    8.811629039000000e+002
    9.273169384000000e+002
    9.286120022000000e+002
    8.043067359000000e+002
    8.118614560000000e+002
    8.462170678000000e+002
    8.360298995000001e+002
    8.038597923000000e+002
    8.439968771000000e+002
    7.368498737000000e+002
    8.250231635000000e+002];

xK=east_x;
yK=north_y;
zK=elev_z;
% End added data

%% Fit: 'untitled fit 1'.
ft = fittype( 'poly11' );
opts = fitoptions( ft );
opts.Robust = 'Bisquare';
opts.Weights = zeros(1,0);
[fitresult, gof] = fit( [xK, yK], zK, ft, opts );

% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, [xK, yK], zK, 'Style', 'PredObs' );
grid on
% Label axes
xlabel( 'xK' );
ylabel( 'yK' );
zlabel( 'zK' );
legend( h, 'untitled fit 1', 'zK vs. xK, yK', 'Location', 'NorthEast' );