Integral of an angle that depends on a view The problem
There are two lines: $l_1$ and $l_2$. For any point $p$ in $\mathbb{R}^3$, we can calculate an angle between a plane that contains $p$ and $l_1$, and a plane that contains $p$ and $l_2$. Let's call this angle $A(p)$. Now what I need is to find a formula for $A(D)$, where $D$ is a planar figure contained by two lines and a
hyperbola. This needs to be implemented in a real-time computer application, so the formula has to be simple in terms of computational complexity. However, I haven't found any solution at all so far, so let's say the formula complexity is not that important for now.
What I have so far
We can always transform the system in a way, that the area $D$ is laying on the $XOY$ plane (so every point that belongs to $D$ will have coordinates
 $[x,y,0]$). Let's describe input lines $l_1$ and $l_2$ as two pairs of a point and a unit direction vector:
$$l_1: [a,b,c] + t[d,e,f]$$
$$l_2: [g,h,i] + u[j,k,l]$$
$$\text{where}: t, u \in \mathbb R$$
$A([x,y,z])$ equals to the angle between planes normals, and the normals can be found using cross product. If we normalize the vectors (dividing by their lengths), the angle can be retrieved by "reversing" the dot product with arcus cosinus:
$$\begin{align}A([x,y,0]) & = \operatorname{acos}\left(\frac{\left[x - a,y - b, -c\right] \times [d,e,f]}{\sqrt{(x - a)^2 + (y - b)^2 + c^2}} \circ \frac{[x - g, y - h, -i] \times [j, k, l]}{\sqrt{(x - g)^2 + (y - h)^2 + i^2}}\right) \\
& = \operatorname{acos} \left( \frac{\left((x - a)(x - g) + (y - b)(y - h) + ci \right) (dj + ek + fl) - \left((x - a)j + (y - b)k\right)\left((x - g)d + (y - h)e\right)}{\sqrt{(x - a)^2 + (y - b)^2 + c^2}\sqrt{(x - g)^2 + (y - h)^2 + i^2}} \right)\end{align}$$
To find the formula for $A(D)$ we simply need to integrate $A([x,y,0])$ over $D$:
$$\iint_D {A(x,y,0) \ \mathrm{d}x \ \mathrm{d}y}$$
I have no ideas of how to solve manually even the first, internal integral and the only CAS I have is Derive 6 - I had entered this integral and left computer running non-stop. After 4 days of
 nothing, I turned it off... And that's where I got so far.
