% 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