Solid Angle of Ellipse¶

1. Problem Description¶

Calculate the (signed) solid angle for an ellipse from the origin \((0,0,0)\)

In short, complete the function s.t.

/**
 * @brief Computes the solid angle for an ellipse from the origin(0,0,0)
 *
 * @param c Center of the ellipse. R^3
 * @param n Normal vector of the ellipse. R^3
 * @param d Direction of the long axis. R^3  ensure that <n, d> = 0
 * @param a Semi-major axis length of the ellipse. R
 * @param b Semi-minor axis length of the ellipse. R ensure that a >= b >= 0
 * @return The solid angle in steradians.
 */
double solid_angle_ellipse(
    const Eigen::Vector3d &c,
    const Eigen::Vector3d &n,
    const Eigen::Vector3d &d,
    const double a,
    const double b
);

2. Math Expression Needed¶

2.1 Problem Simplify¶

We can assume that the Normal vector of the ellipse \(\boldsymbol{n}\) is parallel to the z-axis and the direction of the long axis \(\boldsymbol{d}\) is parallel to the x-axis, because it's easy to apply a coordinate transformation to satisfy this condition.

So, in the following expression, we assume that \(\boldsymbol{n} = (0,0,1)\) and \(\boldsymbol{d} = (1,0,0)\)

Also, we assume that \(\boldsymbol{c} = (p,q,h)\)

2.2 Analytical solution for the coaxial case¶

It's easy to calculate the solid angle using the eletric field. As we suppose a point charge is placed at the origin, and the soliet angle problem can be transformed into a calculation of the electric flux.

So the answer is

where

and

where \(k' \equiv \sqrt{1-k^2}\)

2.3 Analytical case related to sections of a right circular cone¶

We can apply a coordiante transformation such that, in the transformerd frame, the \(z'\)-axis is aligned with the direction from the ellipse center to the observation point, and the \(x'\) axis is aligned with the ellipse's major axis.

After this transformation, the resulting elliptical cone can be characterized by its new effective semi-axes, which can then be used to compute the solid angle.

And here's the new axis length after transformation.

and then, we can apply the conclusion of Chapter 2.2 to calculate the answer, but here's the point, in this occasion, we're not sure that if \(a > b\) or \(a < b\), so, we should first swap \(a\) and \(b\) if \(a < b\) before apply the conclusion.

2.4 General Solution¶

2.4.1 Rotation of the coordiante system¶

In the \((x,y,z)\) coordinate system, the sheared core obeys the equation

we need to rotate the coordinate system \((x,y,z)\) to \((x',y',z')\) where the \(z'\) axis lies along the center line of the right elliptic cone.

This rotation need a right-handed rotation by an angle \(\gamma\) about the \(y\)-axis to give a new coordinate system \((\hat{x},\hat{y},\hat{z})\) and a additional right-handed rotation by an angle \(\psi\) about the \(\hat{x}\) axis to give the new coordinate system \((x',y',z')\)

Then, we get

Substituting Eq.(119) into Eq.(117), expanding and collecting terms gives

where

and where

2.4.2 Solution for the rotation angle¶

For the \(z'\) axis to lie along the axis of the cone, the coefficients \(D_1\) and \(E_1\) must vanish simultaneously and hence

After solving the system of equations, we get

where \(w_j\), \(j \in {1,2,3}\) is the three roots of the function

it is proved from the geometry that all three roots are real number.

only one root that satisfied \(B_1^2-4A_1C_1 < 0\) is the correct root that needed.

2.4.3 Equation of the ellipse in canonical form¶

The value of \(z'\) in Eq.(127) for \(F_1\) is arbitrary and can be conveniently set to unity.

Since \(E_1\) and \(F_1\) is zero now, we can apply another right-handed rotation by an angle \(\lambda\) about the \(z'\) axis to set \(B_1\) to zero.

We denoting the new rotated coordinates by \((\bar{x},\bar{y},\bar{z})\), then do the same steps of the previous part.

Then, we can get the final equation of the ellipse in the form

where

and

as the equation of Eq.(57), we can get the final solid angle which is equal to

3. Results¶

3.1 time Complexity¶

The horizontal axis repersents the number of triangles used to divide the ellipse in the numerical intergration, while the vertical axis shows the average time required to compute the solid angle. The orange curve correspoinds to the numerical-integration approach, whereas the red dashed line denotes the analytical formula(\(1.9123\mu s\)).

3.2 computational error¶

The error was determined by comparing our results with those reported in the paper, which are accurate to 20 significant digits The horizontal axis repersents the number of triangles used to divide the ellipse in the numerical intergration, while the vertical axis shows the average error in the calculated solid angle. The orange curve correspoinds to the numerical-integration approach, whereas the red dahed line denotes the analytical formula(\(9.20974e-15\)).

summery table

3.3 accelerate¶

if we disable the elliptic integral promotion to long double, we can get the table following

if we disable the elliptic integral promotion to double, we can get the table following. But for some extreme samples, analytical fomula error would goto 1e-5, which will maybe not be accpeted.

评论