Consider a point $P~(P_x,P_y,0)$ which lies somewhere in the Cartesian region $x > 0. $
Consider a simple 3D cone surface with half angle $\gamma$, originating at point $\mathcal{O}$ $(0,0,0)$ whose axis is the x-axis $(x,0,0)$.
The cone may be considered as consisting of a bundle of straight line segments originating at point $\mathcal{O}$ and each making an angle $\gamma$ with the x-axis.
I wish to "reflect" the point P perpendicularly across each of these lines to produce a set of points. These points can be connected back to origin point $\mathcal{O}$ to produce a second set of lines which describe a cone-like surface.
I require a formula to calculate the "average" vector V obtained by averaging the x,y,z components of the second set of lines.
I have used basic (high school) trigonometric formulae and simple calculus to derive an expression for the average vector as follows:-
$V_x = P_x (1 - 2 sin^2(\gamma$))
$V_y = - P_y (1 - cos^2(\gamma)sin(\gamma$))
$V_z = 0$.
I have tested the above expression with a simple computer model.
The formula for $V_x$ and $V_z$ seem to be accurate to 1 in 30 trillion, which is nice. But the formula for $V_y$ seems to be in error as accurate to only 1 in 4000. When I try the alternative formula:
$V_y = - P_y (1 - cos^2(\gamma)sin^2(\gamma))$
This agrees with the computer model to within 1 in 40 trillion, which is nice.
But I cannot presently see how to derive the adjusted formula for $V_y$.
I wish to know if this sort of problem is addressed elsewhere and if there is a derivation available.
Tl;dr: my answer is
$$\begin{split} V_x &= P_x(1 - 2\sin^2\gamma),\\ V_y &= -P_y\cos^2\gamma,\\ V_z &= 0. \end{split}$$
Let's take just one line and reflect your point across it. Line is described by its vector $\textbf{a}$ so reflected point $Q$ is given by this formula:
$$\textbf{q}=2\frac{\textbf{p}\cdot\textbf{a}}{|\textbf{a}|^2}\textbf{a} - \textbf{p},$$ where $\textbf{p}$ and $\textbf{q}$ are just position vectors of points $P$ and $Q$ respectively.
Next, let's parametrize $\textbf{a}$. We can do it as follows: $\textbf{a}=(\cos\gamma, \sin\alpha\,\sin\gamma, \cos\alpha\,\sin\gamma)$. It can easily be seen that when $\alpha$ goes from $0$ to $2\pi$, the end of vector $\textbf{a}$ draws exactly a circle on your cone surface. Moreover, $|\textbf{a}|=1$.
Considering that $\textbf{p}\cdot\textbf{a} = P_x\cos\gamma + P_y\sin\alpha\,\sin\gamma$, we have
$$\textbf{q}=\begin{bmatrix} 2(P_x\cos\gamma + P_y\sin\alpha\,\sin\gamma)\cos\gamma -P_x\\ 2(P_x\cos\gamma + P_y\sin\alpha\,\sin\gamma)\sin\alpha\,\sin\gamma -P_y\\ 2(P_x\cos\gamma + P_y\sin\alpha\,\sin\gamma)\cos\alpha\,\sin\gamma \end{bmatrix}.$$
Let's prepare it for the final step:
$$\begin{align} \textbf{q}_x &= \sin\alpha\cdot 2P_y\sin\gamma\,\cos\gamma + P_x(2\cos^2\gamma - 1),\\ \textbf{q}_y &= \frac{1-\cos 2\alpha}{2}\cdot 2P_y\sin^2\gamma + \sin\alpha\cdot2 P_x\cos\gamma\,\sin\gamma - P_y,\\ \textbf{q}_z &= \frac{1}{2}\sin2\alpha\cdot 2P_y\sin^2\gamma + \cos\alpha\cdot 2P_x\sin\gamma\,\cos\gamma. \end{align}$$
Last thing to note is that averaging any of $\sin\alpha,\cos\alpha,\sin2\alpha,\cos2\alpha$ over $[0, 2\pi]$ gives you $0$ because it's the multiple of a period for each of those functions. So
$$\begin{split} V_x &= P_x(2\cos^2\gamma - 1) = P_x(1 - 2\sin^2\gamma),\\ V_y &= P_y\sin^2\gamma - P_y = -P_y\cos^2\gamma,\\ V_z &= 0. \end{split}$$