function [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = struik(dzx,dzy,dzxx,dzyy,dzxy)
% function [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = struik(dzx,dzy,dzxx,dzyy,dzxy)
% Calculates Gaussian (K) and mean curvature (H), most positive and
% least positive principal curvatures (k1 and k2, respectively), and
% the associated direction cosines for k1 and k2 (a1,b1,g1, and 
% a2, b2, and g2, respectively) of a surface Z (Z = xi + yj + zk).
% This code relies upon a classical approach (e.g., Struik, 1961).
% See equations of Struik (1961), p. 74,75,80,81
% The eigenvector trends are found and then used to obtain the direction
% cosines of the eigenvectors in three dimensions.
% This code passed several tests on 6/21/07 with parabolic cylinders, a
% monkey saddle, and an "egg carton" surface (z = cos(x)*cos(y)).
%
% (1)  Z = xi + yj + z(x,y)k, where big Z is the surface to be analyzed
% (2)  dZ/dx = i + (dz/dx)k.  This is a tangent that trends east.
% (3)  dZ/dy = j + (dz/dy)k.  This tangent trends north. 
% (4)  d2Z/dx2 = (d2z/dx2)k,  d2Z/dy2 = (d2z/dy2)k, d2Z/dxdy = (d2z/dxdy)k. 
%
% Input parameters (Do not confuse big Z with little z) 
% dzx,dzy = first partial derivatives of z with respect to x and y
% dzxx,dzyy,dzxy =  second part. derivs. of z with respect to x, xy, and y  
%
% Output parameters 
% K = Gaussian curvature (k1*k2)
% H = Mean curvature (H) (k1 + k2)/2
% k1 = most positive principal curvature
% k2 = most negative principal curvature
% a1,b1,g1 = direction cosines for k1 with respect to the x,y,z axes
% a2,b2,g2 = direction cosines for k2 with respect to the x,y,z axes
%       % See online notes for GG303 for more on direction cosines
%
% Examples
% [X,Y] = meshgrid(-2:0.1:2,-1:0.1:1);
% % Test 1: Parabolic cylinder parallel to x-axis z = x.^2
% z = X.^2; dzx = 2.*X; dzy = 0.*Y;
% dzxx = 2*ones(size(X)); dzxy =  dzy; dzyy = dzy;
% [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = struik(dzx,dzy,dzxx,dzyy,dzxy)
% % Test 2: Parabolic cylinder parallel to y-axis z = y.^2
% z = Y.^2;dzx = 0.*X; dzy = 0.*dzx; 
% dzxx = dzx; dzxy = dzx; dzyy = 2*ones(size(Y));
% [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = struik(dzx,dzy,dzxx,dzyy,dzxy)
% % Test 3: Parabolic cylinder parallel with axis in xy plane
% t = pi/4*ones(size(X));     
% cost = cos(t);  sint = sin(t);
% z = (X.*cost + Y.*sint).^2;
% dzx = 2*(X.*cost + Y.*sint).*cost; dzy = 2*(X.*cost + Y.*sint).*sint;  
% dzxx = 2*cost.*cost; dzxy = 2*sint.*cost; dzyy = 2*sint.*sint;
% [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = struik(dzx,dzy,dzxx,dzyy,dzxy)
% % Test 4: Monkey saddle z = x.^3 - 3xy.^2;
% z = X.^3 - 3*x.*y.^2;
% dzx = 3*X.^2 - 3*Y.^2; dzy = -6*X.*Y;
% dzxx = 6*X; dzxy = -6*Y; dzyy = -6*X;
% % Test 5: Egg carton z = cos(x).*cos(y);
% z = cos(X).*cos(Y);
% dzx = -sin(X).*cos(Y); dzy = -cos(X).*sin(Y);
% dzxx = -cos(X).*cos(Y); dzxy =  sin(X).*sin(Y); dzyy = -cos(X).*cos(Y);

[m,n] = size(dzx);

% Calculate First (E,F,G) and Second (e,f,g) Fundamental Coefficients
% Pollard & Fletcher (2005), among others, use L,M,N instead of e,f,g.
% See Monge Patch formulas of Struik, eq. 2-7, p. 58, and eq. 5-9 on p. 75.
E = 1 + dzx.^2;	        % E = dZx*dZx (dot product of dZx with itself);
F =     dzx.*dzy;	    % F = dZx*dZy (dot product of dZx with dZy);
G = 1 + dzy.^2;	        % G = dZy*dZy (dot product of dZy with itself);
EG_F2 = E.*G - F.^2;
root1 = sqrt(EG_F2);    % root1 is the length of (dZ/dx) x (dZ/dy)
e = dzxx./root1;        f = dzxy./root1;        g = dzyy./root1;

% Calculate Gaussian (K) and mean (H) curvatures.
% See eq. 7-2 and 7-3 of Struik (p. 83) for K and H.
% See p. 282 and 283 of section 14.4 of Gray for k1 and k2.
K = (e.*g - f.^2)./(EG_F2);                     % K = k1.*k2;
H = (E.*g - 2*f.*F + e.*G)./(2*(EG_F2));        % H = (k1+k2)/2;
% r = sqrt(H.^ 2 - K);       % k1 = H + r;       % k2 = H - r;

% Calculate principal directions (a_,b_,g_) & principal curvatures
% tan_theta = trend of principal direction = dy/dx
[theta1,theta2,k1,k2,a1,b1,a2,b2] = initialize(m,n);
% Solve for tan_theta (=dy/dx) in the following equation
% to get the direction cosines (a_,b_) of the principal trends
% a.*tan_theta.^2 + b.*tan_theta + c = 0.  See Struik, p. 80, eq. 6-5a 
% tan_theta1 = -b + sqrt(b.^2 - 4ac)/2a
% tan_theta2 = -b - sqrt(b.^2 - 4ac)/2a
a = F.*g-G.*f;      b = E.*g-G.*e;      c = E.*f-F.*e; 
sqrtb2_4ac = sqrt(b.^2 - 4.*a.*c);

% General case: eigenvectors do not trend parallel to x and y
% For data from surfaces in nature, this should be the only case used
i = find(a~=0);	 
% Find tan_theta (= dy/dx) for directions of k1,k2.  See Struik, eq. 6-5a
theta1(i) = atan( (-b(i) + sqrtb2_4ac(i))./(2.*a(i)) );
theta2(i) = atan( (-b(i) - sqrtb2_4ac(i))./(2.*a(i)) );
a1(i) = cos(theta1(i)); 	% alpha for k1;
b1(i) = sin(theta1(i));     % beta for k1;
a2(i) = cos(theta2(i));     % alpha for k2;
b2(i) = sin(theta2(i));     % beta for k2;
% Calculate k1 and k2 using eq. 6-3 of Struik
k1(i) = (e(i)+2.*f(i).*tan(theta1(i))+g(i).*tan(theta1(i)).^2)...
      ./(E(i)+2.*F(i).*tan(theta1(i))+G(i).*tan(theta1(i)).^2);
k2(i) = (e(i)+2.*f(i).*tan(theta2(i))+g(i).*tan(theta2(i)).^2)...
      ./(E(i)+2.*F(i).*tan(theta2(i))+G(i).*tan(theta2(i)).^2);

% Special cases
% Parametric lines are local lines of curvature, (e.g., surface is locally
% cylindrical with an axis trending parallel to x or y) and umbilic points
% (where  curvatures are the same in every direction [locally spherical])
% This is where F=f=0 and so a=c=0.  
j = find(a==0);    
% Calculate k1 and k2 using eq. 6-8 of Struik
k1(j) = e(j)./E(j);  	% k1 is along the x-direction (Struik, p. 81)
k2(j) = g(j)./G(j);     % k2 is along the y-direction (Struik, p. 81) 

% The results so far do NOT require k1 >= k2.
% Check results to make sure k2 does not exceed k1 and that the proper
% direction cosines and proper principal values are matched up.
index = find(k2>k1);
% Where needed, switch the k1,k2 values; a1,a2 values; and b1,b2 values   
c2 = k1(index);     k1(index) = k2(index);      k2(index) = c2; 
A2 = a1(index);     a1(index) = a2(index);      a2(index) = A2; 
B2 = b1(index);     b1(index) = b2(index);      b2(index) = B2; 
 
% Find the eigenvector z-components in 3D by multiplying the x,y components
% by dz/dx and dz/dy, respectively, and summing.  See (2) and (3) above.
g1 = a1.*dzx + b1.*dzy;         g2 = a2.*dzx + b2.*dzy;

% Now normalize the eigenvector components
length_e1 = sqrt(a1.^2 + b1.^2 + g1.^2);
length_e2 = sqrt(a2.^2 + b2.^2 + g2.^2);
a1 = a1./length_e1;     b1 = b1./length_e1;     g1 = g1./length_e1; 
a2 = a2./length_e2;     b2 = b2./length_e2;     g2 = g2./length_e2; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [m1,m2,m3,m4,m5,m6,m7,m8] = initialize(m,n)
% function [m1,m2,m3,m4,m5,m6,m7,m8] = initialize(m,n)
% Initializes matrices to zeros(m,n) and ones(m,n)
m1 = zeros(m,n);    m2 = m1;	m3 = m1;    m4 = m1;   
m5 = ones(m,n);     m6 = m1;    m7 = m1;    m8 = m5;