Find two tangent points for given circle and two points inside the circle , Lua

I have two Points A and B (given with coordinates ax, ay and bx, by). This two points are inside of given Circle C (given with coordinates cx, cy and radius cr).

How to find points T1 and T2? Points T1 and T2 must be on the tangent between circle C and imaginary circles D1 and D2 (both are on the points A and B).

Example values:

A (30, 30)
B (40, 20)
C (50, 50, 40)

right solution:

T1 (28.6061, 16.202)
T2 (87.3939, 35.798)

I've started with the middle point M: mx, my = (ax+bx)/2, (ay+by)/2

The perpendicular line to line AB:

m1x, m1y = mx + (ay-my), my - (ax-mx)
m2x, m2y = mx + (by-my), my - (bx-mx)

And how to get the circles D1 and D2?

You are almost there, however the algebra is complicated here - so I'll give you just an idea, but not a complete answer. I'll consider a circle with the center in D1.

Note, that distance(A,D1) = distance(B,D1) = distance(D1,T1). You can find a unit vector N, perpendicular to the line AB, and express vector D1 in terms of the midpoint vector M and the vector N, multiplied by a scalar parameter t:

D1 = M + t*N                                        (1)

The distance(C,T1) is known to be equal to cr, and also it's equal to the sum:

distance(C,T1) = distance(C,D1) + distance(D1,T1) = cr

This can be rewritten as:

distance(C,M + t*N) + distance(A,M + t*N) = cr      (2)

This is an equation with a single unknown variable t, which can be solved algebraically (it's not easy) - its solutions will give you two points D1 and D2.

Another option is to solve this equation numerically, incrementing the parameter t by small value and computing a result in a loop. This kind of solvers used to be available in various programming languages (not sure about Lua).

---

ADDITION #1. I used C++ and two open-source libraries to verify the idea above (I don't know Lua). Hopefully you'll be able to convert my code to Lua - it's elementary geometry plus finding a root of a function using bisection method. The code is below:

#include <cmath>
#include <iostream>

#include <boost/math/tools/roots.hpp>
#include <CGAL/Simple_cartesian.h>

namespace bmt = boost::math::tools;

using Kernel = CGAL::Simple_cartesian<double>;
using Circle = Kernel::Circle_2;
using Point = Kernel::Point_2;
using Vector = Kernel::Vector_2;

// ------ return vector, corresponding to the point P
Vector toVector(Point const& P)
{
  return P - CGAL::ORIGIN;
}

// ------ return unit vector, corresponding to the vector V 
Vector unitVector(Vector const& V)
{
  return V / std::sqrt(V.squared_length());
}

// ------ return the end point of the vector with origin in (0,0)
Point endPoint(Vector const& V)
{
  return CGAL::ORIGIN + V;
}

// ------ return the point of intersection of the circle C and its radius, passing through the point P 
Point boundaryPoint(Circle const& C, Point const& P)
{
  return C.center() + unitVector(P - C.center()) * std::sqrt(C.squared_radius());
}

// ------ given points and circle
Point const A{30, 30};
Point const B{40, 20};
Circle const C{{50, 50}, 40 * 40};

// ------ computed midpoint and perpendicular unit vector
auto const M = (toVector(A) + toVector(B)) / 2.0;
auto const N = unitVector((A - B).perpendicular(CGAL::CLOCKWISE)); 

// ------ the root of this function is used to find the center of the circle
double func(double const X)
{
  auto const D = endPoint(M + X * N);
  return
  (
    std::sqrt(CGAL::squared_distance(C.center(), D))
    +
    std::sqrt(CGAL::squared_distance(A, D))
    -
    std::sqrt(C.squared_radius())
  );
}

// ------ this function stops the root search when the root is found
bool tolerance(double const L, double const R)
{
  return (R - L) <= 0.00001;
}

// ------ the Boost Math Tools function 'bisect' is used to find the root in the interval [MIN,MAX] 
Point boundaryPoint(double const MIN, double const MAX)
{
  auto const t = bmt::bisect(func, MIN, MAX, tolerance);
  auto const D = endPoint(M + t.first * N);
  return boundaryPoint(C, D);
}

int main()
{
  auto const T1 = boundaryPoint(-50.0, 0.0);
  std::cout << T1.x() << '\t' << T1.y() << std::endl;
  auto const T2 = boundaryPoint(0.0, 50.0);
  std::cout << T2.x() << '\t' << T2.y() << std::endl;
}

---

Addition #2. This problem can be solved analytically, that is – without numerical methods. I’ll still use the parametric expression (1) for the point D1. The equation (2) can be elaborated a little bit further. Let’s find a point F, lying on the line, passing through points M and D1, and closest to the center C of the big circle. The vector C-F will be orthogonal to the unit vector N. Some notation:

a = half-length of the segment (A,B)
u = length of the segment (C,F)
v = length of the segment (M,F) 
R = radius of the big circle (= your cr) 

The segment (C,T1) is divided by the point D1 into two parts - (C,D1) and (D1,T1). The length of the segment (C,D1) is equal to sqrt(u^2 + (v-t)^2), and the length of the segment (D1,T1) is equal to the radius of the inscribed circle – that is sqrt(a^2 + t^2). So, the equation (2) becomes:

 sqrt(a^2 + t^2) + sqrt(u^2 + (v-t)^2) = R

We can try to solve this equation by hand, but it’s not easy. Fortunately, the symbolic algebra website "PolyMathLove.com" can solve it – please see here. However, I didn’t verify this solution.