%Jay Chapman
%3 point problem solver

% use UTM coordinate system where x,y,z are in meters

disp('**  Enter points as [x y z]')
disp('**  Where east, north and up are positive')
disp('')

pt1 = input('First Point:');
pt2 = input('Second Point:');
pt3 = input('Third Point:');

A = [pt1(1),pt1(2),1;pt2(1),pt2(2),1;pt3(1),pt3(2),1];
z = [pt1(3);pt2(3);pt3(3)];

u = A\z;
mx = u(1,:);      % slope of x
my = u(2,:);      % slope of y
z0 = u(3,:);      % elevation of plane at origin

%display results
fprintf(...             
'\n  **Equation:    %g(x) + %g(y) + %g(z) = 0 \n\n'...
    ,mx,my,z0)

dp = atan(sqrt((mx*mx)+(my*my)));  %dip calculation
dip = dp*(180/pi);

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:*  %g   %g \n\n'...
    ,strk,dip)