Formula for distance between Incenter and Orthocenter.

Here is a try to get the formula quickly.

Disclaimer: I will use more or less standard notations in a triangle $\Delta ABC$. The side lengths are $a,b,c$, the incenter is $I$, the inradius is $r$ (not $\rho$), the circumradius is $R$, the area is $S$. Let $AA'$, $B'$, $CC'$ be the heights intersecting in $H$. The inradius of the $H$-pedal triangle $\Delta A'B'C'$ will be denoted by $r_H$.

The barycentric coordinates of $I$ are $I=[a:b:c]$, explicitly we have (identifying with affixes) $$ (a+b+c)I=aA+bB+cC\ . $$ Then considering vectors: $$ \begin{aligned} a(a+b+c)\overrightarrow{IA} &= a^2\overrightarrow{AA} + ab\overrightarrow{BA} + ac\overrightarrow{CA} \ ,\\ b(a+b+c)\overrightarrow{IB} &= ba\overrightarrow{AB} + b^2\overrightarrow{BB} + bc\overrightarrow{CB} \ ,\\ c(a+b+c)\overrightarrow{IC} &= ca\overrightarrow{AC} + cb\overrightarrow{BC} + c^2\overrightarrow{CC} \ , \end{aligned} $$ so after adding: $$ \tag{$1$} a\overrightarrow{IA}+ b\overrightarrow{IB}+ c\overrightarrow{IC} =0\ . $$

Write now $\overrightarrow{HA}=\overrightarrow{HI}+\overrightarrow{IA}$, take norms, and multiply with $a$. Write the similar equalities for $B$, $C$ instead of $A$ to get: $$ \begin{aligned} aHA^2 &=aHI^2 + aIA^2 + 2a\overrightarrow{HI}\cdot \overrightarrow{IA}\ ,\\ bHB^2 &=bHI^2 + bIB^2 + 2b\overrightarrow{HI}\cdot \overrightarrow{IB}\ ,\\ cHC^2 &=cHI^2 + cIC^2 + 2c\overrightarrow{HI}\cdot \overrightarrow{IC}\ . \end{aligned} $$ Adding, we get rid of the scalar products by $(1)$, and obtain: $$ \tag{$2$} aHA^2+bHB^2+cHC^2 = (a+b+c)HI^2 + \underbrace{aIA^2+bIB^2+cIC^2}_{\text{some constant }K^3}\ . $$

Note: We did not use any special property of $H$, so $(2)$ is valid for any point $X$ instead of $H$, explicitly $$ \tag{$2_X$} aXA^2+bXB^2+cXC^2 = (a+b+c)XI^2 + K^3\ . $$ (Key-word: Leibnitz formula.)

(For $I$ we have used only $I=[a:b:c]$. So we can also replace $I$ by a general point $P=[u:v:w]$, using the corresponding barycentric coordinates $u,v,w$ instead.)

Now compute the relation $a(2_A)+b(2_B)+c(2_C)$ divided by $(a+b+c)$, weights being used to make appear on the R.H.S. the term $\sum a AI^2=K^3$. So after the sum we have on the R.H.S. $K^3+K^3$. On the L.H.S. we get $\frac 1{a+b+c}\sum(abc^2+acb^2)=2abc$. This gives the value for $K^3$, $$ K^3 = abc = 4R\;S=4R\; rp=2Rr\;(a+b+c)\ . $$ Here $p=\frac 12(a+b+c)$ is the half perimeter. Inserting this into $(2)$ we get $$ \tag{$3$} \frac1{a+b+c}(aHA^2+bHB^2+cHC^2) = HI^2 + \underbrace{2Rr}_{abc/(a+b+c)}\ . $$

We want now to show $(!)$ the following similar relation: $$ \tag{$4$} 2r^2 \overset != HI^2 + 2Rr_H\ . $$ Showing it is equivalent to: $$ \tag{$4'$} \frac1{2p}(aHA^2+bHB^2+cHC^2) - \frac{abc}{2p} \overset != 2r^2 - 2Rr_H\ . $$

It is time to write down some trigonometric formulas. $$ \begin{aligned} BC &= a\ ,\\ B'C' &= BC\cos A = a\cos A &&\text{ from }\Delta ABC\sim \Delta AC'B'\ ,\\ AH &= AB'\sin\widehat{AHB'}=c\cos A\sin C= 2R\cos A\ ,\\ HA'&=HB\cos\widehat{AHB'} =2R\cos B\cos C\ ,\\ r_H &=HA'\cdot \frac{B'C'}{BC}=2R\cos A\cos B\cos C &&\text{ from }\Delta HB'C'\sim \Delta HCB\ ,\text{ so}\\ 2Rr_H &= HA\cdot HA'=HB\cdot HB'=HC\cdot HC'\ . \end{aligned} $$

I understand that trigonometry should be avoided, but then computations are hard to write. Above, there are also some alternative relations showing how to avoid - if really needed - trigonometric formulas. At any rate, no trigonometric relations will be used below. (E.g. like writing $\cos^2A$ in terms of $\cos 2A$ or so.)

We want to show $(4')$ and are in position to show this as an algebraic formula in the variables $a,b,c$. For this, we use for $S^2$ Heron, and $R$ is eliminated via $abc=4RS$. The version of $(4')$ suited for this polynomial check is: $$ \tag{$4''$} \frac1{2p}4R^2(a\cos^2A + b\cos^2B + c\cos^2 C) - \frac{abc}{2p} \overset != \frac{2S^2}{p^2} - 4R^2\cos A\cos B\cos C\ . $$ Equivalently: $$ \tag{$4'''$} \frac1{2p}4R^2(a\cos^2A + b\cos^2B + c\cos^2 C) + \frac1{2p}4R^2(2p\; \cos A\cos B\cos C) \overset != \frac{abc+4(p-a)(p-b)(p-c)}{2p} \ . $$ And we multiply with $2p$.

I have to finish, so a computer check is done for the above.

var('a,b,c')
p = (a + b + c)/2
SS = p*(p - a)*(p - b)*(p - c)
RR = (a*b*c/4)^2 / SS
def f(a, b, c):    return (b^2 + c^2 - a^2)/2/b/c
cosA, cosB, cosC = f(a,b,c), f(b,c,a), f(c,a,b)

RHS = a*b*c + 4*(p - a)*(p - b)*(p - c)
LHS = 4*RR*(a*cosA^2 + b*cosB^2 + c*cosC^2 
            + 2*p*cosA*cosB*cosC)

print(bool(LHS == RHS))

And the computer checks the relation. Explicitly:

sage: print( bool(LHS == RHS) )
True
sage: 2*RHS.expand()
-a^3 + a^2*b + a*b^2 - b^3 + a^2*c + b^2*c + a*c^2 + b*c^2 - c^3
sage: 2*LHS.factor().expand()
-a^3 + a^2*b + a*b^2 - b^3 + a^2*c + b^2*c + a*c^2 + b*c^2 - c^3

(I will search for a better finish based on geometric relations.)