Skip to content

Current source modelling

Package Module scubas.sources.current_sources Class Direction, Current

This module models ionospheric line-current sources and computes associated magnetic and electric perturbations.

Unit conventions

Default inputs use km for length and MA for current. Internally these values are converted to SI units (m, A) before field calculations.

Quick start

import numpy as np
from scubas.sources.current_sources import Current, Direction

x = np.linspace(-500, 500, 201)  # km
source = Current(
    h=110.0,
    I=1.0,
    current_direction=Direction.GEO_EAST_WEST,
)
source.compute_B_field(x=x, length_unit="km")

print(source.B_fields[0]["direction"], source.B_field_units)
print(source.E_fields[0]["direction"], source.E_field_units)

API reference

scubas.sources.current_sources.Direction

Bases: Enum

Source code in scubas/sources/current_sources.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
class Direction(Enum):
    GEO_NORTH_SOUTH = 1
    GEO_EAST_WEST = 2
    GEO_VERT_DOWN = 3

    @classmethod
    def enumerate_key(cls, parameter: str = "B"):
        """
        Enumerate possible keys based on parameter
        """
        if cls == Direction.GEO_EAST_WEST:
            key = rf"${parameter}$_y"
        if cls == Direction.GEO_NORTH_SOUTH:
            key = rf"${parameter}$_x"
        if cls == Direction.GEO_VERT_DOWN:
            key = rf"${parameter}$_z"
        return key

enumerate_key(parameter='B') classmethod

Enumerate possible keys based on parameter

Source code in scubas/sources/current_sources.py 
30
31
32
33
34
35
36
37
38
39
40
41
@classmethod
def enumerate_key(cls, parameter: str = "B"):
    """
    Enumerate possible keys based on parameter
    """
    if cls == Direction.GEO_EAST_WEST:
        key = rf"${parameter}$_y"
    if cls == Direction.GEO_NORTH_SOUTH:
        key = rf"${parameter}$_x"
    if cls == Direction.GEO_VERT_DOWN:
        key = rf"${parameter}$_z"
    return key

scubas.sources.current_sources.Current

Bases: object

Source code in scubas/sources/current_sources.py 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
class Current(object):
    def __init__(
        self,
        h: float = 110.0,
        a: float = 0.0,
        I: float = 1.0,
        p: np.array = 0.0 * np.arange(10, 200),
        omega: np.array = 2 * np.pi * np.arange(10, 200),
        site: Site = PROFILES.Quebec,
        length_unit: str = "km",
        current_unit: str = "MA",
        current_direction: Direction = Direction.GEO_EAST_WEST,
        t: int = 0,
        layer: int = 0,
        return_Z: bool = False,
    ):
        """
        Current systems for one time instance / one location
            h: Height given in unit
            a: Width of the current, 'Cauchy Distributed'
            I: Currents in Amps/M Amps
            p: List of skin depths of the wave Z/(j*omega)
            omega: List of angular frequencies in rad/s
            site: Site profile default Quebec
            length_unit: Units for h, p, and a
            current_unit: Units of I
            current_direction: Direction of the currents, in geographic coordinates
            t: Time instance of the current
            layer: Layer where calculate Z and p
        """
        self.length_unit = length_unit
        self.current_unit = current_unit
        self.length_unit_multiplier = 1e3 if length_unit == "km" else 1.0
        self.current_unit_multiplier = 1e6 if current_unit == "MA" else 1.0
        self.h = h * self.length_unit_multiplier  # in meters
        self.p = p * self.length_unit_multiplier  # in meters
        self.a = a * self.length_unit_multiplier  # in meters
        self.I = I * self.current_unit_multiplier  # in Amps
        self.omega = omega
        self.freq = omega / (2 * np.pi)
        self.site = site
        self.current_direction = current_direction
        self.t = t
        # If site is given calculate the skin depth p at top layer
        if site:
            ret = self.site.calcP(self.freq, layer=layer, return_Z=return_Z)
            if return_Z:
                self.p, self.Z_out, self.Z = ret
            else:
                self.p = ret
        return

    def compute_B_field(
        self,
        x: np.array = np.zeros((1)),
        length_unit: str = "km",
    ):
        """
        For a line current along y directed, we can only expect
        magnetic field along x and z directed.
        x: List of distance along x-direction given in unit
        length_unit: Units for x
        """
        length_unit_multiplier = 1e3 if length_unit == "km" else 1.0
        x *= length_unit_multiplier  # to meters
        B_field_directions = [
            Direction.GEO_EAST_WEST,
            Direction.GEO_NORTH_SOUTH,
            Direction.GEO_VERT_DOWN,
        ]
        # Remove 'y' directions only from the magnetic fields
        B_field_directions.remove(self.current_direction)
        logger.info(f"|Current| along : {self.current_direction} is {self.I} A")
        logger.info(
            f"Direction of the produced magnetic fields: {[d for d in B_field_directions]}"
        )
        # Modified height by Cauchy distributed current source
        h, h_mod, factor = (
            self.a + self.h,
            self.a + self.h + (2 * self.p),
            (C.mu_0 * self.I / (2 * np.pi)),
        )
        self.B_fields = [
            dict(
                direction=B_field_directions[0],
                total=factor
                * ((h / (h**2 + x**2)) + (h_mod / (h_mod**2 + x**2))),
                external=factor * (h / (h**2 + x**2)),
                internal=factor * (h_mod / (h_mod**2 + x**2)),
            ),
            dict(
                direction=B_field_directions[1],
                total=-1
                * factor
                * ((x / (h**2 + x**2)) - (x / (h_mod**2 + x**2))),
                external=-1 * factor * (x / (h**2 + x**2)),
                internal=factor * (x / (h_mod**2 + x**2)),
            ),
        ]
        for b in self.B_fields:
            logger.warning(
                f"Magnetic field B along {b['direction']} is {b['total'].min()},{b['total'].max()}"
            )
        self.__compute_E_field__(x, "m")
        self.B_field_units = "T"
        return

    def __compute_E_field__(
        self,
        x: np.array = np.zeros((1)),
        length_unit: str = "km",
    ):
        """
        For a line current along y directed, we can only expect
        magnetic field along y directed.
        """
        length_unit_multiplier = 1e3 if length_unit == "km" else 1.0
        x *= length_unit_multiplier  # to meters
        E_field_directions = [self.current_direction]
        logger.info(f"Direction of the E_field: {E_field_directions}")
        # Modified height by Cauchy distributed current source
        h = self.a + self.h
        self.E_fields = [
            dict(
                direction=E_field_directions[0],
                total=-1j
                * (self.omega * C.mu_0 * self.I / (2 * np.pi))
                * np.log(
                    np.sqrt((h + 2 * self.p) ** 2 + x**2) / np.sqrt(h**2 + x**2)
                ),
            ),
        ]
        for e in self.E_fields:
            logger.warning(
                f"Magnetic field E along {e['direction']} is {e['total'].min()},{e['total'].max()}"
            )
        self.E_field_units = "V/m"
        return

