15284544330

Committed 27 May 2025 07:45PM UTC coverage: 47.35%. First build

Build # 15284544330

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #51 from victoriyaforsythe/release_candidate

Release candidate

Run Details

398 of 572 new or added lines in 8 files covered. (69.58%)

956 of 2019 relevant lines covered (47.35%)

3.31 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

75.13

/PyIRI/igrf_library.py

#!/usr/bin/env python
# ----------------------------------------------------------
# Distribution statement A. Approved for public release.
# Distribution is unlimited.
# This work was supported by the Office of Naval Research.
# ----------------------------------------------------------
"""This library contains components for IGRF.

References
----------
Alken et al. (2021). International geomagnetic reference field: the
thirteenth generation. Earth, Planets and Space, 73(49),
doi:10.1186/s40623-020-01288-x.

"""

import numpy as np
import os
from scipy import interpolate


def verify_inclination(inc):
    """Test the validity of an inclination input.

    Parameters
    ----------
    inc : int, float, or array-like
        Magnetic inclination in degrees.

    Returns
    -------
    inc : int, float, or array-like
        Magnetic inclination in degrees

    Raises
    ------
    ValueError
        If inclination is not between -90 and 90 degrees

    """
    # Ensure the inclination to be tested is a numpy array
    test_inc = np.asarray(inc)
    if test_inc.shape == ():
        test_inc = np.asarray([inc])

    # Get a mask that specifies if each inclination is realistic
    good_inc = (test_inc <= 90.0) & (test_inc >= -90.0)

    # Return original input if all are good, correct or raise error otherwise
    if good_inc.all():
        return inc
    else:
        # Adjust values within a small tolerance
        test_inc[(test_inc > 90.0) & (test_inc <= 90.1)] = 90.0
        test_inc[(test_inc < -90.0) & (test_inc >= -90.1)] = -90.0
        good_inc = (test_inc <= 90.0) & (test_inc >= -90.0)

        if good_inc.all():
            # Return the adjusted input using the original type
            if isinstance(inc, (float, int, np.floating, np.integer)):
                return test_inc[0].astype(type(inc))
            elif isinstance(inc, list):
                return list(test_inc)
            else:
                return test_inc
        else:
            raise ValueError('unrealistic magnetic inclination value(s)')


def inc2modip(inc, alat):
    """Calculate modified dip angle from magnetic inclination.

    Parameters
    ----------
    inc : array-like
        Magnetic inclination in degrees.
    alat : array-like
        Flattened array of latitudes in degrees.

    Returns
    -------
    modip_deg : array-like
        Modified dip angle in degrees.

    Notes
    -----
    This function calculates modified dip angle for a given inclination and
    geographic latitude.

    """
    inc = verify_inclination(inc)

    alat_rad = np.deg2rad(alat)
    rad_arg = np.deg2rad((inc / np.sqrt(np.cos(alat_rad))))
    modip_rad = np.arctan(rad_arg)
    modip_deg = np.rad2deg(modip_rad)

    return modip_deg


def inc2magnetic_dip_latitude(inc):
    """Calculate magnetic dip latitude from magnetic inclination.

    Parameters
    ----------
    inc : array-like
        Magnetic inclination in degrees.

    Returns
    -------
    magnetic_dip_latitude : array-like
        Magnetic dip latitude in degrees.

    Notes
    -----
    This function calculates magnetic dip latitude.

    """
    inc = verify_inclination(inc)
    arg = 0.5 * (np.tan(np.deg2rad(inc)))
    magnetic_dip_latitude = np.rad2deg(np.arctan(arg))

    return magnetic_dip_latitude


def inclination(coeff_dir, date_decimal, alon, alat, alt=300., only_inc=True):
    """Calculate magnetic inclination using IGRF13.

    Parameters
    ----------
    coeff_dir : str
        Where IGRF13 coefficients are located.
    date_decimal : float
        Decimal year
    alon : array-like
        Flattened array of geographic longitudes in degrees.
    alat : array-like
        Flattened array of geographic latitudes in degrees.
    alt : flt
        Altitude in km, default is 300 km.
    only_inc : bool
        Only output the inclination, otherwise output the magnetic field
        information as well (default=True)

    Returns
    -------
    inc : array-like
        Magnetic inclination in degrees.
    X : array-like (optional)
        North component of the magnetic field in nT.
    Y : array-like (optional)
        East component of the magnetic field in nT.
    Z : array-like (optional)
        Vertical component of the magnetic field in nT.
    dec : array-like (optional)
        Declination of the magnetic field in degrees.
    hoz : array-like (optional)
        Horizontal intensity of the magnetic field in nT.
    eff : array-like (optional)
        Total intensity of the magnetic field in nT.

    Raises
    ------
    IOError
        If unable to find IGRF coefficient file

    Notes
    -----
    This code is a slight modification of the IGRF13 pyIGRF release,
    https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
    The reading of the IGRF coefficient file was modified to speded up the
    process, and the main code was simlified.

    """
    # Set altitude
    aalt = np.full(shape=alon.shape, fill_value=alt)

    # Open IGRF file and read to array the main table
    igrf_file = os.path.join(coeff_dir, 'IGRF', 'IGRF13.shc')
    if not os.path.isfile(igrf_file):
        raise IOError('unable to find IGRF coefficient file: {:}'.format(
            igrf_file))

    with open(igrf_file, mode='r') as fopen:
        file_array = np.genfromtxt(fopen, delimiter='', skip_header=5)

    # Create the time array of years that match the years in the file
    # (check if file is change to the newer version)
    igrf_time = np.arange(1900, 2025 + 5, 5)

    # Exclude first 2 columns, these are the m and n indecies
    igrf_coeffs = file_array[:, 2:]

    # Maximum degree of polynomials
    nmax = 13

    # Colatitude (from 0 to 180)
    colat = 90.0 - alat

    # Compute geocentric colatitude and radius from geodetic colatitude
    # and height
    alt, colat, sd, cd = gg_to_geo(aalt, colat)

    # Interpolate coefficients from 5 year slots to given decimal year
    igrf_extrap = interpolate.interp1d(igrf_time, igrf_coeffs,
                                       fill_value='extrapolate')
    coeffs = igrf_extrap(date_decimal)

    # Compute the main field B_r, B_theta and B_phi value for the location(s)
    br, bt, bp = synth_values(coeffs.T, alt, colat, alon, nmax)

    # Rearrange to X, Y, Z components
    X = -bt
    Y = bp
    Z = -br

    # Rotate back to geodetic coords if needed
    t = X
    X = X * cd + Z * sd
    Z = Z * cd - t * sd

    # Compute the four non-linear components
    dec, hoz, inc, eff = xyz2dhif(X, Y, Z)

    if only_inc:
        out = inc
    else:
        out = (inc, X, Y, Z, dec, hoz, eff)

    return out


