a <- (runif(2000) \* (40-10) + 10) \* pi/180


b <- (runif(2000) \* (40-10) + 10) \* pi/180


c <- (runif(2000) \* (40-10) + 10) \* pi/180


This chooses 2000 (pseudo-)random numbers in $(0,1)$ and then scatters them between 10 and 40, then multiplies by $\pi/180$, and does this three times independently, and calls the resulting vectors a, b, and c. Then


y <- (cos(a) + cos(b) + cos(c) + cos(a)\*cos(b)\*cos(c))/(sin(a) + sin(b) + sin(c) + sin(a)\*sin(b)\*sin(c))


z <- 1/tan((a+b+c)/3)


The above is as much of my code as I'll post at this late hour.


Plot z against y and it looks like a straight line! A trigonometric identity? No, it's not, as we shall see. Find the ordinary least-squares fit: $z=f + gy$. It turns out $f=0.2949023$ and $g = 0.7159759$ (results will vary if you use pseudorandomness rather than, as one might well do, chosing numbers systematically). Those numbers are weird.


Now look at the residuals from this least-squares fit. Plot $y$ on the horizontal axis and the residuals on the vertical axis. If this were really a trigonometric identity, you'd expect beensy residuals that differ from $0$ only because numerical computation is imperfect. But: The overwhelming majority of the residuals are between $\pm0.01$, most of the rest between $\pm0.02$, and only two of them, this time, bigger in absolute value than $0.03$. And the plot shows a pretty pattern, which is definitely neither a straight line nor a curve.


Since, it's *not* really a trigonometric identity, why do we need this microscope to see its deviation from straightness, and what's up with those funny coefficients, the slope and the $z$-intercept?



trigonometric non-identity
So sorry, I thought it was an identity but I typed lhs in place of rhs in one place. I probably wouldn't have answered otherwise. Now that I have, let me say that the linear fit does seem good. I can't say how unexpected that is, but perhaps not very. Here is a stab at an explanation along the lines of my previous silly answer. The short version is that the $y,z$ above, for the range given, are not *that* far off $\frac{4}{a+b+c}$ and $\frac{3}{a+b+c}$ so it is perhaps not amazing that $z \approx0.75y.$


You have $y=\frac{C(a)+C(b)+C(c)+C(a)C(b)C(c)}{S(a)+S(b)+S(c)+S(a)S(b)S(c)}$ and $z=\frac{C((a+b+c)/3)}{S((a+b+c)/3)}$ for $C(u)=\cos(u)$ and $S(u)=\sin(u)$ over the range $\frac{\pi}{18} \le u \le \frac{2\pi}{9} \lt 0.7.$


If we instead use the linear approximations $C(u)=1$ and $S(u)=u$ we are comparing $y=\frac{3+1}{a+b+c+abc}$ to $z=\frac{3}{a+b+c}.$ If you run your random experiments on those two I think that you will find the correlation coefficients similar and the residuals smaller. Those linear approximations are getting kind of off near the upper end of the range but at the middle $\theta=5\pi/36\approx 0.436,$ we have $\cos(\theta) \approx 0.906$ and $\sin(\theta) \approx 0.423.$ 


Keep that $S$ but use instead $C(u)=1-\frac{u^2}{2}$ to get $y=\frac{4-(a^2+b^2+c^2)+(a^2b^2+a^2c^2+b^2c^2)/4-a^2b^2c^2/8}{a+b+b+abc} \approx \frac{4-(a^2+b^2+c^2)}{a+b+b+abc}$ and $z=\frac{3-(a+b+c)^2/6}{a+b+c}.$ With those the correlation coefficients are again about the same and the residuals larger than in the linear case but smaller than for $\cos,\sin$


Since I worked it out: with $S(u)=u-u^3/6$, $C(u)=1-u^2/2$ 
$\frac{C(a)+C(b)+C(c)+C(a)C(b)C(c)}{S(a)+S(b)+S(c)+S(a)S(b)S(c)}=\frac{27(32-8(a^2+b^2+c^2)+2(a^2b^2+a^2c^2+b^2c^2)-a^2b^2c^2)}{216(a+b+c)+36(6abc-a^3-b^3-c^3)-36abc(a^2+b^2+c^2)+6(a^3+b^3+c^3)(abc)-a^3b^3c^3} $


and $\frac{C((a+b+c)/3)}{S((a+b+c)/3)}=\frac{9(18-(a+b+c)^2)}{(a+b+c)(54-(a+b+c)^2)}.$


I'll suggest that the higher order terms I dropped are not so significant both because the variables (with average value under $0.5$) occur to higher powers and because the coefficients are smaller. Higher order Taylor polynomials would have used $u^4/24$ and $u^5/120.$ 


I stuck to the form given although I might have preferred cross multiplying to have polynomials. It might be better to use polynomials centered at $u=\frac{5\pi}{36}$ which seem to look like $C(u) \approx 0.9986+0.0131u-0.5454u^2+0.0704u^4$ and $S(u)\approx -0.0005+1.0044u-0.0136u^2-0.1511u^3.$ That I leave to anyone who cares to check it out.


I'd also prefer to avoid the singularity at $0$ and work with the reciprocals. Then for $y=\frac{\sin a+\sin b +\sin c +\sin a\sin b \sin c}{\cos a+\cos b +\cos c +\cos a\cos b \cos c}$ and $z=\tan \frac{a+b+c}{3}$ the best linear fit for experiments like yours is $z \approx 1.021y+0.067.$


In that spirit I'll mentions that when $a=b=c=x$ we are comparing (the reciprocals) of $\tan(x)$ and $\frac{3\sin(x)+\sin^3(x)}{3\cos(x)+\cos^3(x)}=\tan(x)\frac{3+4\tan^2(x)}{4+3\tan^2(x)}.$ Those are equal at $x=0$ and $x=\frac{\pi}{4}.$ The derivatives differ more at those points than for the internal subinterval $(\frac{\pi}{18},\frac{2\pi}{9}).$ The exactly center interval of width $\frac{\pi}{6}$ is $(\frac{\pi}{24},\frac{5\pi}{24}).$ Perhaps the correlation would be even better, especially centered at $\frac{\pi}{8}$. This is only a downward shift of $\frac{\pi}{72}$ so it may not make much difference. Again, I don't plan to check.


