% Jay Chapman, 2016
        % modified from a code first put together by 
        % Dr. Ron Bruhn (emeritus, Univ. of Utah) and his students

% Can be cited with:
%    Chapman, J.B., Pavlis, T.L., Bruhn, R.L., Worthington,
%           L.L., Gulick, S.P.S., and Berger, A.L., 2012, Structural
%           relationships in the eastern syntaxis of the St. Elias
%           orogen, Alaska: Geosphere, v.8, p.105-126.
        
 
clear all
close all


% ***** Input DEM Resolution Here  *******
posting = 30;   %DEM posting (meters)
                % generally the same as horiz. resolution or pixel size
                % common resolutions below:
                % 1 arc second = 30m
                % 1/3 arc second = 10m

                
 % *** no need to edit any other code below for specific DEMs ******
                
[fn,pn] = uigetfile('*.mat','DEM file: ');
eval(['load ' fn]);

dem = AAA;
[nrow,ncol] = size(dem);

east = 0:posting:posting*(ncol-1);
north = 0:posting:posting*(nrow-1);

[easting,northing] = meshgrid(east,north);

Hfig1 = figure(1);
surf(easting,northing,dem); %plots dem elevation values as surface
view(0,90); %set the dem view to "map view"
axis equal; %sets data units equal in x and y directions
shading interp; %smooths color transition in dem
xlabel('Easting (m)')
ylabel('Northing (m)')
hold on
contour3(easting,northing,dem,10,'k'); %adds contour lines to figure
hold off


%User chooses a point on dem to have plane pass through
 figure(1)
 [x,y] = ginput(1);
 xpt = x;
 ypt = y;
 zpt = interp2(easting,northing,dem,x,y); %interprets elevation from dem

 
 
%Prompt user for the strike and dip of the plane, right hand rule
prompt={'Enter strike (right hand rule): ','Enter dip:'};
   name='Intersecting Plane';
   numlines=1;
   defaultanswer={'0','0'};
   
  answer=inputdlg(prompt,name,numlines,defaultanswer);
  strike = str2num(answer{1});
  dip = str2num(answer{2});


% Use the strike and dip to compute the apparent dips in the east and north
% directions. Note that dip angles are signed to indicate east (+) vs west (-) and
% north (+) vs south (-)

% Calculates the dip direction from the strike azimuth
% East/North dip angles positive, West/South = negative
   if strike <= 270
       dip_dir = strike + 90;
   else
       dip_dir = strike + 90 - 360;
   end

 % Calculates the angle between dip direction and east and north
 % axes along with the sign of the east and west components.
if dip_dir <= 90
    dip_trig_angle = 90 - dip_dir;
elseif dip_dir <= 180 & dip_dir > 90
      dip_trig_angle = 270+180-dip_dir;
elseif dip_dir <=270 & dip_dir > 180
     dip_trig_angle = 360+90-dip_dir;
else
    dip_trig_angle = 360+90-dip_dir;
end
 
% Take cos/sin to determine apparent dip
east_cos = cosd(dip_trig_angle);
north_cos = sind(dip_trig_angle);
adip_east = atand(tand(dip)*east_cos);
adip_north = atand(tand(dip)*north_cos);


% Compute the intersection of the plane on the east and north axis
east_intercept = xpt + zpt/tand(adip_east);
north_intercept = ypt + zpt/tand(adip_north);



% Define 3 points in vector that make up plane
% selected point=row 1, east axis point=row 2, north axis point = row3
points = [xpt ypt zpt; east_intercept ypt 0; xpt north_intercept+.001 0];


% Calculate coefficients (ABCD) of plane equation (Ax+By+Cz = D)
po = points(1,:);
p1 = points(2,:);
p2 = points(3,:);

a1 = p1(1) - po(1);
a2 = p1(2) - po(2);
a3 = p1(3) - po(3);

b1 = p2(1) - po(1);
b2 = p2(2) - po(2);
b3 = p2(3) - po(3);

A = (a2*b3 - a3*b2);
B = (b1*a3 - a1*b3);
C = (a1*b2 - a2*b1); 
D = (a3*b2 - a2*b3)*po(1) + (a1*b3 - a3*b1)*po(2) + (b1*a2 - a1*b2)*po(3);

%Calulate the surface of the plane (like a DEM for the plane)
zg = -(A*easting + B*northing + D)/C;


%Calculate difference between plane and DEM, where difference = 0 is where
%   the two surfaces intersect
zr = zg - dem;

%Plot intersection on original figure
  figure(1)
  hold on
 contact = [0 5]; % The contour intervals that bound the contact within 5 meters.
[cresidual,Hresidual] = contour3(easting,northing,zr,contact);
[zres] = interp2(easting,northing,dem,cresidual(1,:),cresidual(2,:));
Hplot_contact = plot3(cresidual(1,:),cresidual(2,:),zres+2,'.w');
  hold off
  
  
 


