public static void functionCircleCircle(double[] x, ref double func, object obj)
{
    Circle[] circles = (Circle[])obj;
    double theta = x[0];
    double phi = x[1];
    Vector w1 = Math.Cos(theta) * circles[0].U + Math.Sin(theta) * circles[0].V;
    Vector w2 = Math.Cos(phi) * circles[1].U + Math.Sin(phi) * circles[1].V;
    double d = (circles[1].Center - circles[0].Center + circles[1].Radius * w2 - circles[0].Radius * w1).Norm;
    func = d * d;
}

// define the circles
// first case
Circle c1 = new Circle(new Point(0, 0, 0), 10, new Vector(0, 0, 1));   // black circle
Circle c2 = new Circle(new Point(0, 0, 30), 10, new Vector(1, 0, 0));   // green circle

// second case
//Circle c1 = new Circle(new Point(0, 0, 0), 10, new Vector(0, 0, 1));   // black circle
//Circle c2 = new Circle(new Point(0, 0, 30), 10, new Vector(0, 1, 0));   // red circle

double[] x = new double[] { 0, 0 };   // initial guess
double epsg = 1.0e-6;
double epsf = 0;
double epsx = 0;
double diffstep = 1.0e-6;
int maxits = 0;
alglib.minlbfgsstate state;
alglib.minlbfgsreport rep;
alglib.minlbfgscreatef(1, x, diffstep, out state);
alglib.minlbfgssetcond(state, epsg, epsf, epsx, maxits);
alglib.minlbfgsoptimize(state, functionCircleCircle, null, circles);
alglib.minlbfgsresults(state, out x, out rep);

Point p1 = circles[0].Eval(x[0]);   // eval optimized angle on first circle
Point p2 = circles[1].Eval(x[1]);   // eval optimized angle on second circle
LineSegment segment = new LineSegment(p1, p2);   // this *should* be the shortest line between the two circles

// Output
Console.WriteLine("Theta = {0}\tPhi = {1}", x[0], x[1]);
Console.WriteLine("Length = {0}", segment.Length);

DrawCircle(circles[0]);
DrawCircle(circles[1]);
DrawLineSegment(segment);