def gg_to_geo(h, gdcolat):
    """Compute geocentric colatitude and radius from geodetic colat and height.

    Parameters
    ----------
    h : array-like
        Altitude in kilometers.
    gdcolat : array-like
        Geodetic colatitude in degrees.

    Returns
    -------
    rad : array-like
        Geocentric radius in kilometers.
    thc : array-like
        Geocentric colatitude in degrees.
    sd : array-like
        Rotate B_X to gd_lat.
    cd : array-like
        Rotate B_Z to gd_lat.

    Notes
    -----
    IGRF-13 Alken et al., 2021.

    References
    ----------
    Equations (51)-(53) from "The main field" (chapter 4) by Langel,
    R. A. in: "Geomagnetism", Volume 1, Jacobs, J. A., Academic Press,
    1987.

    Malin, S.R.C. and Barraclough, D.R., 1981. An algorithm for
    synthesizing the geomagnetic field. Computers & Geosciences, 7(4),
    pp. 401-405.

    """
    # Use WGS-84 ellipsoid parameters
    eqrad = 6378.137  # equatorial radius
    flat = 1. / 298.257223563
    plrad = eqrad * (1. - flat)  # polar radius
    ctgd = np.cos(np.deg2rad(gdcolat))
    stgd = np.sin(np.deg2rad(gdcolat))
    a2 = eqrad * eqrad
    a4 = a2 * a2
    b2 = plrad * plrad
    b4 = b2 * b2
    c2 = ctgd * ctgd
    s2 = 1. - c2
    rho = np.sqrt(a2 * s2 + b2 * c2)
    rad = np.sqrt(h * (h + 2. * rho) + (a4 * s2 + b4 * c2) / rho**2)
    cd = (h + rho) / rad
    sd = (a2 - b2) * ctgd * stgd / (rho * rad)
    cthc = ctgd * cd - stgd * sd

    # Clip cthc to [-1, 1] not to cause problems in arccos
    cthc = np.clip(cthc, -1.0, 1.0)

    # Arccos returns values in [0, pi] so we need to convert to degrees
    thc = np.rad2deg(np.arccos(cthc))  # arccos returns values in [0, pi]

    return rad, thc, sd, cd


def geo_to_gg(radius, theta):
    """Compute geodetic colatitude and vertical height above the ellipsoid.

    Parameters
    ----------
    radius : array-like
        Geocentric radius in kilometers.
    theta : array-like
        Geocentric colatitude in degrees.

    Returns
    -------
    height : array-like
        Altitude in kilometers.
    beta : array-like
        Geodetic colatitude.

    Notes
    -----
    IGRF-13 Alken et al., 2021.
    Compute geodetic colatitude and vertical height above the ellipsoid from
    geocentric radius and colatitude. Round-off errors might lead to a
    failure of the algorithm especially but not exclusively for points close
    to the geographic poles. Corresponding geodetic coordinates are returned
    as NaN.

    References
    ----------
    Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
    geodetic coordinates", IEEE Transactions on Aerospace and Electronic
    Systems}, 1994, vol. 30, num. 3, pp. 957-961.

    """
    # Use WGS-84 ellipsoid parameters
    a = 6378.137  # equatorial radius
    b = 6356.752  # polar radius
    a2 = a**2
    b2 = b**2
    e2 = (a2 - b2) / a2  # squared eccentricity
    e4 = e2 * e2
    ep2 = (a2 - b2) / b2  # squared primed eccentricity
    r = radius * np.sin(np.radians(theta))
    z = radius * np.cos(np.radians(theta))
    r2 = r**2
    z2 = z**2
    F = 54. * b2 * z2
    G = r2 + (1. - e2) * z2 - e2 * (a2 - b2)
    c = e4 * F * r2 / G**3
    s = (1. + c + np.sqrt(c**2 + 2 * c))**(1. / 3.)
    P = F / (3 * (s + 1. / s + 1.)**2 * G**2)
    Q = np.sqrt(1. + 2 * e4 * P)
    r0 = (-P * e2 * r / (1. + Q)
          + np.sqrt(0.5 * a2 * (1. + 1. / Q) - P * (1. - e2) * z2
                    / (Q * (1. + Q)) - 0.5 * P * r2))
    U = np.sqrt((r - e2 * r0)**2 + z2)
    V = np.sqrt((r - e2 * r0)**2 + (1. - e2) * z2)
    z0 = b2 * z / (a * V)
    height = U * (1. - b2 / (a * V))
    beta = 90. - np.degrees(np.arctan2(z + ep2 * z0, r))

    return height, beta


