documentclass[crop,tikz,margin=10pt]{standalone}
usepackage{tikz}
usetikzlibrary{intersections}
usetikzlibrary{spy}
usepackage{pgfplots}
begin{document}
begin{tikzpicture}[
    spy using outlines = {circle,size=3cm,magnification=5,connect spies},
]
begin{axis}
addplot+[domain = 0:2*pi] expression {sin(deg(x))};
coordinate (spy point) at (axis cs: 0, 0);
coordinate (magnifying glass) at (rel axis cs: -0.4, 0.2);
end{axis}

spy on (spy point) in node at (magnifying glass);
end{tikzpicture}
begin{tikzpicture}
defradiusa{0.3}
defradiusb{3}

defxa{5}
defya{3}

defxb{0}
defyb{0}

coordinate (magnifying glass) at (xa, ya);
coordinate (spy point) at (xb, yb);

pgfmathsetmacroxp{(xb * radiusa - xa * radiusb) / (radiusa - radiusb)}
pgfmathsetmacroyp{(yb * radiusa - ya * radiusb) / (radiusa - radiusb)}
pgfmathsetmacrodistancea{sqrt((xp - xa) * (xp - xa) + (yp - ya) * (yp - ya) - radiusa * radiusa))}
pgfmathsetmacrodistanceb{sqrt((xp - xb) * (xp - xb) + (yp - yb) * (yp - yb) - radiusb * radiusb))}
pgfmathsetmacrodenoma{(xp - xa)*(xp - xa) + (yp - ya)*(yp - ya)}
pgfmathsetmacrodenomb{(xp - xb)*(xp - xb) + (yp - yb)*(yp - yb)}

pgfmathsetmacroxc{(radiusa * radiusa * (xp - xa) + radiusa * (yp - ya) * distancea) / denoma + xa}
pgfmathsetmacroyc{(radiusa * radiusa * (yp - ya) - radiusa * (xp - xa) * distancea) / denoma + ya}

pgfmathsetmacroxe{(radiusa * radiusa * (xp - xa) - radiusa * (yp - ya) * distancea) / denoma + xa}
pgfmathsetmacroye{(radiusa * radiusa * (yp - ya) + radiusa * (xp - xa) * distancea) / denoma + ya}

pgfmathsetmacroxd{(radiusb * radiusb * (xp - xb) + radiusb * (yp - yb) * distanceb) / denomb + xb}
pgfmathsetmacroyd{(radiusb * radiusb * (yp - yb) - radiusb * (xp - xb) * distanceb) / denomb + yb}

pgfmathsetmacroxf{(radiusb * radiusb * (xp - xb) - radiusb * (yp - yb) * distanceb) / denomb + xb}
pgfmathsetmacroyf{(radiusb * radiusb * (yp - yb) + radiusb * (xp - xb) * distanceb) / denomb + yb}


draw (magnifying glass) circle(radiusa);
draw (spy point) circle(radiusb);

% draw (xa, ya) node[scale=3, green] {.};
% draw (xb, yb) node[scale=3, green] {.};
% draw (xp, yp) node[scale=3, blue] {.};
% draw (xc, yc) node[scale=3, red] {.};
% draw (xd, yd) node[scale=3, red] {.};
% draw (xe, ye) node[scale=3, red] {.};
% draw (xf, yf) node[scale=3, red] {.};

% draw (xa, ya) -- (xp, yp);
% draw (xb, yb) -- (xp, yp);

draw (xc, yc) -- (xd, yd);
draw (xe, ye) -- (xf, yf);
end{tikzpicture}
end{document}

