% written by Jay Chapman, 2009?
% Code to find best-fit plane to a series of x,y,z points

% Can cite with:

%  Chapman, J.B., Worthington, L.L., Pavlis, T.L., Bruhn, R.L.,
%       and Gulick, S.P.S., 2011, The Suckling Hills fault, Kayak
%       Island zone, and accretion of the Yakutat microplate, Alaska:
%       Tectonics, v.30.

clear all
close all

%  ********  Text File Import *************
%format of text file is 1st column=x, 2nd column=y, 3rd column = z
%X and Y data in UTM coordinates, Z is elevation in meters
% no headers in the text file

[filename,pathname] = uigetfile('*.txt','Select XYZ point data');
full = fullfile(pathname,filename);
A = importdata(full);  % import data

B = ones(length(A),3);  %matrix of all ones
B(:,1) = A(:,1);  %add in X and Y data
B(:,2) = A(:,2);

Z = A(:,3);  %elevation data

u = B\Z;  %perform inversion
mx = u(1,:);       % slope of x
my = u(2,:);      % slope of y
z0 = u(3,:);      % elevation of plane at origin

         
% **Equation:    %g(x) + %g(y) + %g(z) = 0
%     g(x)=mx,  g(y)=my,  g(z)=z0
 
dp = atan(sqrt((mx*mx)+(my*my)));  %dip calculation
dip = dp*(180/pi);  %radians
 
st  = atan2(mx,my);         % strike calculation (true dip +90)
str = st*(180/pi)+90;       % gives degrees from north before +90
 
if (str < 0)
    strk = 360+str;
else
    strk = str;
end;   

fprintf(...                 %display results
'\n  Strike and Dip:\n  %g   %g \n\n'...
    ,strk,dip)