compute_B_field(x=np.zeros(1), length_unit='km')

For a line current along y directed, we can only expect magnetic field along x and z directed. x: List of distance along x-direction given in unit length_unit: Units for x

Source code in scubas/sources/current_sources.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def compute_B_field(
    self,
    x: np.array = np.zeros((1)),
    length_unit: str = "km",
):
    """
    For a line current along y directed, we can only expect
    magnetic field along x and z directed.
    x: List of distance along x-direction given in unit
    length_unit: Units for x
    """
    length_unit_multiplier = 1e3 if length_unit == "km" else 1.0
    x *= length_unit_multiplier  # to meters
    B_field_directions = [
        Direction.GEO_EAST_WEST,
        Direction.GEO_NORTH_SOUTH,
        Direction.GEO_VERT_DOWN,
    ]
    # Remove 'y' directions only from the magnetic fields
    B_field_directions.remove(self.current_direction)
    logger.info(f"|Current| along : {self.current_direction} is {self.I} A")
    logger.info(
        f"Direction of the produced magnetic fields: {[d for d in B_field_directions]}"
    )
    # Modified height by Cauchy distributed current source
    h, h_mod, factor = (
        self.a + self.h,
        self.a + self.h + (2 * self.p),
        (C.mu_0 * self.I / (2 * np.pi)),
    )
    self.B_fields = [
        dict(
            direction=B_field_directions[0],
            total=factor
            * ((h / (h**2 + x**2)) + (h_mod / (h_mod**2 + x**2))),
            external=factor * (h / (h**2 + x**2)),
            internal=factor * (h_mod / (h_mod**2 + x**2)),
        ),
        dict(
            direction=B_field_directions[1],
            total=-1
            * factor
            * ((x / (h**2 + x**2)) - (x / (h_mod**2 + x**2))),
            external=-1 * factor * (x / (h**2 + x**2)),
            internal=factor * (x / (h_mod**2 + x**2)),
        ),
    ]
    for b in self.B_fields:
        logger.warning(
            f"Magnetic field B along {b['direction']} is {b['total'].min()},{b['total'].max()}"
        )
    self.__compute_E_field__(x, "m")
    self.B_field_units = "T"
    return