Magnetic Navigation

                        Geographic North
                              │
                              │ D (Declination)
                              │╱
                    ─────────────────────  Horizontal plane
                             ╲│
                              │
                              │ I (Inclination/Dip)
                              ▼
                           F (Total)

     X = H cos(D)  (North component)
     Y = H sin(D)  (East component)
     Z            (Down component)

Your compass reads: 045° (magnetic)
Declination: 12° East

True Bearing = 045° + 12° = 057° True

If plotting on map (True → Magnetic):
Map bearing: 057° True
Compass setting: 057° - 12° = 045° Magnetic

Your compass reads: 045° (magnetic)
Declination: 16° West

True Bearing = 045° - 16° = 029° True

If plotting on map (True → Magnetic):
Map bearing: 029° True
Compass setting: 029° + 16° = 045° Magnetic

           20°W    15°W    10°W    5°W    0°    5°E    10°E   15°E   20°E
             │       │       │      │      │      │       │      │      │
   ┌─────────┼───────┼───────┼──────┼──────┼──────┼───────┼──────┼──────┤
   │         │       │       │      │      │      │       │      │      │
   │ Pacific │       │       │   Agonic    │      │       │      │      │
   │ Northwest│      │       │    Line     │      │       │      │      │
   │         │       │       │      │      │      │       │      │      │
   │      Seattle    │       │      │   Chicago  │       │  NYC │      │
   │         │       │       │      │      │      │       │ Boston     │
   │    San Francisco│       │  Denver     │      │       │      │      │
   │         │       │       │      │      │      │       │      │      │
   │    Los Angeles  │       │      │ Dallas     │  Atlanta      │      │
   │      (12°E)     │       │      │      │      │   (6°W)│      │      │
   │         │       │       │      │      │      │       │      │      │
   └─────────┴───────┴───────┴──────┴──────┴──────┴───────┴──────┴──────┘
                                    ↑
                              Agonic Line
                         (0° declination)

Year    Agonic Line Position (approximate US crossing)
1900    Through Maine, passing west of NYC
1950    Through Ohio, Great Lakes
2000    Through Michigan, Mississippi
2025    Through Louisiana, Wisconsin
2050    Expected through Texas, Minnesota (projected)

MAGNETIC NORTH DECLINATION AT CENTER OF SHEET
GN     MN
 ╲    ╱
  ╲  ╱
   ╲╱
   │
   │ True North

Magnetic declination: 12°30' East
Annual change: 0°6' West

For 2025 (if map dated 2020):
Years elapsed: 5
Change: 5 × 6' = 30' = 0.5° West
Current declination: 12.5° - 0.5° = 12.0° East

#!/usr/bin/env python3
"""
IGRF-14 Magnetic Field Calculator
Military-grade implementation using spherical harmonics
"""

import numpy as np
from scipy.special import lpmv  # Associated Legendre polynomials

def schmidt_quasi_normalize(n, m, P_nm):
    """
    Apply Schmidt quasi-normalization to Legendre polynomial.

    Schmidt normalization: P_n^m(cos θ) * S_n^m
    where S_n^m = sqrt((2 - δ_m0) * (n-m)! / (n+m)!)
    """
    from math import factorial
    if m == 0:
        return P_nm
    else:
        S = np.sqrt(2 * factorial(n - m) / factorial(n + m))
        return S * P_nm

def igrf_field(lat, lon, alt_km, year, coeffs):
    """
    Calculate magnetic field components using IGRF model.

    Parameters:
        lat: Geodetic latitude (degrees)
        lon: Longitude (degrees)
        alt_km: Altitude above WGS84 ellipsoid (km)
        year: Decimal year (e.g., 2025.5)
        coeffs: Dictionary of Gauss coefficients

    Returns:
        X, Y, Z: Magnetic field components (nT)
    """
    # Convert to geocentric coordinates
    a = 6378.137  # WGS84 semi-major axis (km)
    f = 1/298.257223563  # WGS84 flattening

    lat_rad = np.radians(lat)
    lon_rad = np.radians(lon)

    # Geocentric radius and co-latitude
    sin_lat = np.sin(lat_rad)
    cos_lat = np.cos(lat_rad)

    # Earth radius at latitude (simplified)
    r = a + alt_km
    theta = np.pi/2 - lat_rad  # Co-latitude

    # Reference radius for IGRF
    a_ref = 6371.2  # km

    # Initialize field components
    X, Y, Z = 0.0, 0.0, 0.0

    N_max = 13  # IGRF degree

    for n in range(1, N_max + 1):
        for m in range(0, n + 1):
            # Get coefficients (interpolate for year)
            g = coeffs.get((n, m, 'g'), 0)
            h = coeffs.get((n, m, 'h'), 0) if m > 0 else 0

            # Associated Legendre polynomial
            P = lpmv(m, n, np.cos(theta))
            P = schmidt_quasi_normalize(n, m, P)

            # Derivative of P with respect to theta
            if theta != 0 and theta != np.pi:
                dP = (n * np.cos(theta) * P -
                      (n + m) * lpmv(m, n-1, np.cos(theta))) / np.sin(theta)
            else:
                dP = 0

            # Radial scaling
            ratio = (a_ref / r) ** (n + 2)

            # Accumulate components
            cos_m_lon = np.cos(m * lon_rad)
            sin_m_lon = np.sin(m * lon_rad)

            X += ratio * (g * cos_m_lon + h * sin_m_lon) * dP
            Y += ratio * m * (g * sin_m_lon - h * cos_m_lon) * P / np.sin(theta) if theta != 0 else 0
            Z -= ratio * (n + 1) * (g * cos_m_lon + h * sin_m_lon) * P

    return X, Y, Z