def synth_values(coeffs, radius, theta, phi, nmax=None, nmin=1, grid=False):
    """Compute radial, colatitude and azimuthal field components.

    Parameters
    ----------
    coeffs : array-like
        Coefficients of the spherical harmonic expansion. The last
        dimension is equal to the number of coefficients, `N` at the grid
        points.
    radius : array-like
        Array containing the radius in km.
    theta : array-like
        Array containing the colatitude in degrees.
    phi : array-like
        Array containing the longitude in degrees.
    nmax : int or NoneType
        Maximum degree up to which expansion is to be used, if None it will be
        specified by the ``coeffs`` variable.  However, smaller values are
        also valid if specified. (default=None)
    nmin : int
        Minimum degree from which expansion is to be used. Note that it will
        just skip the degrees smaller than ``nmin``, the whole sequence of
        coefficients 1 through ``nmax`` must still be given in ``coeffs``.
        (default=1)
    grid : bool
        If ``True``, field components are computed on a regular grid. Arrays
        ``theta`` and ``phi`` must have one dimension less than the output
        grid since the grid will be created as their outer product
        (default=False).

    Returns
    -------
    B_radius : array-like
        Radial field components in km.
    B_theta : array-like
        Colatitude field components in degrees.
    B_phi : array-like
        Azimuthal field components in degrees.

    Raises
    ------
    ValueError
        If an inappropriate input value is supplied

    Notes
    -----
    by IGRF-13 Alken et al., 2021
    Based on chaosmagpy from Clemens Kloss (DTU Space, Copenhagen)
    Computes radial, colatitude and azimuthal field components from the
    magnetic potential field in terms of spherical harmonic coefficients.
    A reduced version of the DTU synth_values chaosmagpy code.

    References
    ----------
    Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
    geodetic coordinates", IEEE Transactions on Aerospace and Electronic
    Systems}, 1994, vol. 30, num. 3, pp. 957-961.

    """
    # Ensure ndarray inputs
    coeffs = np.array(coeffs, dtype=float)
    radius = np.array(radius, dtype=float) / 6371.2  # Earth's average radius
    theta = np.array(theta, dtype=float)
    phi = np.array(phi, dtype=float)
    if (np.amin(theta) < 0.0) | (np.amax(theta) > 180.0):
        raise ValueError('Colatitude outside bounds [0, 180].')

    if nmin < 1:
        raise ValueError('Only positive nmin allowed.')

    # Handle optional argument: nmax
    nmax_coeffs = int(np.sqrt(coeffs.shape[-1] + 1) - 1)  # degrees
    if nmax is None:
        nmax = nmax_coeffs
    elif nmax <= 0:
        raise ValueError('Only positive, non-zero nmax allowed.')

    if nmax > nmax_coeffs:
        nmax = nmax_coeffs

    if nmax < nmin:
        raise ValueError(f'Nothing to compute: nmax < nmin ({nmax} < {nmin}.)')

    # Manually broadcast input grid on surface
    if grid:
        theta = theta[..., None]  # First dimension is theta
        phi = phi[None, ...]  # Second dimension is phi

    # Get shape of broadcasted result
    try:
        broad = np.broadcast(radius, theta, phi,
                             np.broadcast_to(0, coeffs.shape[:-1]))
    except ValueError:
        raise ValueError(''.join(['Cannot broadcast grid shapes (excl. last ',
                                  'dimension of coeffs):\nradius: ',
                                  repr(radius.shape), '\ntheta: ',
                                  repr(theta.shape), '\nphi: ', repr(phi.shape),
                                  '\ncoeffs: ', repr(coeffs.shape)]))
    grid_shape = broad.shape

    # Initialize radial dependence given the source
    r_n = radius**(-(nmin + 2))

    # Compute associated Legendre polynomials as (n, m, theta-points)-array
    Pnm = legendre_poly(nmax, theta)

    # Save sinth for fast access
    sinth = Pnm[1, 1]

    # Calculate cos(m*phi) and sin(m*phi) as (m, phi-points)-array
    phi = np.radians(phi)
    cos_mp = np.cos(np.multiply.outer(np.arange(nmax + 1), phi))
    sin_mp = np.sin(np.multiply.outer(np.arange(nmax + 1), phi))

    # allocate arrays in memory
    B_radius = np.zeros(grid_shape)
    B_theta = np.zeros(grid_shape)
    B_phi = np.zeros(grid_shape)
    num = nmin**2 - 1
    for n in range(nmin, nmax + 1):
        B_radius += (n + 1) * Pnm[n, 0] * r_n * coeffs[..., num]
        B_theta += -Pnm[0, n + 1] * r_n * coeffs[..., num]
        num += 1
        for m in range(1, n + 1):
            B_radius += ((n + 1) * Pnm[n, m] * r_n
                         * (coeffs[..., num] * cos_mp[m]
                            + coeffs[..., num + 1] * sin_mp[m]))
            B_theta += (-Pnm[m, n + 1] * r_n
                        * (coeffs[..., num] * cos_mp[m]
                           + coeffs[..., num + 1] * sin_mp[m]))
            with np.errstate(divide='ignore', invalid='ignore'):
                # Handle poles using L'Hopital's rule
                div_Pnm = np.where(theta == 0., Pnm[m, n + 1],
                                   Pnm[n, m] / sinth)
                div_Pnm = np.where(theta == np.degrees(np.pi),
                                   -Pnm[m, n + 1], div_Pnm)
            B_phi += (m * div_Pnm * r_n
                      * (coeffs[..., num] * sin_mp[m]
                         - coeffs[..., num + 1] * cos_mp[m]))
            num += 2
        r_n = r_n / radius  # Equivalent to r_n = radius**(-(n+2))

    return B_radius, B_theta, B_phi