documentclass[crop,tikz,margin=10pt]{standalone}
usepackage{tikz}
usetikzlibrary{intersections}
usetikzlibrary{spy}
usepackage{pgfplots}
makeatletter
newcommandxcoord[2][center]{{%
    pgfpointanchor{#2}{#1}%
    pgfmathparse{pgf@x/pgf@xx}%
    pgfmathprintnumber{pgfmathresult}%
}}
newcommandycoord[2][center]{{%
    pgfpointanchor{#2}{#1}%
    pgfmathparse{pgf@y/pgf@yy}%
    pgfmathprintnumber{pgfmathresult}%
}}
makeatother
begin{document}
begin{tikzpicture}
begin{scope}[
    spy using outlines = {circle,size=3cm,magnification=5},
]
begin{axis}
addplot+[domain = 0:2*pi] expression {sin(deg(x))};
coordinate (spy point) at (axis cs: 0, 0);
coordinate (magnifying glass) at (rel axis cs: -0.4, 0.2);
coordinate (a) at (rel axis cs: 0.5, 1.1);
end{axis}
node at (a) {spy point: xcoord{spy point}, ycoord{spy point}, glass: xcoord{magnifying glass}, ycoord{magnifying glass}};

spy on (spy point) in node at (magnifying glass);
end{scope}
coordinate (c) at (-1.6589690159337693, 0.10010961563788023);
coordinate (d) at (0.7862061968132461, 2.6420219231275763);
coordinate (e) at (-2.962541913303562, 2.6233998438799935);
coordinate (f) at (0.5254916173392875, 3.1466799687759988);
draw (c) -- (d);
draw (e) -- (f);
end{tikzpicture}
end{document}

from math import sqrt
import matplotlib.pyplot as plt


def main():
    radiusa = 1.5
    radiusb = 1.5 / 5
    xa = -2.74
    ya = 1.14
    xb = 0.57
    yb = 2.85

    figure, ax = plt.subplots()
    circlea = plt.Circle((xa, ya), radiusa, color='C0')
    circleb = plt.Circle((xb, yb), radiusb, color='C0')
    ax.add_artist(circlea)
    ax.add_artist(circleb)

    (xc, yc), (xd, yd), (xe, ye), (xf, yf) = compute(xa, ya, radiusa, xb, yb, radiusb)

    ax.plot([xc, xd], [yc, yd], color='C0')
    ax.plot([xe, xf], [ye, yf], color='C0')

    ax.set_xlim(
        min(xa - radiusa, xb - radiusb),
        max(xa + radiusa, xb + radiusb),
    )
    ax.set_ylim(
        min(ya - radiusa, yb - radiusb),
        max(ya + radiusa, yb + radiusb),
    )

    print("\coordinate (c) at ({}, {});".format(xc, yc))
    print("\coordinate (d) at ({}, {});".format(xd, yd))
    print("\coordinate (e) at ({}, {});".format(xe, ye))
    print("\coordinate (f) at ({}, {});".format(xf, yf))

    plt.show()


def compute(xa, ya, radiusa, xb, yb, radiusb):
    xp = (xb * radiusa - xa * radiusb) / (radiusa - radiusb)
    yp = (yb * radiusa - ya * radiusb) / (radiusa - radiusb)
    distancea = sqrt((xp - xa) * (xp - xa) + (yp - ya) * (yp - ya) - radiusa * radiusa)
    distanceb = sqrt((xp - xb) * (xp - xb) + (yp - yb) * (yp - yb) - radiusb * radiusb)
    denoma = (xp - xa)*(xp - xa) + (yp - ya)*(yp - ya)
    denomb = (xp - xb)*(xp - xb) + (yp - yb)*(yp - yb)

    xc = (radiusa * radiusa * (xp - xa) + radiusa * (yp - ya) * distancea) / denoma + xa
    yc = (radiusa * radiusa * (yp - ya) - radiusa * (xp - xa) * distancea) / denoma + ya

    xe = (radiusa * radiusa * (xp - xa) - radiusa * (yp - ya) * distancea) / denoma + xa
    ye = (radiusa * radiusa * (yp - ya) + radiusa * (xp - xa) * distancea) / denoma + ya

    xd = (radiusb * radiusb * (xp - xb) + radiusb * (yp - yb) * distanceb) / denomb + xb
    yd = (radiusb * radiusb * (yp - yb) - radiusb * (xp - xb) * distanceb) / denomb + yb

    xf = (radiusb * radiusb * (xp - xb) - radiusb * (yp - yb) * distanceb) / denomb + xb
    yf = (radiusb * radiusb * (yp - yb) + radiusb * (xp - xb) * distanceb) / denomb + yb

    return (xc, yc), (xd, yd), (xe, ye), (xf, yf)


if __name__ == '__main__':
    main()