%Sheep Mountain Anticline Example
close all
clear
clc

% Load Data
load SMA_rot_50m.mat
dx=50;
dy=50;
[n,h]=size(Z);
Z(1)=609.12;

%Spectral Filtering to Remove Small Scale Surface Noise
Z2=[Z fliplr(Z)]; Z3=[Z2; flipud(Z2)];
[m,o]=size(Z3);
fx=(1:(o/2))/o/dx;
fx=[0, fx, -fliplr(fx(1:end-1))];
fy=(1:(m/2))/m/dy;
fy=[0, fy, -fliplr(fy(1:end-1))];
A=fft2(Z3);
[Fx,Fy]=meshgrid(fx,fy);
clear Z2 Z3
s=.002;	% Spatial Frequency Threshold
imask=sqrt(Fx.^2+Fy.^2)<s;
As=imask.*A;
Zs=real(ifft2(As));
A=fft2(Zs);
Zs=Zs(1:n,1:h);
clear A Fx Fy Ar As
Z=Zs;

Kt=0;
warning off
curvature
warning on
sizeZ=n*h;
space=10;
sizedec=sizeZ/space;
index=(1:sizedec)*space;

% Plot Kmax Magnitudes and Orientations
figure
surf(X,Y,Z,Kmax);
caxis([min(min(Kmin)) max(max(Kmax))]);
Xquiv=X(index);
Yquiv=Y(index);
Zquiv=Z(index)+15;
Kmaxxquiv=Kmaxx(index);
Kmaxyquiv=Kmaxy(index);
Kmaxzquiv=Kmaxz(index);
hold on
quiver3(Xquiv,Yquiv,Zquiv,Kmaxxquiv,Kmaxyquiv,Kmaxzquiv,'k.');
hold off
colorbar

% Plot Kmin Magnitudes and Orientations
figure
surf(X,Y,Z,Kmin);
caxis([min(min(Kmin)) max(max(Kmax))]);

Kminxquiv=Kminx(index);
Kminyquiv=Kminy(index);
Kminzquiv=Kminz(index);
hold on
quiver3(Xquiv,Yquiv,Zquiv,Kminxquiv,Kminyquiv,Kminzquiv,'k.');
hold off

% Plot Geologic Curvature with Various Curvature Thresholds

figure
surf(X,Y,Z,GC);
caxis([-4 3]);
view([-83,11])


Kt=0.0001;
warning off
curvature
warning on

figure
surf(X,Y,Z,GC);
view([-83,11])
caxis([-4 3]);


Kt=0.0005;
warning off
curvature
warning on

figure
surf(X,Y,Z,GC);
view([-83,11])
caxis([-4 3]);


Kt=0.00105;
warning off
curvature
warning on

figure
surf(X,Y,Z,GC);
view([-83,11])
caxis([-4 3]);