Navigation Mathematics

        C
       /|
      / |
   c /  | a
    /   |
   /    |
  /_____|
 A   b   B

Where:
- A, B, C are angles
- a, b, c are opposite sides
- C = 90° (right angle)

        North Pole
           /\
          /  \
         /    \
        /      \
    A  /________\ B
      \          /
       \   C    /
        \      /
         \    /
          \  /
           \/
         Equator

import math

# Convert to radians
lat1 = math.radians(33.9425)
lon1 = math.radians(-118.4081)
lat2 = math.radians(40.6413)
lon2 = math.radians(-73.7781)

dlat = lat2 - lat1
dlon = lon2 - lon1

# Haversine
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
R = 6371  # km

distance = R * c
print(f"Distance: {distance:.2f} km")  # 3974.21 km

      N
      │
  LAX → → → → → → → → JFK
  65° → 67° → 70° → ... → 90°

Initial bearing: ~66° (ENE)
Final bearing: ~90° (E)

haversine() {
    # Calculate great-circle distance using Haversine formula
    # Usage: haversine <lat1> <lon1> <lat2> <lon2> [unit]
    # Units: km (default), m, nm, mi
    # Example: haversine 33.9425 -118.4081 40.6413 -73.7781 km
    local lat1=$1 lon1=$2 lat2=$3 lon2=$4 unit=${5:-km}

    awk -v lat1="$lat1" -v lon1="$lon1" -v lat2="$lat2" -v lon2="$lon2" -v unit="$unit" '
    BEGIN {
        PI = 3.14159265358979323846

        # Convert to radians
        lat1_rad = lat1 * PI / 180
        lon1_rad = lon1 * PI / 180
        lat2_rad = lat2 * PI / 180
        lon2_rad = lon2 * PI / 180

        dlat = lat2_rad - lat1_rad
        dlon = lon2_rad - lon1_rad

        # Haversine formula
        a = sin(dlat/2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(dlon/2)^2
        c = 2 * atan2(sqrt(a), sqrt(1-a))

        # Select radius based on unit
        if (unit == "km") R = 6371.0088
        else if (unit == "m") R = 6371008.8
        else if (unit == "nm") R = 3440.065
        else if (unit == "mi") R = 3958.761
        else R = 6371.0088

        distance = R * c
        printf "%.4f\n", distance
    }'
}

    """
    Convert geodetic coordinates to ECEF (Earth-Centered, Earth-Fixed)

    Uses WGS84 ellipsoid parameters.

    Formula:
        N = a / sqrt(1 - e² * sin²(φ))
        X = (N + h) * cos(φ) * cos(λ)
        Y = (N + h) * cos(φ) * sin(λ)
        Z = (N(1 - e²) + h) * sin(φ)
    """
    lat_rad = math.radians(point.lat)
    lon_rad = math.radians(point.lon)

    sin_lat = math.sin(lat_rad)
    cos_lat = math.cos(lat_rad)
    sin_lon = math.sin(lon_rad)
    cos_lon = math.cos(lon_rad)

    # Radius of curvature in prime vertical
    N = WGS84_A / math.sqrt(1 - WGS84_E2 * sin_lat ** 2)

    x = (N + point.height) * cos_lat * cos_lon
    y = (N + point.height) * cos_lat * sin_lon
    z = (N * (1 - WGS84_E2) + point.height) * sin_lat

    return ECEF(x, y, z)


def ecef_to_geographic(ecef: ECEF) -> GeoPoint:
    """
    Convert ECEF to geodetic coordinates using Bowring's iterative method
    """
    x, y, z = ecef.x, ecef.y, ecef.z

    # Longitude is straightforward
    lon = math.atan2(y, x)

    # Iterative solution for latitude
    p = math.sqrt(x**2 + y**2)
    lat = math.atan2(z, p * (1 - WGS84_E2))  # Initial approximation

    for _ in range(10):  # Usually converges in 2-3 iterations
        sin_lat = math.sin(lat)
        N = WGS84_A / math.sqrt(1 - WGS84_E2 * sin_lat**2)
        lat_new = math.atan2(z + WGS84_E2 * N * sin_lat, p)