def legendre_poly(nmax, theta):
    """Calculate associated Legendre polynomials P(n,m).

    Parameters
    ----------
    nmax : int
        Maximum degree up to which expansion, with a minimum value of one.
    theta : array-like
        Array containing the colatitude in degrees.

    Returns
    -------
    Pnm : array-like
        Evaluated values and derivatives.

    Raises
    ------
    IndexError
        If nmax is less than 1

    Notes
    -----
    by IGRF-13 Alken et al., 2021
    Returns associated Legendre polynomials P(n,m) (Schmidt
    quasi-normalized), and the derivative :d P(n,m) / d theta evaluated
    at theta.

    References
    ----------
    Langel, R. A. "Chapter four: Main field." Geomagnetism, edited by JA Jacobs
    1 (1987).
    Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
    geodetic coordinates", IEEE Transactions on Aerospace and Electronic
    Systems}, 1994, vol. 30, num. 3, pp. 957-961.

    """
    costh = np.cos(np.radians(theta))
    sinth = np.sqrt(1 - costh**2)
    Pnm = np.zeros((nmax + 1, nmax + 2) + costh.shape)
    Pnm[0, 0] = 1  # is copied into trailing dimenions
    Pnm[1, 1] = sinth  # write theta into trailing dimenions via broadcasting
    rootn = np.sqrt(np.arange(2 * nmax**2 + 1))

    # Recursion relations after Langel "The Main Field" (1987),
    # eq. (27) and Table 2 (p. 256)
    for m in range(nmax):
        Pnm_tmp = rootn[m + m + 1] * Pnm[m, m]
        Pnm[m + 1, m] = costh * Pnm_tmp
        if m > 0:
            Pnm[m + 1, m + 1] = sinth * Pnm_tmp / rootn[m + m + 2]
        for n in np.arange(m + 2, nmax + 1):
            d = n * n - m * m
            e = n + n - 1
            Pnm[n, m] = ((e * costh * Pnm[n - 1, m]
                          - rootn[d - e] * Pnm[n - 2, m]) / rootn[d])

    # dP(n,m) = Pnm(m,n+1) is the derivative of P(n,m) vrt. theta
    Pnm[0, 2] = -Pnm[1, 1]
    Pnm[1, 2] = Pnm[1, 0]
    for n in range(2, nmax + 1):
        n = int(n)
        Pnm[0, n + 1] = - np.sqrt((n * n + n) / 2.) * Pnm[n, 1]
        Pnm[1, n + 1] = ((np.sqrt(2.0 * (n * n + n)) * Pnm[n, 0]
                          - np.sqrt((n * n + n - 2)) * Pnm[n, 2]) / 2.0)
        for m in np.arange(2, n):
            m = int(m)
            Pnm_part1 = np.sqrt((n + m) * (n - m + 1)) * Pnm[n, m - 1]
            Pnm_part2 = np.sqrt((n + m + 1.0) * (n - m)) * Pnm[n, m + 1]
            Pnm[m, n + 1] = 0.5 * Pnm_part1 - Pnm_part2
        Pnm[n, n + 1] = np.sqrt(2. * n) * Pnm[n, n - 1] / 2.0

    return Pnm


def xyz2dhif(x, y, z):
    """Calculate declination, intensity, inclination of mag field.

    Parameters
    ----------
    x : array-like
        North component of the magnetic field in nT.
    y : array-like
        East component of the magnetic field in nT.
    y : array-like
        Vertical component of the magnetic field in nT.

    Returns
    -------
    dec : array-like
        Declination of the magnetic field in degrees.
    hoz : array-like
        Horizontal intensity of the magnetic field in nT.
    inc : array-like
        Inclination of the magnetic field in degrees.
    eff : array-like
        Total intensity of the magnetic filed in nT.

    Notes
    -----
    by IGRF-13 Alken et al., 2021
    Calculate D, H, I and F from (X, Y, Z)
    Based on code from D. Kerridge, 2019.

    """
    hsq = x * x + y * y
    hoz = np.sqrt(hsq)
    eff = np.sqrt(hsq + z * z)
    dec = np.rad2deg(np.arctan2(y, x))
    inc = np.rad2deg(np.arctan2(z, hoz))

    return dec, hoz, inc, eff

