The following sample code demonstrates how to do the above analytically.

C++

VB.NET (Rhino 4)

C# (Rhino 4)

 Public override IRhinoCommand.result RunCommand(IRhinoCommandContext context)
 {
   MRhinoGetObject go = new MRhinoGetObject();
   go.SetCommandPrompt("Select Surface to evaluate");
   go.SetGeometryFilter(IRhinoGetObject.GEOMETRY_TYPE_FILTER.surface_object);
   IRhinoGet.result res = go.GetObjects(1, -1);
   if (res == IRhinoGet.result.@object)
   {
     IRhinoObjRef objref = go.Object(0);
     IOnSurface surf = objref.Surface();
     if (surf != null)
     {
       double u = 0; double v = 0; double gauss = 0; double mean = 0;
       // the magnitude of the max and min curvatures
       double k_max = 0; double k_min = 0;
       On3dVector du = new On3dVector(); On3dVector dv = new On3dVector(); 
       On3dVector duu = new On3dVector(); On3dVector duv = new On3dVector(); 
       On3dVector dvv = new On3dVector(); On3dVector n = new On3dVector();
       // the principal directions
       On3dVector Kmax = new On3dVector(); On3dVector Kmin = new On3dVector();  


       // select a point on the surface by constraining GetPoint to the surface 
       MRhinoGetPoint gp = new MRhinoGetPoint();
       gp.Constrain(surf);
       gp.GetPoint();
       if (gp.CommandResult()!= IRhinoCommand.result.success)
         return gp.CommandResult();


       On3dPoint pt_start = gp.Point();
       if (surf.GetClosestPoint(pt_start, ref u, ref v) == false)
         return IRhinoCommand.result.failure;
       if (surf.Ev2Der(u, v, ref pt_start, ref du, ref dv, ref duu, ref duv, ref dvv)==false)
         return IRhinoCommand.result.failure;


       // find the unit normal at the selected point
       n = OnUtil.ON_CrossProduct(du, dv);
       n.Unitize();


       //  Do the calculation (below)
       EvPrincipalCurvature( du, dv, duu, duv, dvv, n,
                             ref gauss,
                             ref mean,
                             ref k_max,ref k_min,ref Kmax,ref Kmin);


       // Create the end points for the principal curvature lines 
       // from the start point + the direction vector
       On3dPoint[] pt_K = new On3dPoint[2];
       for (int i = 0; i < 2; i++)
         pt_K[i] = new On3dPoint (pt_start);
       On3dPoint.op_AdditionAssignment(pt_K[0],Kmax);
       On3dPoint.op_AdditionAssignment(pt_K[1], Kmin);


       // Create the lines representing the direction vectors
       OnLine[] princDir = new OnLine[2];
       for (int i = 0; i < 2; i++)
       {
         princDir[i] = new OnLine(pt_start, pt_K[i]);
         context.m_doc.AddCurveObject(princDir[i]);
       }
       RhUtil.RhinoApp().Print( "Gaussian curvature = "  + gauss
                                 + "  Mean curvature = " + mean+"\r\n");
       context.m_doc.Redraw();
     }
   }
   return IRhinoCommand.result.success;
 }
 /// <summary>
 /// solves for the principal curvatures and directions of those curvatures
 /// returns false if the point is an Umbillic and the Principal Directions are not unique
 /// </summary>
 /// <param name="du">first partial</param> 
 /// <param name="dv">first partial</param> 
 /// <param name="duu">second partial</param> 
 /// <param name="duv">second mixed partial</param> 
 /// <param name="dvv">second partial</param> 
 /// <param name="n">normal vector</param> 
 /// <param name="gauss">(out)gaussian curvature</param> 
 /// <param name="mean">(out)mean curvature</param> 
 /// <param name="k1">(out)maximum principal curvature</param> 
 /// <param name="k2">(out)minimum principal curvature</param> 
 /// <param name="e1">(out)direction of max curvature</param> 
 /// <param name="e2">(out)direction of min curvature</param> 
 /// <returns>true if the point is not an umbillic</returns>
 private bool EvPrincipalCurvature(
    On3dVector du, On3dVector dv, 
    On3dVector duu, On3dVector duv, On3dVector dvv, On3dVector n,
    ref double gauss, ref double mean,
    ref double k1, ref double k2, ref On3dVector e1, ref On3dVector e2
    )
 {
   double tol = .0000001;  // tolerance for testing the Umbillic
   n.Unitize();            // ensure that we are using the unit normal
   double E = du * du;     // components of the first fundamental form
   double F = du * dv;
   double G = dv * dv;
   double L = duu * n;     // components of the second fundamental form
   double M = duv * n;
   double N = dvv * n;


   // calculate the gaussian and mean curvatures
   gauss = (L * N - M * M) / (E * G - F * F);                  
   mean = (G * L - 2 * F * M + E * N) / (2 * (E * G - F * F));  


   // calculate the principal curvatures (k1 = max curvature)
   k1 = mean + Math.Sqrt(mean * mean - gauss);
   k2 = mean - Math.Sqrt(mean * mean - gauss);
   if (Math.Abs(k1 - k2) < tol)  // test if the point is an Umbillic
       return false;
   if (k1 < k2)
   {
     double temp = k1;
     k1 = k2;
     k2 = temp;
   }


   // calculate the components of the principal directions in the tangent space
   double[] a = new double[2] { (k1 * F - M) / (L - k1 * E),
                                (k2 * F - M)/(L - k2 * E) };
   double[] b = new double[2] { 1.0, 1.0 };
   for (int i=0; i<2; i++){        // normalize them
     double mag = Math.Sqrt(a[i] * a[i] + b[i] * b[i]);
     a[i] = a[i] / mag;
     b[i] = b[i] / mag;
   }


   // create the principal direction unit vectors
   e1 = du*a[0] + dv*b[0];
   e1.Unitize();
   e2 = du * a[1] + dv * b[1];
   e2.Unitize();
   return true;
   }
 }