def calculate_elements(X, Y, Z):
    """Calculate all seven magnetic elements from X, Y, Z."""
    H = np.sqrt(X**2 + Y**2)
    F = np.sqrt(X**2 + Y**2 + Z**2)
    D = np.degrees(np.arctan2(Y, X))  # Declination
    I = np.degrees(np.arctan2(Z, H))  # Inclination

    return {
        'X': X,      # nT
        'Y': Y,      # nT
        'Z': Z,      # nT
        'H': H,      # nT
        'F': F,      # nT
        'D': D,      # degrees
        'I': I       # degrees
    }

     160°W   120°W   80°W    40°W    0°     40°E    80°E   120°E   160°E
       │       │       │       │      │       │       │       │       │
   12°N┼───────┼───────┼───────┼──────┼───────┼───────┼───────┼───────┼
       │       │       │       │      │       │       │       │       │
       │       │       │    ╱──┼──────┼───────┼──╲    │       │       │
       │  ╱────┼───────┼───╱   │      │       │   ╲───┼───────┼────╲  │
       │ ╱     │       │  ╱    │      │       │    ╲  │       │     ╲ │
    0°N┼───────┼───────┼───────┼──────┼───────┼───────┼───────┼───────┼
       │ ╲     │       │  ╲    │      │       │    ╱  │       │     ╱ │
       │  ╲────┼───────┼───╲   │      │       │   ╱───┼───────┼────╱  │
   12°S┼───────┼───────┼───────┼──────┼───────┼───────┼───────┼───────┼
       │       │       │       │      │       │       │       │       │
                                       ↑
                           Magnetic Equator (I = 0°)
                      Maximum deviation: ~12° from geographic equator

     Direction of travel arrow
              ↑
    ┌─────────┴─────────┐
    │                   │
    │    ╭─────────╮    │
    │    │    N    │    │ ← Rotating bezel (declination adjustment)
    │    │  ╱   ╲  │    │
    │    │ W     E │    │ ← Magnetic needle (red end = North)
    │    │  ╲   ╱  │    │
    │    │    S    │    │
    │    ╰─────────╯    │
    │                   │ ← Baseplate (transparent for map work)
    │  ──────────────── │ ← Orienting lines
    │                   │
    └───────────────────┘

        Sighting wire
             │
    ┌────────┼────────┐
    │ Cover  │        │
    │  ╲     │     ╱  │
    │   ╲    │    ╱   │
    │    ╲   │   ╱    │
    │     ╲  │  ╱     │ ← Lens for reading dial
    │      ╲ │ ╱      │
    │       ╲│╱       │
    │        N        │
    │     W     E     │ ← Floating dial (mils or degrees)
    │        S        │
    └─────────────────┘

1. Locate declination adjustment screw (small screw on bezel)
2. Using screwdriver or key:
   - For EAST declination: Rotate orienting arrow EAST (clockwise)
   - For WEST declination: Rotate orienting arrow WEST (counter-clockwise)
3. Set to your local declination value
4. Now compass reads TRUE bearings directly

1. Read magnetic bearing from compass
2. Apply declination mentally:
   - East declination: ADD to get True
   - West declination: SUBTRACT to get True

Example (12°E declination):
Compass reads: 045° Magnetic
True bearing: 045° + 12° = 057° True

# Using geomag library
pip install geomag

python3 << 'EOF'
import geomag
from datetime import date

# McAllen, TX coordinates
lat, lon = 26.2034, -98.2300
today = date.today()

declination = geomag.declination(lat, lon, today)
print(f"Declination at McAllen, TX: {declination:.2f}°")
print(f"Direction: {'East' if declination > 0 else 'West'}")
EOF

source examples/navigation/nav-calc.sh

# Convert magnetic to true
apply_declination 45 12 E   # 57.00 (12°E declination)
apply_declination 45 16 W   # 29.00 (16°W declination)

    T   V   M   D   C
    ↓   ↓   ↓   ↓   ↓
  TRUE → MAGNETIC → COMPASS
       ↑ add West  ↑ add West
       ↓ subtract  ↓ subtract

Going DOWN (T to C): Add West, Subtract East
Going UP (C to T):   Subtract West, Add East

1. Position vessel at a known location with reference marks
2. Use pelorus or gyro compass as reference for TRUE headings
3. Steady on each cardinal and intercardinal heading (8 points minimum)
4. Record compass reading vs true heading at each point
5. Calculate deviation: δ = True - Compass - Declination
6. Construct deviation table/curve

Standard headings: 000°, 045°, 090°, 135°, 180°, 225°, 270°, 315°

1. Grid the area with 100m spacing
2. Record magnetic elements at each point
3. Compare with WMM/IGRF predictions
4. Map residuals (observed - predicted)
5. Contour anomaly regions
6. Mark anomaly zones on charts/maps

1. Plant a straight stick vertically
2. Mark shadow tip position at time T1
3. Wait 15-30 minutes
4. Mark shadow tip position at time T2
5. Draw line from T1 to T2
6. This line runs approximately East-West
7. Perpendicular to this line = North-South axis

Accuracy: ±5° (improves with longer wait time)
Note: More accurate near noon when shadow moves fastest

1. Point hour hand at the sun
2. Bisect angle between hour hand and 12 o'clock
3. This line points approximately South (before noon) or North (after noon)

Correction for Daylight Saving Time:
- Use 1 o'clock instead of 12 o'clock

Accuracy: ±15° (depends on latitude and time of year)