I would appreciate any help on solving this problem.
 A: This is not a solution, but it might inspire someone (who's better in geometry) to find the final answer. Previously, I solved similar problem in $\mathbb{R}^2$:
I needed to find occludance at point E, which was equal to the length of an arc on the unit circle (with center in E), created by projecting occluder on
that circle. To calculate length of the arc, I needed to find angle alpha (arc length = arc radius * arc angle = 1 * alpha):
Click here for image (I have < 10 reputation points so I can't embed images yet...)
So I needed to find an angle beetween vectors $[P_x - E_x, P_y - E_y]$ and $[K_x - E_x, K_y - E_y]$. At first, I did that using the dot product:
$$A(x,0) = \operatorname{acos}\left(\frac{\overline{P} - \overline{E}}{|\overline{P} - \overline{E}|} \circ \frac{\overline{K} - \overline{E}} {|\overline{K} - \overline{E}|}\right)$$
Let's say that this angle formula is the final occludance at point E (one still needs to perform some scaling by $2\pi$ afterwards, but that's not the point here and let's just skip it).
Now what I needed was to find occludance not at a single point E, but for some whole segment of coordinates $[0;0]$ - $[L;0]$. So I needed to integrate $A(x)$ with lower bound limit = $0$ and the higher one = $L$. Let's give one letter labels to every vector coordinate ($m = P_x$; $n = P_y$; $k = K_x$; $l = K_y$) and after separating them, the result is:
$$\int_0^L{\operatorname{acos}\left(\frac{(x - k)(x - m) + ln}{\sqrt{(x-k)^2+l^2}\sqrt{(x-m)^2+n^2}}\right)dx}$$
It wasn't that hard to solve this integral, but the final computational complexity was to high in practice.
Then I found, that the whole problem can be solved without vector math, just bare trigonometry as soon as possible. Instead of taking $A(x)$ from the dot product, I took one angle that was between $Y$ axis and $\overline{EK}$, and second one that was between $Y$ axis and $\overline{EP}$, and then summed them up (that's possible because tan() is continuous in range $(-\frac{\pi}{2} ; \frac{\pi}{2})$), or rather calculate the difference between them because tan is $0$ on Y axis, $>0$ on the left side and $<0$ on the right side:
Click here for another image
$\operatorname{tan}(\operatorname{FirstAngle}) = \frac{k}{l}$, so $\operatorname{atan}(\frac{k}{l}) = \operatorname{FirstAngle}$. Analogicaly we can get the SecondAngle. Now our unit circle moves on $X$ axis from $0$ to $L$, so we need to adjust coordinates $k$ and $m$ by $L$. And that's how I got to a much better formula:
$$\int_0^L{\left(\operatorname{arctan}\left(\frac{k - x}{l}\right)-\operatorname{arctan}\left(\frac{m - x}{n}\right)\right)dx}$$
These two arcus tangens integrals are extremely simple to solve and their solutions are very easy in terms of comutational complexity.
Like I said in some comment above: having some kind of reference hyperplane (in this case $Y$ axis) and eliminating the square root from denominator were the "tricks" to success. I feel that similar "tricks" need to be apply in the 3D variant, but I can't find them.
A: I will refer constantly to the graphical representation at the bottom. It gathers the main ideas together with notations used hereafter.
One can see that reference axes are such that the $x$ axis is directed by the shortest distance segment and $y$ and $z$ axis have been chosen in such a way that the angles made by lines $L_1$ and $L_2$ are symmetrical with respect to plane $x0y$. 
Consider a generic point $M_0=(x_0,y_0,z_0)$. Let us call (P) the plane (in green on the figure) with equation $x=x_0$. In plane (P), "moving" from $Q_0(x_0,0,0)$ to $M_0(x_0,y_0,z_0)$ can be decomposed into a red "move" and a blue "move". In the red move, the first plane doesn't change ; in the blue move the second plane doesn't change. 
Let $c=\cos(\alpha)$ and $s=\sin(\alpha)$, with $u_1=(0,c,s)$ and $u_2=(0,-c,s)$ (unit vectors). 
Let $\vec{i}$ be the unit directing vector of $x$ axis.
An immediate computation gives 
$$\tag{1}c_1=\dfrac{z_0c+y_0s}{2sc} \ \ \text{and} \ \ c_2=\dfrac{z_0c-y_0s}{2sc}.$$
Thus, we can write:
$$\vec{Q_0P_0}=c_1\vec{u_2} \ \ \text{and} \ \ \vec{P_0M_0}=c_2\vec{u_1}.$$
Let us be more precise. The normal vector $N_2$ to the plane defined by $M_0$ and $L_2$ is:
$$N_2=\vec{u_2} \times \vec{A_2M_0} = \vec{u_2} \times (\vec{A_2Q_0}+\vec{Q_0P_0} + \vec{Q_0M_0})= $$
$$N_2=\vec{u_2} \times \vec{A_2Q_0}+ \vec{u_2} \times \vec{Q_0P_0} + \underbrace{ \vec{u_2} \times \vec{Q_0M_0}}_{0}$$
$$N_2=(a+x_0)\vec{u_2} \times \vec{i} +c_1\vec{u_2} \times \vec{u_1}$$
$$\tag{i}N_2=(x_0+a)\underbrace{\vec{u_2} \times \vec{i}}_{\text{unit vector}} -2c_1cs \ \vec{i}$$
For similar reasons, the normal vector $N_1$ to the plane defined by $M_0$ and $L_1$ is:
$$\tag{ii}N_1=(x_0-a) \vec{u_1} \times \vec{i} +2c_2cs \ \vec{i}$$
It suffices now to use formulas (i) and (ii) in order to have a vectorial expression of the looked for angle: 


*

*either as the arcsine of the cross product $N_1/\|N_1\| \times N_2/\|N_2\|.$ 

*or as the arccosine of the dot product $N_1/\|N_1\| \bullet N_2/\|N_2\|.$
Remark 1: One could arrive at formulas (i) and (ii) by using the following analytical computations:


*

*the plane defined by $(x_0,y_0,z_0)$ and line $L_1$ has the following equation:


$$\begin{vmatrix}a & a & x_0 & x\\
    0 & c & y_0 & y\\
    0 & s & z_0 & z\\
    1 & 1 & 1 & 1\end{vmatrix}=0 \iff acz_0 - acz + asy - asy_0 - cxz_0 + cx_0z + sxy_0 - sx_0y=0 $$


*

*the plane defined by $(x_0,y_0,z_0)$ and line $L_2$ has the following equation:


$$\begin{vmatrix}-a & -a & x_0 & x\\
    0 & -c & y_0 & y\\
    0 & s & z_0 & z\\
    1 & 1 & 1 & 1\end{vmatrix}=0 \iff 
acz_0 - acz - asy + asy_0 + cxz_0 - cx_0z + sxy_0 - sx_0y=0$$
From these equations, it is easy to extract the normal vectors to these two planes:
$N_1=\begin{pmatrix}sy_0-cz_0\\as-sx_0\\cx_0-ac\end{pmatrix}=(-a+x_0)\begin{pmatrix}0\\c\\s\end{pmatrix}\times\begin{pmatrix}1\\0\\0\end{pmatrix}+2c_1cs\begin{pmatrix}1\\0\\0\end{pmatrix}$ which, using the expression of $c_1$ given in (1) is identical to (ii) (up to a global sign change).
and $N_2=\begin{pmatrix}sy_0+cz_0\\-as-sx_0\\-cx_0-ac\end{pmatrix}=(-a-x_0)\begin{pmatrix}0\\-c\\s\end{pmatrix}\times\begin{pmatrix}1\\0\\0\end{pmatrix}+2c_2cs\begin{pmatrix}1\\0\\0\end{pmatrix}$ which is identical to (i) (up to a global sign change).
2) It is interesting to note that, in formulas (i) and (ii), $2cs=sin(2\alpha)$.

