Coverage for magnetismi / model / cof.py: 97.39%
99 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-10 18:29:16 +00:00
1import importlib.resources
2import math
3import os
4import pathlib
5from typing import Any, ContextManager, no_type_check
7from magnetismi import ENCODING
9CS = ('n', 'm', 'gnm', 'hnm', 'dgnm', 'dhnm')
10STORE_PATH = pathlib.Path('magnetismi', 'model')
11TWIN_END_TOKEN = '9' * 48
13YEARS_COVERED = tuple(y for y in range(2020, 2030 + 1))
14MODEL_FROM_YEAR = {y: '2020' if y < 2026 else '2025' for y in range(2020, 2030 + 1)}
17def _normalize_path(path: Any) -> str:
18 """Normalize a path by ensuring it is a string.
20 If the resulting string contains path separators, an exception is raised.
21 """
22 str_path = str(path)
23 parent, file_name = os.path.split(str_path)
24 if parent: 24 ↛ 25line 24 didn't jump to line 25 because the condition on line 24 was never true
25 raise ValueError(f'{path!r} must be only a file name')
26 return file_name
29@no_type_check
30def _path(package, resource) -> ContextManager[pathlib.Path]:
31 """A context manager providing a file path object to the resource.
33 If the resource does not already exist on its own on the file system,
34 a temporary file will be created. If the file was created, the file
35 will be deleted upon exiting the context manager (no exception is
36 raised if the file was deleted prior to the context manager
37 exiting).
38 """
39 return importlib.resources._common.as_file(importlib.resources._common.files(package) / _normalize_path(resource))
42class Coefficients:
43 """The WMM coefficients relevant to the year given."""
45 model = {} # type: ignore
47 def load(self, model_year: int) -> None:
48 """Load the model."""
49 model_resource = f'wmm-{model_year}.txt'
50 with _path(__package__, model_resource) as model_path:
51 with open(model_path, 'rt', encoding=ENCODING) as handle:
52 recs = [line.strip().split() for line in handle if line.strip() and line.strip() != TWIN_END_TOKEN]
54 epoc, model, model_date_string = recs[0]
55 self.model['epoch'] = float(epoc)
56 self.model['model_code'] = model
57 self.model['model_date_string'] = model_date_string
59 self.model['coeffs'] = [
60 {c: int(r[s]) if c in ('n', 'm') else float(r[s]) for s, c in enumerate(CS)}
61 for r in recs[1:]
62 if len(r) == len(CS)
63 ]
65 def _zeroes(self, count: int) -> list[float]:
66 """Short and sweet."""
67 return [0.0 for _ in range(count)]
69 def evaluate(self) -> None:
70 """Ja, ja, ja!"""
71 self.maxord = self.maxdeg = 12
72 self.tc = [self._zeroes(13) for _ in range(14)]
73 self.sp = self._zeroes(14)
74 self.cp = self._zeroes(14)
75 self.cp[0] = 1.0
76 self.pp = self._zeroes(13)
77 self.pp[0] = 1.0
78 self.p = [self._zeroes(14) for _ in range(14)]
79 self.p[0][0] = 1.0
80 self.dp = [self._zeroes(13) for _ in range(14)]
81 self.a = 6378.137
82 self.b = 6356.7523142
83 self.re = 6371.2
84 self.a2 = self.a * self.a
85 self.b2 = self.b * self.b
86 self.c2 = self.a2 - self.b2
87 self.a4 = self.a2 * self.a2
88 self.b4 = self.b2 * self.b2
89 self.c4 = self.a4 - self.b4
91 self.c = [self._zeroes(14) for _ in range(14)]
92 self.cd = [self._zeroes(14) for _ in range(14)]
94 for wmmnm in self.model['coeffs']:
95 m = wmmnm['m']
96 n = wmmnm['n']
97 gnm = wmmnm['gnm']
98 hnm = wmmnm['hnm']
99 dgnm = wmmnm['dgnm']
100 dhnm = wmmnm['dhnm']
101 if m <= n: 101 ↛ 94line 101 didn't jump to line 94 because the condition on line 101 was always true
102 self.c[m][n] = gnm
103 self.cd[m][n] = dgnm
104 if m != 0:
105 self.c[n][m - 1] = hnm
106 self.cd[n][m - 1] = dhnm
108 # Convert Schmidt normalized Gauss coefficients to unnormalized
109 self.snorm = [self._zeroes(13) for _ in range(13)]
110 self.snorm[0][0] = 1.0
111 self.k = [self._zeroes(13) for _ in range(13)]
112 self.k[1][1] = 0.0
113 self.fn = [float(n) for n in range(14) if n != 1]
114 self.fm = [float(n) for n in range(13)]
115 for n in range(1, self.maxord + 1):
116 self.snorm[0][n] = self.snorm[0][n - 1] * (2.0 * n - 1) / n
117 j = 2.0
118 # for (m=0,D1=1,D2=(n-m+D1)/D1;D2>0;D2--,m+=D1):
119 m = 0
120 D1 = 1
121 D2 = (n - m + D1) / D1
122 while D2 > 0:
123 self.k[m][n] = (((n - 1) * (n - 1)) - (m * m)) / ((2.0 * n - 1) * (2.0 * n - 3.0))
124 if m > 0:
125 flnmj = ((n - m + 1.0) * j) / (n + m)
126 self.snorm[m][n] = self.snorm[m - 1][n] * math.sqrt(flnmj)
127 j = 1.0
128 self.c[n][m - 1] = self.snorm[m][n] * self.c[n][m - 1]
129 self.cd[n][m - 1] = self.snorm[m][n] * self.cd[n][m - 1]
130 self.c[m][n] = self.snorm[m][n] * self.c[m][n]
131 self.cd[m][n] = self.snorm[m][n] * self.cd[m][n]
132 D2 = D2 - 1
133 m = m + D1
135 def __init__(self, year: int):
136 """Initialize the mode from file in case it covers the requested year."""
137 if year not in YEARS_COVERED:
138 raise ValueError(f'requested year ({year}) not within {YEARS_COVERED}')
140 self.model = {
141 'year_requested': year,
142 'model_id': MODEL_FROM_YEAR[year],
143 }
144 self.load(self.model['model_id'])
145 self.evaluate()