Finding the width of the envelope around an uncertainty ellipse

Written by Vince Cronin in Mathematica 8, starting on August 16, 2011;  revised 23 August 2011

Problem

We are given the vertical uncertainty (a) and the horizontal uncertainty (b) in km where, in general, a ≥ b.  Uncertainties a and b are along the major and minor axes of the uncertainty ellipse, respectively.  We are also given the slope of the nodal plane (δ) in degrees.

Imagine a line that is inclined at an angle δ relative to horizontal and that is tangent to the uncertainty ellipse.

Question:  What is the shortest distance between the tangent line and the center of the ellipse?

Approach:  Find the point (D) at which the line is tangent to the ellipse.  If the center of the ellipse is labeled point O, let line segment DO be the hypotenuse of a right triangle.  The second side of the triangle through vertex D is coincident with the tangent line.  The third side of the triangle, from a point P along the tangent line to point O, forms a right angle with the tangent line.  The length of side PO is the half - width of the envelope around the uncertainty ellipse.

Input

The user must supply input data in the colored boxes below.  Variable a is the vertical uncertainty of the focal location, in kilometers.

a=3;

Variable b is the horizontal uncertainty of the focal location, in kilometers.

b=2;

Variable deltaIn is the dip angle of the nodal plane, in degrees.

deltaIn=60;

Description of solution

Some useful relationships

Variable a is the length of the semi - major axis of the uncertainty ellipse, and b is the length of the semi - minor axis of the ellipse.  The foci of the ellipse are a distance of 2 c apart along the major axis, where

c = (Sqrt[(a^2)-(b^2)])

c=(Sqrt[(a^2)- (b^2)]);

The ellipse is the set of all points (x, y) for which

((x^2)/(b^2))+((y^2)/(a^2))=1

The parametric equations of the ellipse are

x = b Cos[θ]  and  y = a Sin[θ],

where θ is measured anticlockwise from the positive x axis;  that is, θ is measured from horizontal.

Recall the Law of Sines

(a/(Sin[A]))=(b/(Sin[B]))=(c/(Sin[C]))

for a triangle with side lengths a, b and c and opposite interior angles of A, B and C, respectively.

Computation for the typical case in which the vertical uncertainty is greater than or equal to the horizontal uncertainty

Refer to the illustration below for a graphical definition of the point and angle labels.

The input dip angle was expressed in degrees.  Convert degrees to radian measure.

δ=deltaIn ((π/180));

We know the length of two sides of the triangle f-g-h:  side fg = 2c and side gh = 2 a.  The angle between the lines fh and fg is π-δ and is known from input data.  From the Law of Sines, we can determine the angle α formed between lines gh and fh.

((2 a)/(Sin[π-δ]))=((2 c)/(Sin[α]))

α=ArcSin[((2 c Sin[π-δ])/(2 a))];

The length of side fh is computed using the Law of Sines

((2 a)/(Sin[π-δ]))=(fh/(Sin[δ-α]))

fh=((2 a Sin[δ-α])/(Sin[π-δ]));

Now focusing on triangle d - f - h, we know that this is an isoceles triangle (length dh = length df).   The angle α (computed previously) is between lines dh and fh, and also between lines hf and df.  We can compute the length of side df using the Law of Sines.

(fh/(Sin[π-(2 α)]))=(df/(Sin[α]))

df=((fh Sin[α])/(Sin[π-(2 α) ]));

The three interior angles of triangle d - f - g are :  2 α between sides gd and fd, (π-(δ+α)) between sides df and gf, and (π-(ϕ+(π-(δ+α)))) between sides fg and dg.  More simply put, the interior angles are 2 α, (π - δ - α), and (δ - α).

(df/(Sin[δ-α]))=(dg/(Sin[π-δ-α]))

dg=((df Sin[π-δ-α])/(Sin[δ-α]));

The distance from the tangent point (d) to the y axis is called xd

xd=dg Sin[δ-α];

The distance from the tangent point (d) to the x axis is called yd

yd=(dg Cos[δ-α])-c;

We shall call the center of the ellipse point o, as in “origin of the coordinate system.”  The distance from point o to the tangent point d is called dLength

dLength=(Sqrt[(xd^2)+ (yd^2)]);

We also label one last point (p) at the intersection of the tangent line through point d with the line through o that is perpendicular to the tangent line.  We find the half width of the uncertainty envelope by analyzing the right triangle with apices at points o, d and p.

The angle between (1) the line segment between point o and point d, and (2) the horizontal (x) axis is given by

ArcTan[(yd/xd)]

The angle between the tangent line and the horizontal (x) axis is the same as the dip angle of the nodal plane:  δ.  The angle between op and the horizontal (x) axis is

((π/2)-δ)

The angle between sides od and op is given by

ArcTan[(yd/xd)] - ((π/2)-δ)

The distance between the tangent line and point o is tha half - width of the uncertainty envelope.

envelopeHalfWidth=dLength Cos[ArcTan[(yd/xd)]-((π/2)-δ)];

N[%]

The old estimate involved the determination of the half width of the envelope around a cylindrically shaped uncertainty volume (Cronin, V.S., Millard, M., Seidman, L, and Bayliss, B., 2008, The Seismo-Lineament Analysis Method [SLAM] -- A Reconnaissance Tool to Help Find Seismogenic Faults: Environmental and Engineering Geoscience, v. 14, no. 3, p. 199-219).

oldEstimate=(a Cos[δ])+(b Sin[δ]);

N[%]

Module

The following module embeds the code described above.  Three values must be passed to the module:  in1 is the vertical uncertainty in km, which must be greater than or equal to the horizontal uncertainty;  in2 is the horizontal uncertainty in km; and in3 is the dip angle of the nodal plane in degrees.  The module returns the half width of the uncertainty volume, measured between parallel enveloping surfaces that are tangent to opposite sides of the uncertainty ellipsoid and that have the same dip angle as the nodal plane.

findHalfWidth[in1_,in2_,in3_]:=Module[{m1,m2,m3,m4,m5,m6,m7,m8,m9,result},m1=(Sqrt[(in1^2)-(in2^2)]);m2=in3 ((π/180));m3=ArcSin[((2 m1 Sin[π-m2])/(2 a))];m4=((2 a Sin[m2-m3])/(Sin[π-m2]));m5=((m4 Sin[m3])/(Sin[π-(2 m3)]));m6=((m5 Sin[π-m2-m3])/(Sin[m2-m3]));m7=m6 Sin[m2-m3];m8=(m6 Cos[m2-m3])-m1;m9=(Sqrt[(m7^2)+ (m8^2)]);result=m9 Cos[ArcTan[(m8/m7)]-((π/2)-m2)];result];

Output

The answer is given in kilometer units.

answer=N[findHalfWidth[a,b,deltaIn]]

References

The Ellipse, accessed 16 August 2011 online via http://home.scarlet.be/~ping1339/ellipse.htm