function [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = smthseq(dzx,dzy,dzxx,dzyy,dzxy)
% function [K,H,k1,k2,a1,b1,g1,a2,b2,g2] = smthseq(dzx,dzy,dzxx,dzyy,dzxy)
% Calculates Gaussian (K) and mean curvature (H), most positive and
% least positve 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 + z(x,y)k),
% where elevations (z) are given at points on a regular x,y grid.
% The code is based on course notes by Smith and Sequin (S&S):
% http://www.cs.berkeley.edu/%7Esequin/CS284/TEXT/diffgeom.pdf
% These notes present a particularly clear and complete treatment of the
% eigenvectors (directions of principal curvatures), something lacking in
% most treatments I have seen.  The principal directions should be 
% orthogonal, but this is hard to see in most treatments.
%
% A key part of understanding the principal directions is understanding
% the reference frame they are described in.
% For example, classical treatments (e.g., Struik, 1961; Lipschutz, 1969)
% deal with the projections of the principal directions onto a parameter
% plane (commonly a horizontal plane in frame 'a' below), but these
% projected directions in the parameter plane usually are not orthogonal.
% Additionally, because the eigenvalues and eigenvectors are obtained from
% separate equations, one must match the correct the correct eigenvalue to
% the correct eigenvector.  So this approach is prone to error.
% The equations of Struik on p. 74,75,80,81, are still helpful to .
% With the shape operator approach of Gray (1991), the principal curvature
% directions at a point are calculated in the tangent plane to the surface
% using basis vectors in frame 'b'that, in general, are not orthogonal, and
% the matrix of shape operator coefficients is not symmetric.  
% Some kind of conversion is required in both the classical and shape
% operator treatments to yield principal curvatures directions
% that are orthogonal in both the tangent plane and in three dimensions.
% Smith and Sequin provide a clear treatment of that.  They relate results
% in frames 'b' and 'c'.  In their frame 'c' formalation the reference
% frame is orthonormal, the matrix for which the eigenvectors are found is
% symmetric, so the eigenvectors must be orthogonal; this makes sense.
% Important note: eq. 54 of S&S for matrix II_shat contains a key error:
% the (1,1) term needs to be multiplied by L; this is corrected here.
% This code passed several tests on 6/20/07 with parabolic cylinders, a
% monkey saddle, and an "egg carton" surface (z = cos(x)*cos(y)).
%
% Three reference frames come into play:
% a) Global orthonormal frame: x (east), y(north), z (up)
% b) Local frame dZ/dx (easterly), dZ/dy(northerly), n (surface normal)
% c) Local orthonormal frame: s (easterly), t(tangent), n (surface normal)
% (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. 
% (5)  N = (dZ/dx)x(dZ/dx) = -i(dz/dx) -j(dz/dy) + k. 
% (6)  n = N/|N| = (-i(dz/dx)-j(dz/dy)+k)/sqrt(1+dzdx^2+dzdy^2)
% (7)  s = dZ/dx/|dZ/dx| = (-i(dz/dx)-j(dz/dy)+k)/(1+dzdx^2+dzdy^2)
% (8)  t = (n)x(s)
% In (5), (6), and (8) the 'x' between parenthese [)x(] means cross product
% A final point deals with notation of the coefficients of the fundamental
% forms.  In nearly all treatments I have seen, the coefficients of the
% first fundamental form are E,F, and G.  For the coefficients of the
% second fundametal form, either e, f, and G or L, M, and N are used.  
% This code uses e,f, and g.
%
%
% 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 Form coefficients
% Pollard & Fletcher (2005), among others, use L,M,N instead of e,f,g.
% See Monge Patch formulas (Struik, eq. 2-7, p. 58) and (1) and (2) above.
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);
rootE = sqrt(E);        % rootE = length of vector dZx.
EG_F2 = E.*G - F.^2;    % EG_F2 = dot (dZx x dZy, dZx x dZy)
root1 = sqrt(EG_F2);    % root1 is the length of N; see eq. (5) above

% Establish the orthonormal reference frame unit basis vectors n, s, and t.
% The xyz components of these unit vectors are direction cosines.
% Find the unit normal (n) to the surface.  See eq. (6) above
nx = -dzx./root1;      ny = -dzy./root1;     nz = 1./root1;  
% Find the unit basis tangent vector s = dZx/|dZx|; see eq. (7) above.
% s =  dot(dZx/|dZx|,i)i + 0j +  dot(dZx/|dZx|,k)k
% Its components [dot(s,i), dot(s,j), dot(s,k)] are its direction cosines.
sx = 1./rootE;          sy = zeros(m,n);        sz = dzx./rootE;	
% Find the unit basis tangent vector t = n x s by crossing n and s.
tx = ny.*sz - nz.*sy;   ty = nz.*sx - nx.*sz;   tz = nx.*sy - ny.*sx; 

% Set coefficients of the Second (e,f,g) Fundamental Form.  
% Smith and Sequin use (L,M,N) instead of (e,f,g).
% See eq. 5-9 on p. 75 of Struik.
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 (1991) 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 & principal curvatures
% First form the symmetric matrix II_shat of Smith and Sequin, eq. 54
% This gives matrix II components in the orthonormal stn reference frame
% NOTE! The entry in the first row, first column of II_shat should be
% e./E instead of 1/E.  This corrects the error in eq. (54) of S&S.
II_shat_11 = e./E;      
II_shat_12 = (f - (F.*e./E))./root1;
% II_shat_21= II_shat_12;
II_shat_22 = (1./EG_F2).*( (F.^2.*e./E) - 2*F.*f + E.*g);
% Find the eigenvectors and eigenvalues of the symmetric matrix II_shat
% The eigenvector components of II_shat are in the tangent plane of s and t
% NOTE: eigenvectors returned by eig2d are of unit length.
[k1,k2,e1s,e1t,e2s,e2t] = ...
  eig2d(II_shat_11(:),II_shat_12(:),II_shat_12(:),II_shat_22(:));
% Reshape the values returned from eig2d
k1 = reshape(k1,m,n);	e1s = reshape(e1s,m,n);     e1t = reshape(e1t,m,n);
k2 = reshape(k2,m,n);	e2s = reshape(e2s,m,n);     e2t = reshape(e2t,m,n);

% Finally, rotate the principal directions from the stn frame to xyz.
% The n-components of the eigenvectors are zero and don't contribute.
a1 = sx.*e1s + tx.*e1t;     % cosine of angle between e1 and x-axis  
b1 = sy.*e1s + ty.*e1t;     % cosine of angle between e1 and y-axis   
g1 = sz.*e1s + tz.*e1t;     % cosine of angle between e1 and z-axis
a2 = sx.*e2s + tx.*e2t;     % cosine of angle between e2 and x-axis   
b2 = sy.*e2s + ty.*e2t;     % cosine of angle between e2 and y-axis   
g2 = sz.*e2s + tz.*e2t;     % cosine of angle between e2 and z-axis