Range Rings and Range Ellipses

[AGP]

Based on the discussion on range rings from this post https://sourceforge.net/p/geographiclib/discussion/1026621/thread/1df0a86b/ I have been able to recreate geodesic circles that are dead-on when compared to ArcGIS.

What I am researching now is how to extend that to "range ellipses". A description here http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/table-to-ellipse.htm. If we take this geodetic ellipse on the surface of the earth and generalize it to have a constant radius then we will have the range ring. The azimuth I think I can do by doing what I do for range circles but then I should use

For az = 0 To 360 Step N
azCalc = az + azTilt
Use azCalc and a length L to calculate end point
Next az

What I am not able to do is come up with an equation on how to compute the length L at an azCalc. Would it be based on simple ellipse formulas? or do geodesic calculations affect it? My instinct says yes the geodesics will come into play.

Any help or guidance is appreciated.

[Charles Karney]

For a plane ellipse, the polar representation relative to its center is

r(theta) = a*b / sqrt(b^2 * cos(theta)^2 + a^2 * sin(theta)^2)

see https://en.wikipedia.org/wiki/Ellipse#Polar_form_relative_to_center

Whether or not it is sensible to generalize this to a spheroid is up to you. But certainly this definition reduces to a geodesic circle when a = b.

My own preference would be to define the ellipse in terms of its foci: the sum of the distances to the two foci is constant for each point on the ellipse. In general, this will give you a different curve compared to the polar representation. However you will recover the geodesic circle when the two foci are coincident.

[AGP]

OK let me give this some thought and i will see what I can do with it.

[AGP]

Here is where I am with this challenge. I assume all my coordinates are in WGS84 and the input is ptStart and ptEnd. The I proceed as follows in .NET using the wrapper classes:

        'assume a WGS84 datum
        Dim geod As Geodesic = New Geodesic 'WGS84
        Dim s12 As Double, azi1 As Double, azi2 As Double
        Dim a12 As Double = geod.Inverse(ptStart.Y, ptStart.X, ptEnd.Y, ptEnd.X, s12, azi1, azi2)
        Dim geoLine As GeodesicLine = New GeodesicLine(geod, ptStart.Y, ptStart.X, azi1, Mask.ALL)

        'calculate geodesic distance half-way between the two points in meters. this will eventually be our semi-major axis.
        Dim halfDistMeters As Double = 0.5 * s12

        'compute the half-way point along the geodesic
        Dim xLonMid As Double, yLatMid As Double
        geoLine.Position(halfDistMeters, yLatMid, xLonMid)
        Dim ptMid As Point = New Point(xLonMid, yLatMid)

        'calculate the azimuth between the mid point to destination in degrees
        Dim azMidToPoa As Double = GeodesicAzimuthInDegrees(ptMid, ptEnd)

        'for great circle the flattening constant is zero. The flattening constant is f=(a-b)/a. if the semi-major axis, a, and the semi-minor axis, b,
        'are equal (as they would be in a perfect sphere) then that becomes f=(a-a)/a=0/a=0
        Dim a As Double = halfDistMeters
        Dim f As Double = 0.6 'large to test the flattening of the ellipse. change to desired factor.
        Dim b As Double = a * (1 - f)

        'we will hard code the angle sweep to be from 0 to 360 by 2 degrees. in the future we can probably make this user configurable.
        Dim radiusInMeters As Double, yLatPerim As Double, xLonPerim As Double
        For az As Double = 0 To 360 Step 2
            'radiusInMeters = GetDistAtTheta(az, a, b, azPodToPoa) 'this does not work
            radiusInMeters = GetDistAtTheta(az, a, b, azMidToPoa) 'this does work
            geod.Direct(ptMid.Y, ptMid.X, az, radiusInMeters, yLatPerim, xLonPerim)

            Dim ptPerim As Point = New Point(xLonPerim, yLatPerim)
        Next az

Public Function GeodesicAzimuthInDegrees(ByVal pt1 As Point, ByVal pt2 As Point) As Double
        'assume a WGS84 datum
        Dim geod As Geodesic = New Geodesic ' WGS84
        Dim s12gd As Double, azi1gd As Double, azi2gd As Double
        Dim a12gd As Double = geod.Inverse(pt1.Y, pt1.X, pt2.Y, pt2.X, s12gd, azi1gd, azi2gd)
        If azi1gd < 0 Then azi1gd += 360 'adjust for negative angle and always make positive
        Return azi1gd
End Function

Private Function GetDistAtTheta(ByVal azInDeg As Double, ByVal a As Double, ByVal b As Double, ByVal angRotationInDeg As Double) As Double
    Dim thetaInDeg As Double = AzimuthToPolar(azInDeg) - AzimuthToPolar(angRotationInDeg)

    'In polar coordinates, with the origin at the center of the ellipse and with the angular coordinate theta measured from the major axis,
    'the ellipse's equation is: r(theta) = a*b / sqrt(b^2 * cos(theta)^2 + a^2 * sin(theta)^2)
    Dim thetaInRad As Double = Math.PI * thetaInDeg / 180.0
    Dim r As Double = (a * b) / Math.Sqrt(b ^ 2 * Math.Cos(thetaInRad) ^ 2 + a ^ 2 * Math.Sin(thetaInRad) ^ 2)
    Return r
End Function

'https://math.stackexchange.com/questions/926226/conversion-from-azimuth-to-counterclockwise-angle
Private Function AzimuthToPolar(ByVal azInDeg As Double) As Double
    'For azimuth a, compute h=450−a; if h>360, return h−360; otherwise return h.
    Dim h As Double = 450 - azInDeg
    If h > 360 Then
        Return h - 360
    Else
        Return h
    End If
End Function

[AGP]

Here are a couple of samples:
MIA to ORD
BNA to CUN
LAX to YYZ

[Charles Karney]

OK, this looks fine. However, I remind you that you've only implemented one possible definition of an "ellipse" on a general surface. Which is the best definition to use will depend on your application.

[AGP]

Completely understand. I took the specific problem of a range circle and tried to generalize it to a range ellipse. I tried to incorporate geodesics to come up with a suitable definition. For my purposes I am trying to minimize the points used in a path finding algorithm. The flatness of the ellipse will determine how "far out" we will search for connection points. Many thanks for your help.