A hard geometry problem involving harmonic divisions 
Let acute triangle $ABC$. Let $A_1$ and $A_2$ the intersections of the circle of diameter $(BC)$ and the altitude from $A$ to $BC$ ($A_1$ is closer to $A$ than $A_2$). Similarily define points $B_1$, $B_2$, $C_1$, $C_2$. Let $A'$ the intersection of $B_1C_2$ and $B_2C_1$. Similarily define points $B'$ and $C'$. Prove that $AA'$, $BB'$ and $CC'$ concur.


My idea: Maybe it helps to see that $(A,H;A_1,A_2)$ (where $H$ is the orthocenter of triangle $ABC$) is a harmonic division. Also $H$ is the radical center of the 3 drawn circles, so $B_1C_1B_2C_2$ is cylcic (by power of a point). It is obvious by the figure that $A'$ must lie on $BC$. But I don't know how to prove it. Can you please help me? Thanks in advance!
If it helps, the problem comes from a Romanian book about harmonic divisions, but the solution is ommited.
 A: We solve the problem using trilinear coordinates. The altitude $AD$ is the set of points whose coordinates $x:y:z$ satisfy
$$y\cos B=z\cos C$$
The circle with diameter $BC$ is defined analogously, the points satisfying
$$yz=x(x\cos A-y\cos B-z\cos C)$$
(See exercise (174) here for the reference.) Arbitrarily setting $x=1$ (since trilinear coordinates are ratios) and then solving for $y,z$ gives us the coordinates of $A_1$ and $A_2$:
$$A_{1,2}=1:
-\cos C\pm\sqrt{\frac{\cos C}{\cos B}(\cos A+\cos B\cos C)}:
-\cos B\pm\sqrt{\frac{\cos B}{\cos C}(\cos A+\cos B\cos C)}$$
The plus sign gives $A_1$ and the minus sign gives $A_2$; $B_1,B_2,C_1,C_2$ can be obtained by cyclically permuting $A,B,C$ in the above equation.
Now associate the vector $(u,v,w)^T$ with both the point at coordinates $u:v:w$ and the line $ux+vy+wz=0$. It is well-known that the line through points $P_1$ and $P_2$ is $(\mathbf P_1×\mathbf P_2)\cdot(x,y,z)^T=0$ and that the intersection of lines $l_1$ and $l_2$ is $\mathbf l_1×\mathbf l_2$. Based on this, the intersection of the lines $B_1C_2$ and $C_1B_2$ is
$$A'=(\mathbf B_1×\mathbf C_2)×(\mathbf C_1×\mathbf B_2)$$
$$=0:(\cos A\cos C+\cos B)\sqrt{\cos C(\cos A\cos B+\cos C)}:(\cos A\cos B+\cos C)\sqrt{\cos B(\cos A\cos C+\cos B)}$$
Thus $A'$ lies on $BC$ as you suspected. The line $AA'$ then has normal vector $\mathbf l_A=\mathbf A'×(1,0,0)^T$, and similarly for $\mathbf l_B=BB'$ and $\mathbf l_C=CC'$ by cyclically permuting $A,B,C$; the determinant of the matrix formed by these three vectors is
$$\begin{vmatrix}\mathbf l_A&\mathbf l_B&\mathbf l_C\end{vmatrix}=0$$
Hence the lines $AA',BB',CC'$ concur, as was required to be shown, at the point with trilinear coordinates
$$X=\sqrt{1+\frac{\cos B\cos C}{\cos A}}:\sqrt{1+\frac{\cos C\cos A}{\cos B}}:\sqrt{1+\frac{\cos A\cos B}{\cos C}}$$
$$=\frac1{a\sqrt{b^2+c^2-a^2}}:\frac1{b\sqrt{c^2+a^2-b^2}}:\frac1{c\sqrt{a^2+b^2-c^2}}$$
$$=\frac1{\sqrt{a\cos A}}:\frac1{\sqrt{b\cos B}}:\frac1{\sqrt{c\cos C}}$$

Here is the SymPy code I used to derive all the expressions above:
#!/usr/bin/env python3
from sympy import *
cA, cB, cC = symbols('cA cB cC', positive=True) # cos A, cos B, cos C
x, y, z = symbols('x y z', real=True)

def cycB(p): # ABC -> BCA
    q = p.subs({cA: cB, cB: cC, cC: cA}, simultaneous=True)
    return Matrix([q[2], q[0], q[1]])
def cycC(p): # ABC -> CAB
    q = p.subs({cA: cC, cB: cA, cC: cB}, simultaneous=True)
    return Matrix([q[1], q[2], q[0]])

f1 = y*cB - z*cC
f2 = cA - y*cB - z*cC - y*z
sols = solve([f1, f2], [y, z])
A1 = Matrix([1, sols[1][0].expand(), sols[1][1].expand()])
A2 = Matrix([1, sols[0][0].expand(), sols[0][1].expand()])
print("A1 =", A1)
print("A2 =", A2)
B1 = cycB(A1)
B2 = cycB(A2)
C1 = cycC(A1)
C2 = cycC(A2)
Ap = simplify(  B1.cross(C2).cross(B2.cross(C1))  ) # A'
Ap *= sqrt(cA*cB*cC)/2
print("A' =", Ap)
lA = Ap.cross(Matrix([1, 0, 0]))
lB = cycB(lA)
lC = cycC(lA)
D = Matrix([lA.T, lB.T, lC.T])
pprint(D)
print("det(D) =", D.det()) # 0

X = D.nullspace()[0] * sqrt(cA*cB + cC) / sqrt(cC)
a, b, c = symbols('a b c', positive=True)
X = X.subs(cA, (b**2+c**2-a**2)/(2*b*c))
X = X.subs(cB, (c**2+a**2-b**2)/(2*c*a))
X = X.subs(cC, (a**2+b**2-c**2)/(2*a*b))
Delta = sqrt(-(a - b - c)*(a - b + c)*(a + b - c))*sqrt(a + b + c)/sqrt(2) # area of triangle
X = factor(X, deep=True) / Delta
print("X =", X.simplify())