1	#!/usr/bin/env python
2	# ----------------------------------------------------------
3	# Distribution statement A. Approved for public release.
4	# Distribution is unlimited.
5	# This work was supported by the Office of Naval Research.
6	# ----------------------------------------------------------
7	"""This library contains components for IGRF.
8
9	References
10	----------
11	Alken et al. (2021). International geomagnetic reference field: the
12	thirteenth generation. Earth, Planets and Space, 73(49),
13	doi:10.1186/s40623-020-01288-x.
14
15	"""
16
17	import numpy as np	7✔
18	import os	7✔
19	from scipy import interpolate	7✔
20
21
22	def verify_inclination(inc):	7✔
23	"""Test the validity of an inclination input.
24
25	Parameters
26	----------
27	inc : int, float, or array-like
28	Magnetic inclination in degrees.
29
30	Returns
31	-------
32	inc : int, float, or array-like
33	Magnetic inclination in degrees
34
35	Raises
36	------
37	ValueError
38	If inclination is not between -90 and 90 degrees
39
40	"""
41	# Ensure the inclination to be tested is a numpy array
42	test_inc = np.asarray(inc)	7✔
43	if test_inc.shape == ():	7✔
NEW 44	test_inc = np.asarray([inc])	×
45
46	# Get a mask that specifies if each inclination is realistic
47	good_inc = (test_inc <= 90.0) & (test_inc >= -90.0)	7✔
48
49	# Return original input if all are good, correct or raise error otherwise
50	if good_inc.all():	7✔
51	return inc	7✔
52	else:
53	# Adjust values within a small tolerance
NEW 54	test_inc[(test_inc > 90.0) & (test_inc <= 90.1)] = 90.0	×
NEW 55	test_inc[(test_inc < -90.0) & (test_inc >= -90.1)] = -90.0	×
NEW 56	good_inc = (test_inc <= 90.0) & (test_inc >= -90.0)	×
57
NEW 58	if good_inc.all():	×
59	# Return the adjusted input using the original type
NEW 60	if isinstance(inc, (float, int, np.floating, np.integer)):	×
NEW 61	return test_inc[0].astype(type(inc))	×
NEW 62	elif isinstance(inc, list):	×
NEW 63	return list(test_inc)	×
64	else:
NEW 65	return test_inc	×
66	else:
NEW 67	raise ValueError('unrealistic magnetic inclination value(s)')	×
68
69
70	def inc2modip(inc, alat):	7✔
71	"""Calculate modified dip angle from magnetic inclination.
72
73	Parameters
74	----------
75	inc : array-like
76	Magnetic inclination in degrees.
77	alat : array-like
78	Flattened array of latitudes in degrees.
79
80	Returns
81	-------
82	modip_deg : array-like
83	Modified dip angle in degrees.
84
85	Notes
86	-----
87	This function calculates modified dip angle for a given inclination and
88	geographic latitude.
89
90	"""
91	inc = verify_inclination(inc)	7✔
92
93	alat_rad = np.deg2rad(alat)	7✔
94	rad_arg = np.deg2rad((inc / np.sqrt(np.cos(alat_rad))))	7✔
95	modip_rad = np.arctan(rad_arg)	7✔
96	modip_deg = np.rad2deg(modip_rad)	7✔
97
98	return modip_deg	7✔
99
100
101	def inc2magnetic_dip_latitude(inc):	7✔
102	"""Calculate magnetic dip latitude from magnetic inclination.
103
104	Parameters
105	----------
106	inc : array-like
107	Magnetic inclination in degrees.
108
109	Returns
110	-------
111	magnetic_dip_latitude : array-like
112	Magnetic dip latitude in degrees.
113
114	Notes
115	-----
116	This function calculates magnetic dip latitude.
117
118	"""
119	inc = verify_inclination(inc)	7✔
120	arg = 0.5 * (np.tan(np.deg2rad(inc)))	7✔
121	magnetic_dip_latitude = np.rad2deg(np.arctan(arg))	7✔
122
123	return magnetic_dip_latitude	7✔
124
125
126	def inclination(coeff_dir, date_decimal, alon, alat, alt=300., only_inc=True):	7✔
127	"""Calculate magnetic inclination using IGRF13.
128
129	Parameters
130	----------
131	coeff_dir : str
132	Where IGRF13 coefficients are located.
133	date_decimal : float
134	Decimal year
135	alon : array-like
136	Flattened array of geographic longitudes in degrees.
137	alat : array-like
138	Flattened array of geographic latitudes in degrees.
139	alt : flt
140	Altitude in km, default is 300 km.
141	only_inc : bool
142	Only output the inclination, otherwise output the magnetic field
143	information as well (default=True)
144
145	Returns
146	-------
147	inc : array-like
148	Magnetic inclination in degrees.
149	X : array-like (optional)
150	North component of the magnetic field in nT.
151	Y : array-like (optional)
152	East component of the magnetic field in nT.
153	Z : array-like (optional)
154	Vertical component of the magnetic field in nT.
155	dec : array-like (optional)
156	Declination of the magnetic field in degrees.
157	hoz : array-like (optional)
158	Horizontal intensity of the magnetic field in nT.
159	eff : array-like (optional)
160	Total intensity of the magnetic field in nT.
161
162	Raises
163	------
164	IOError
165	If unable to find IGRF coefficient file
166
167	Notes
168	-----
169	This code is a slight modification of the IGRF13 pyIGRF release,
170	https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
171	The reading of the IGRF coefficient file was modified to speded up the
172	process, and the main code was simlified.
173
174	"""
175	# Set altitude
176	aalt = np.full(shape=alon.shape, fill_value=alt)	7✔
177
178	# Open IGRF file and read to array the main table
179	igrf_file = os.path.join(coeff_dir, 'IGRF', 'IGRF13.shc')	7✔
180	if not os.path.isfile(igrf_file):	7✔
NEW 181	raise IOError('unable to find IGRF coefficient file: {:}'.format(	×
182	igrf_file))
183
184	with open(igrf_file, mode='r') as fopen:	7✔
185	file_array = np.genfromtxt(fopen, delimiter='', skip_header=5)	7✔
186
187	# Create the time array of years that match the years in the file
188	# (check if file is change to the newer version)
189	igrf_time = np.arange(1900, 2025 + 5, 5)	7✔
190
191	# Exclude first 2 columns, these are the m and n indecies
192	igrf_coeffs = file_array[:, 2:]	7✔
193
194	# Maximum degree of polynomials
195	nmax = 13	7✔
196
197	# Colatitude (from 0 to 180)
198	colat = 90.0 - alat	7✔
199
200	# Compute geocentric colatitude and radius from geodetic colatitude
201	# and height
202	alt, colat, sd, cd = gg_to_geo(aalt, colat)	7✔
203
204	# Interpolate coefficients from 5 year slots to given decimal year
205	igrf_extrap = interpolate.interp1d(igrf_time, igrf_coeffs,	7✔
206	fill_value='extrapolate')
207	coeffs = igrf_extrap(date_decimal)	7✔
208
209	# Compute the main field B_r, B_theta and B_phi value for the location(s)
210	br, bt, bp = synth_values(coeffs.T, alt, colat, alon, nmax)	7✔
211
212	# Rearrange to X, Y, Z components
213	X = -bt	7✔
214	Y = bp	7✔
215	Z = -br	7✔
216
217	# Rotate back to geodetic coords if needed
218	t = X	7✔
219	X = X * cd + Z * sd	7✔
220	Z = Z * cd - t * sd	7✔
221
222	# Compute the four non-linear components
223	dec, hoz, inc, eff = xyz2dhif(X, Y, Z)	7✔
224
225	if only_inc:	7✔
226	out = inc	7✔
227	else:
NEW 228	out = (inc, X, Y, Z, dec, hoz, eff)	×
229
230	return out	7✔
231
232
233	def gg_to_geo(h, gdcolat):	7✔
234	"""Compute geocentric colatitude and radius from geodetic colat and height.
235
236	Parameters
237	----------
238	h : array-like
239	Altitude in kilometers.
240	gdcolat : array-like
241	Geodetic colatitude in degrees.
242
243	Returns
244	-------
245	rad : array-like
246	Geocentric radius in kilometers.
247	thc : array-like
248	Geocentric colatitude in degrees.
249	sd : array-like
250	Rotate B_X to gd_lat.
251	cd : array-like
252	Rotate B_Z to gd_lat.
253
254	Notes
255	-----
256	IGRF-13 Alken et al., 2021.
257
258	References
259	----------
260	Equations (51)-(53) from "The main field" (chapter 4) by Langel,
261	R. A. in: "Geomagnetism", Volume 1, Jacobs, J. A., Academic Press,
262	1987.
263
264	Malin, S.R.C. and Barraclough, D.R., 1981. An algorithm for
265	synthesizing the geomagnetic field. Computers & Geosciences, 7(4),
266	pp. 401-405.
267
268	"""
269	# Use WGS-84 ellipsoid parameters
270	eqrad = 6378.137 # equatorial radius	7✔
271	flat = 1. / 298.257223563	7✔
272	plrad = eqrad * (1. - flat) # polar radius	7✔
273	ctgd = np.cos(np.deg2rad(gdcolat))	7✔
274	stgd = np.sin(np.deg2rad(gdcolat))	7✔
275	a2 = eqrad * eqrad	7✔
276	a4 = a2 * a2	7✔
277	b2 = plrad * plrad	7✔
278	b4 = b2 * b2	7✔
279	c2 = ctgd * ctgd	7✔
280	s2 = 1. - c2	7✔
281	rho = np.sqrt(a2 * s2 + b2 * c2)	7✔
282	rad = np.sqrt(h * (h + 2. * rho) + (a4 * s2 + b4 * c2) / rho**2)	7✔
283	cd = (h + rho) / rad	7✔
284	sd = (a2 - b2) * ctgd * stgd / (rho * rad)	7✔
285	cthc = ctgd * cd - stgd * sd	7✔
286
287	# Clip cthc to [-1, 1] not to cause problems in arccos
288	cthc = np.clip(cthc, -1.0, 1.0)	7✔
289
290	# Arccos returns values in [0, pi] so we need to convert to degrees
291	thc = np.rad2deg(np.arccos(cthc)) # arccos returns values in [0, pi]	7✔
292
293	return rad, thc, sd, cd	7✔
294
295
296	def geo_to_gg(radius, theta):	7✔
297	"""Compute geodetic colatitude and vertical height above the ellipsoid.
298
299	Parameters
300	----------
301	radius : array-like
302	Geocentric radius in kilometers.
303	theta : array-like
304	Geocentric colatitude in degrees.
305
306	Returns
307	-------
308	height : array-like
309	Altitude in kilometers.
310	beta : array-like
311	Geodetic colatitude.
312
313	Notes
314	-----
315	IGRF-13 Alken et al., 2021.
316	Compute geodetic colatitude and vertical height above the ellipsoid from
317	geocentric radius and colatitude. Round-off errors might lead to a
318	failure of the algorithm especially but not exclusively for points close
319	to the geographic poles. Corresponding geodetic coordinates are returned
320	as NaN.
321
322	References
323	----------
324	Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
325	geodetic coordinates", IEEE Transactions on Aerospace and Electronic
326	Systems}, 1994, vol. 30, num. 3, pp. 957-961.
327
328	"""
329	# Use WGS-84 ellipsoid parameters
330	a = 6378.137 # equatorial radius	×
331	b = 6356.752 # polar radius	×
332	a2 = a**2	×
333	b2 = b**2	×
334	e2 = (a2 - b2) / a2 # squared eccentricity	×
335	e4 = e2 * e2	×
336	ep2 = (a2 - b2) / b2 # squared primed eccentricity	×
337	r = radius * np.sin(np.radians(theta))	×
338	z = radius * np.cos(np.radians(theta))	×
339	r2 = r**2	×
340	z2 = z**2	×
341	F = 54. * b2 * z2	×
342	G = r2 + (1. - e2) * z2 - e2 * (a2 - b2)	×
343	c = e4 * F * r2 / G**3	×
344	s = (1. + c + np.sqrt(c*2 + 2 c))**(1. / 3.)	×
345	P = F / (3 * (s + 1. / s + 1.)*2 G**2)	×
346	Q = np.sqrt(1. + 2 * e4 * P)	×
347	r0 = (-P * e2 * r / (1. + Q)	×
348	+ np.sqrt(0.5 * a2 * (1. + 1. / Q) - P * (1. - e2) * z2
349	/ (Q * (1. + Q)) - 0.5 * P * r2))
350	U = np.sqrt((r - e2 * r0)**2 + z2)	×
351	V = np.sqrt((r - e2 * r0)*2 + (1. - e2) z2)	×
352	z0 = b2 * z / (a * V)	×
353	height = U * (1. - b2 / (a * V))	×
354	beta = 90. - np.degrees(np.arctan2(z + ep2 * z0, r))	×
355
356	return height, beta	×
357
358
359	def synth_values(coeffs, radius, theta, phi, nmax=None, nmin=1, grid=False):	7✔
360	"""Compute radial, colatitude and azimuthal field components.
361
362	Parameters
363	----------
364	coeffs : array-like
365	Coefficients of the spherical harmonic expansion. The last
366	dimension is equal to the number of coefficients, `N` at the grid
367	points.
368	radius : array-like
369	Array containing the radius in km.
370	theta : array-like
371	Array containing the colatitude in degrees.
372	phi : array-like
373	Array containing the longitude in degrees.
374	nmax : int or NoneType
375	Maximum degree up to which expansion is to be used, if None it will be
376	specified by the ``coeffs`` variable. However, smaller values are
377	also valid if specified. (default=None)
378	nmin : int
379	Minimum degree from which expansion is to be used. Note that it will
380	just skip the degrees smaller than ``nmin``, the whole sequence of
381	coefficients 1 through ``nmax`` must still be given in ``coeffs``.
382	(default=1)
383	grid : bool
384	If ``True``, field components are computed on a regular grid. Arrays
385	``theta`` and ``phi`` must have one dimension less than the output
386	grid since the grid will be created as their outer product
387	(default=False).
388
389	Returns
390	-------
391	B_radius : array-like
392	Radial field components in km.
393	B_theta : array-like
394	Colatitude field components in degrees.
395	B_phi : array-like
396	Azimuthal field components in degrees.
397
398	Raises
399	------
400	ValueError
401	If an inappropriate input value is supplied
402
403	Notes
404	-----
405	by IGRF-13 Alken et al., 2021
406	Based on chaosmagpy from Clemens Kloss (DTU Space, Copenhagen)
407	Computes radial, colatitude and azimuthal field components from the
408	magnetic potential field in terms of spherical harmonic coefficients.
409	A reduced version of the DTU synth_values chaosmagpy code.
410
411	References
412	----------
413	Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
414	geodetic coordinates", IEEE Transactions on Aerospace and Electronic
415	Systems}, 1994, vol. 30, num. 3, pp. 957-961.
416
417	"""
418	# Ensure ndarray inputs
419	coeffs = np.array(coeffs, dtype=float)	7✔
420	radius = np.array(radius, dtype=float) / 6371.2 # Earth's average radius	7✔
421	theta = np.array(theta, dtype=float)	7✔
422	phi = np.array(phi, dtype=float)	7✔
423	if (np.amin(theta) < 0.0) \| (np.amax(theta) > 180.0):	7✔
424	raise ValueError('Colatitude outside bounds [0, 180].')	×
425
426	if nmin < 1:	7✔
427	raise ValueError('Only positive nmin allowed.')	×
428
429	# Handle optional argument: nmax
430	nmax_coeffs = int(np.sqrt(coeffs.shape[-1] + 1) - 1) # degrees	7✔
431	if nmax is None:	7✔
432	nmax = nmax_coeffs	×
433	elif nmax <= 0:	7✔
NEW 434	raise ValueError('Only positive, non-zero nmax allowed.')	×
435
436	if nmax > nmax_coeffs:	7✔
437	nmax = nmax_coeffs	×
438
439	if nmax < nmin:	7✔
440	raise ValueError(f'Nothing to compute: nmax < nmin ({nmax} < {nmin}.)')	×
441
442	# Manually broadcast input grid on surface
443	if grid:	7✔
NEW 444	theta = theta[..., None] # First dimension is theta	×
NEW 445	phi = phi[None, ...] # Second dimension is phi	×
446
447	# Get shape of broadcasted result
448	try:	7✔
449	broad = np.broadcast(radius, theta, phi,	7✔
450	np.broadcast_to(0, coeffs.shape[:-1]))
451	except ValueError:	×
452	raise ValueError(''.join(['Cannot broadcast grid shapes (excl. last ',	×
453	'dimension of coeffs):\nradius: ',
454	repr(radius.shape), '\ntheta: ',
455	repr(theta.shape), '\nphi: ', repr(phi.shape),
456	'\ncoeffs: ', repr(coeffs.shape)]))
457	grid_shape = broad.shape	7✔
458
459	# Initialize radial dependence given the source
460	r_n = radius**(-(nmin + 2))	7✔
461
462	# Compute associated Legendre polynomials as (n, m, theta-points)-array
463	Pnm = legendre_poly(nmax, theta)	7✔
464
465	# Save sinth for fast access
466	sinth = Pnm[1, 1]	7✔
467
468	# Calculate cos(mphi) and sin(mphi) as (m, phi-points)-array
469	phi = np.radians(phi)	7✔
470	cos_mp = np.cos(np.multiply.outer(np.arange(nmax + 1), phi))	7✔
471	sin_mp = np.sin(np.multiply.outer(np.arange(nmax + 1), phi))	7✔
472
473	# allocate arrays in memory
474	B_radius = np.zeros(grid_shape)	7✔
475	B_theta = np.zeros(grid_shape)	7✔
476	B_phi = np.zeros(grid_shape)	7✔
477	num = nmin**2 - 1	7✔
478	for n in range(nmin, nmax + 1):	7✔
479	B_radius += (n + 1) * Pnm[n, 0] * r_n * coeffs[..., num]	7✔
480	B_theta += -Pnm[0, n + 1] * r_n * coeffs[..., num]	7✔
481	num += 1	7✔
482	for m in range(1, n + 1):	7✔
483	B_radius += ((n + 1) * Pnm[n, m] * r_n	7✔
484	* (coeffs[..., num] * cos_mp[m]
485	+ coeffs[..., num + 1] * sin_mp[m]))
486	B_theta += (-Pnm[m, n + 1] * r_n	7✔
487	* (coeffs[..., num] * cos_mp[m]
488	+ coeffs[..., num + 1] * sin_mp[m]))
489	with np.errstate(divide='ignore', invalid='ignore'):	7✔
490	# Handle poles using L'Hopital's rule
491	div_Pnm = np.where(theta == 0., Pnm[m, n + 1],	7✔
492	Pnm[n, m] / sinth)
493	div_Pnm = np.where(theta == np.degrees(np.pi),	7✔
494	-Pnm[m, n + 1], div_Pnm)
495	B_phi += (m * div_Pnm * r_n	7✔
496	* (coeffs[..., num] * sin_mp[m]
497	- coeffs[..., num + 1] * cos_mp[m]))
498	num += 2	7✔
499	r_n = r_n / radius # Equivalent to r_n = radius**(-(n+2))	7✔
500
501	return B_radius, B_theta, B_phi	7✔
502
503
504	def legendre_poly(nmax, theta):	7✔
505	"""Calculate associated Legendre polynomials P(n,m).
506
507	Parameters
508	----------
509	nmax : int
510	Maximum degree up to which expansion, with a minimum value of one.
511	theta : array-like
512	Array containing the colatitude in degrees.
513
514	Returns
515	-------
516	Pnm : array-like
517	Evaluated values and derivatives.
518
519	Raises
520	------
521	IndexError
522	If nmax is less than 1
523
524	Notes
525	-----
526	by IGRF-13 Alken et al., 2021
527	Returns associated Legendre polynomials P(n,m) (Schmidt
528	quasi-normalized), and the derivative :d P(n,m) / d theta evaluated
529	at theta.
530
531	References
532	----------
533	Langel, R. A. "Chapter four: Main field." Geomagnetism, edited by JA Jacobs
534	1 (1987).
535	Zhu, J., "Conversion of Earth-centered Earth-fixed coordinates to
536	geodetic coordinates", IEEE Transactions on Aerospace and Electronic
537	Systems}, 1994, vol. 30, num. 3, pp. 957-961.
538
539	"""
540	costh = np.cos(np.radians(theta))	7✔
541	sinth = np.sqrt(1 - costh**2)	7✔
542	Pnm = np.zeros((nmax + 1, nmax + 2) + costh.shape)	7✔
543	Pnm[0, 0] = 1 # is copied into trailing dimenions	7✔
544	Pnm[1, 1] = sinth # write theta into trailing dimenions via broadcasting	7✔
545	rootn = np.sqrt(np.arange(2 * nmax**2 + 1))	7✔
546
547	# Recursion relations after Langel "The Main Field" (1987),
548	# eq. (27) and Table 2 (p. 256)
549	for m in range(nmax):	7✔
550	Pnm_tmp = rootn[m + m + 1] * Pnm[m, m]	7✔
551	Pnm[m + 1, m] = costh * Pnm_tmp	7✔
552	if m > 0:	7✔
553	Pnm[m + 1, m + 1] = sinth * Pnm_tmp / rootn[m + m + 2]	7✔
554	for n in np.arange(m + 2, nmax + 1):	7✔
555	d = n * n - m * m	7✔
556	e = n + n - 1	7✔
557	Pnm[n, m] = ((e * costh * Pnm[n - 1, m]	7✔
558	- rootn[d - e] * Pnm[n - 2, m]) / rootn[d])
559
560	# dP(n,m) = Pnm(m,n+1) is the derivative of P(n,m) vrt. theta
561	Pnm[0, 2] = -Pnm[1, 1]	7✔
562	Pnm[1, 2] = Pnm[1, 0]	7✔
563	for n in range(2, nmax + 1):	7✔
564	n = int(n)	7✔
565	Pnm[0, n + 1] = - np.sqrt((n * n + n) / 2.) * Pnm[n, 1]	7✔
566	Pnm[1, n + 1] = ((np.sqrt(2.0 * (n * n + n)) * Pnm[n, 0]	7✔
567	- np.sqrt((n * n + n - 2)) * Pnm[n, 2]) / 2.0)
568	for m in np.arange(2, n):	7✔
569	m = int(m)	7✔
570	Pnm_part1 = np.sqrt((n + m) * (n - m + 1)) * Pnm[n, m - 1]	7✔
571	Pnm_part2 = np.sqrt((n + m + 1.0) * (n - m)) * Pnm[n, m + 1]	7✔
572	Pnm[m, n + 1] = 0.5 * Pnm_part1 - Pnm_part2	7✔
573	Pnm[n, n + 1] = np.sqrt(2. * n) * Pnm[n, n - 1] / 2.0	7✔
574
575	return Pnm	7✔
576
577
578	def xyz2dhif(x, y, z):	7✔
579	"""Calculate declination, intensity, inclination of mag field.
580
581	Parameters
582	----------
583	x : array-like
584	North component of the magnetic field in nT.
585	y : array-like
586	East component of the magnetic field in nT.
587	y : array-like
588	Vertical component of the magnetic field in nT.
589
590	Returns
591	-------
592	dec : array-like
593	Declination of the magnetic field in degrees.
594	hoz : array-like
595	Horizontal intensity of the magnetic field in nT.
596	inc : array-like
597	Inclination of the magnetic field in degrees.
598	eff : array-like
599	Total intensity of the magnetic filed in nT.
600
601	Notes
602	-----
603	by IGRF-13 Alken et al., 2021
604	Calculate D, H, I and F from (X, Y, Z)
605	Based on code from D. Kerridge, 2019.
606
607	"""
608	hsq = x * x + y * y	7✔
609	hoz = np.sqrt(hsq)	7✔
610	eff = np.sqrt(hsq + z * z)	7✔
611	dec = np.rad2deg(np.arctan2(y, x))	7✔
612	inc = np.rad2deg(np.arctan2(z, hoz))	7✔
613
614	return dec, hoz, inc, eff	7✔

victoriyaforsythe / PyIRI / 15284544330

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous