xx:=(d^2-rb^2+ra^2)/(2*d)

yy:=ra*Sqrt[1-(xx/ra)^2]

xx//OutputForm

Show[ParametricPlot[{ra*Cos[t],ra*Sin[t]}/.{ra->2},{t,0,2*π}],ParametricPlot[{rb*Cos[t]+d,rb*Sin[t]}/.{rb->4,d->5},{t,0,2*π}],ParametricPlot[{xx,t}/.{ra->2,rb->4,d->5},{t,-2,2}],ParametricPlot[{t,yy}/.{ra->2,rb->4,d->5},{t,-2,2}]]

vrs:={ra->3,rb->2.5,d->4}

Show[ParametricPlot[{ra*Cos[t],ra*Sin[t]}/.vrs,{t,0,2*π}],ParametricPlot[{rb*Cos[t]+d,rb*Sin[t]}/.vrs,{t,0,2*π}],ParametricPlot[{xx,t}/.vrs,{t,-2,2}],ParametricPlot[{t,yy}/.vrs,{t,-2,2}]]

Manipulate[Show[ParametricPlot[{ra*Cos[t],ra*Sin[t]}/.vrs,{t,0,2*π}],ParametricPlot[{yy*Cos[t]*Cos[θ]+xx*Sin[θ],yy*Sin[t]}/.vrs,{t,0,2*π}],PlotRange->{{0,4},{-4,4}}],{θ,0,π/2}]

soln=Solve[{ra*Cos[θ]+rbx*Cos[ϕ]==d,ra*Sin[θ]==rby*Sin[ϕ]},{θ,ϕ}]

qa=(Sqrt[(ra^4) (rbx^4)+(d^2) (ra^2) (rbx^2) (rby^2)-(ra^4) (rbx^2) (rby^2)-(ra^2) (rbx^4) (rby^2)+(ra^2) (rbx^2) (rby^4)])

qb=-d ra (rby^2)

qc=(ra^2) (rbx^2)-(ra^2) (rby^2)

qd=((2 d ra (rby^2) (Sqrt[(ra^2) (rbx^2) ((ra^2) (rbx^2)+(d^2) (rby^2)-(ra^2) (rby^2)-(rbx^2) (rby^2)+(rby^4))]))/((ra^2) (rbx^2)-(ra^2) (rby^2)))

qe=-(d^2) (rby^2)-(ra^2) (rby^2)+(rbx^2) (rby^2)-((2 (d^2) (ra^2) (rby^4))/((ra^2) (rbx^2)-(ra^2) (rby^2)))

qf=((ra (Sqrt[(ra^2) (rbx^2) ((ra^2) (rbx^2)+(d^2) (rby^2)-(ra^2) (rby^2)-(rbx^2) (rby^2)+(rby^4))]))/((ra^2) (rbx^2)-(ra^2) (rby^2)))

qg=d+((d (ra^2) (rby^2))/((ra^2) (rbx^2)-(ra^2) (rby^2)))

soln[[1,1,2,1,1]]==ArcTan[(qb-qa)/qc,-Sqrt[qe-qd]/Sqrt[qc]]

soln[[1,2,2,1,1]]==ArcTan[(qg+qf)/rbx,(ra/rby)*(-Sqrt[qe-qd]/Sqrt[qc])]

soln[[2,1,2,1,1]]==ArcTan[(qb-qa)/qc,Sqrt[qe-qd]/Sqrt[qc]]

soln[[2,2,2,1,1]]==ArcTan[(qg+qf)/rbx,(ra/rby)*(Sqrt[qe-qd]/Sqrt[qc])]

soln[[3,1,2,1,1]]==ArcTan[(qb+qa)/qc,-Sqrt[qe+qd]/Sqrt[qc]]

soln[[3,2,2,1,1]]==ArcTan[(qg-qf)/rbx,(ra/rby)*(-Sqrt[qe+qd]/Sqrt[qc])]

soln[[4,1,2,1,1]]==ArcTan[(qb+qa)/qc,Sqrt[qe+qd]/Sqrt[qc]]

soln[[4,2,2,1,1]]==ArcTan[(qg-qf)/rbx,(ra/rby)*(Sqrt[qe+qd]/Sqrt[qc])]

FullSimplify[qa]

FullSimplify[qb]

FullSimplify[qc]

FullSimplify[qd]

FullSimplify[qe]

FullSimplify[qf]

FullSimplify[qg]

FullSimplify[qa/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qb/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qc/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qd/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qe/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qf/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

FullSimplify[qg/.{rbx->iy*Cos[θ],rby->iy,d->ix*Sin[θ]}]

Solve[{ra*Cos[θ]+iy*Cos[p]*Cos[ϕ]==ix*Sin[p],ra*Sin[θ]==iy*Sin[ϕ]}/.{ix->xx,iy->yy},{θ,ϕ}]//FullSimplify

((2 d ra Sin[p])/((d^2)+(ra^2)-(rb^2)))*xx

xx

yy

psix[θ_]:=ArcSec[ra*Sin[θ]/xx]

alph[θ_]:=ArcTan[(ra*Cos[psix[θ]]-xx*Sin[θ])/(Cos[θ]),ra*Sin[psix[θ]]]

Manipulate[Show[ParametricPlot[{ra*Cos[t],ra*Sin[t]}/.vrs,{t,psix[θ]/.vrs,2*π-psix[θ]/.vrs}],ParametricPlot[{yy*Cos[t]*Cos[θ]+xx*Sin[θ],yy*Sin[t]}/.vrs,{t,2*π-alph[θ]/.vrs,2*π+alph[θ]/.vrs}],PlotRange->{{0,4},{-4,4}}],{θ,(ArcTan[xx/yy]+0.000000001)/.vrs,π/2}]