%Calculate First Partials (dZx, dZy)
[dZx,dZy]=gradient(Z,dx,dy);

%Calculate Second Partials (dZxx, dZxy, dZyy)
[dZxx,dZxy]=gradient(dZx,dx,dy);
[dZxy,dZyy]=gradient(dZy,dx,dy);

%Calculate First Fundamental Coefficiants (E, F, G)
E=1+dZx.^2;
F=dZx.*dZy;
G=1+dZy.^2;

%Calculate Surface Unit Normal Componants (Nx, Ny, NZ)
normN=-sqrt(dZx.^2 +dZy.^2+1);
Nx=-dZx./normN;
Ny=-dZy./normN;
Nz=ones(size(dZx))./normN;

%Calculate Second Fundamental Coefficiants (L, M, N)
L=dZxx./normN;
M=dZxy./normN;
N=dZyy./normN;

%Calculate Prinicpal Curvature Directions and Magnitudes
% (Kminx, Kminy, Kminz, Kmin; Kmaxx, Kmaxy, Kmaxz, Kmax)
[n,h]=size(Z);
Km1=zeros(n,h);
Km2=zeros(n,h);
K1x=zeros(n,h);
K1y=zeros(n,h);
K2x=zeros(n,h);
K2y=zeros(n,h);

for i=1:n
    for j=1:h
        I=[E(i,j) F(i,j);F(i,j) G(i,j)];
        II=[L(i,j) M(i,j);M(i,j) N(i,j)];
        SO=inv(I)*II;
        A=max(max(isnan(SO)));
        if A==0
            [Kd,Km]=eig(SO);
        elseif A==1
            Kd=[NaN NaN;NaN NaN];
            Km=[NaN NaN;NaN NaN];
        end
        %Principle Curvature A
        KmA(i,j)=Km(1,1); KdA1(i,j)=Kd(1,1); KdA2(i,j)=Kd(2,1);
        %Principle Curvature B
        KmB(i,j)=Km(2,2); KdB1(i,j)=Kd(1,2); KdB2(i,j)=Kd(2,2);    
    end
end


% Convert Principle Directions to Tangent Space (3D)
KAx=KdA1;
KAy=KdA2;
KAz=KdA1.*dZx+KdA2.*dZy;
normA = sqrt(KAx.^2+KAy.^2+KAz.^2);
KAx = KAx./normA;
KAy = KAy./normA;
KAz = KAz./normA;

KBx=KdB1;
KBy=KdB2;
KBz=KdB1.*dZx+KdB2.*dZy;
normB = sqrt(KBx.^2+KBy.^2+KBz.^2);
KBx = KBx./normB;
KBy = KBy./normB;
KBz = KBz./normB;

%Sort to Kmax and Kmin
i=find(KmB>KmA);
Kmax=KmA;
Kmin=KmB;
Kmaxx=KAx; Kmaxy=KAy; Kmaxz=KAz;
Kminx=KBx; Kminy=KBy; Kminz=KBz;

Kmax(i)=KmB(i);
Kmin(i)=KmA(i);
Kmaxx(i)=KBx(i); Kmaxy(i)=KBy(i); Kmaxz(i)=KBz(i);
Kminx(i)=KAx(i); Kminy(i)=KAy(i); Kminz(i)=KAz(i);

%Apply Curvature Threshold
i=find(abs(Kmin)<Kt);
Kmin(i)=0;
i=find(abs(Kmax)<Kt);
Kmax(i)=0;

% Calculate Gaussian Curvature
Kg=Kmax.*Kmin;
i=find(Kg<0);
j=find(Kg>0);
Kg(i)=-1.*ones(size(i));
Kg(j)=ones(size(j));

% Calculate Mean Curvature
Km=0.5.*(Kmax+Kmin);
i=find(Km<0);
j=find(Km>0);
Km(i)=-1.*ones(size(i));
Km(j)=ones(size(j));

% Find Saddles
Kst=Km*2;
i=find(abs(Kst)<Kt);
	
%Calculate Geologic Curvature (GC)
% -4=perfect saddle -2.5=dome -1.5=antiform -0.5 antiformal-saddle
% 0.2=plane 0.75=synformal saddle 2=synform 3=basin
GC=Kg;
%Synformal Saddle (0.75)
i=find(Kg<0 & Km<0);
GC(i)=0.75.*ones(size(i));
%Saddle (-4)
i=find(Kg<0 & abs(Kmax+Kmin)<Kt);
GC(i)=-4.*ones(size(i));
%Antiformal Saddle (-0.5)
i=find(Kg<0 & Km>0);
GC(i)=-0.5.*ones(size(i));
%Dome (-2.5)
i=find(Kg>0 & Km>0);
GC(i)=-2.5.*ones(size(i));
%Plane (0.20)
i=find(Kg==0 & Km==0);
GC(i)=0.20.*ones(size(i));
%Synform (2)
i=find(Kg==0 & Km<0);
GC(i)=2*ones(size(i));
%Antiform (-1.5)
i=find(Kg==0 & Km>0);
GC(i)=-1.5.*ones(size(i));
%Basin (3)
i=find(Kg>0 & Km<0);
GC(i)=3.*ones(size(i));

% clear Kg Km i j KmA KmB KBx KBy KBz KAx KAy KAz Km1 Km2 K1x K1y K2x K2y ...
%     normN I II A SO Kd KdA1 KdA2 KdB1 KdB2