Coverage for magnetismi/magnetismi.py: 79.12%
142 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-04 20:21:25 +00:00
1"""Adapted from C code as distributed per https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml."""
3import argparse
4import datetime as dti
5import logging
6import math
7from dataclasses import dataclass
8from typing import no_type_check
10from magnetismi import ABS_LAT_DD_ANY_ARCITC, DEFAULT_MAG_VAR, FEET_TO_KILOMETER, fractional_year_from_date, log
11from magnetismi.model.cof import MODEL_FROM_YEAR, YEARS_COVERED, Coefficients
14@dataclass
15class MagneticVector:
16 dec: float
17 dip: float
18 ti: float
19 bh: float
20 bx: float
21 by: float
22 bz: float
23 lat: float
24 lon: float
25 alt: float
26 time: float
29class Model:
30 """Yes."""
32 def __init__(self, year: int):
33 """Maybe."""
34 self.coefficients = Coefficients(year)
36 @no_type_check
37 def at(self, lat_dd: float, lon_dd: float, alt_ft: float = 0.0, date: dti.date = dti.date.today()):
38 """Get the values at ..."""
39 if date.year not in YEARS_COVERED or self.coefficients.model['model_id'] != MODEL_FROM_YEAR[date.year]: 39 ↛ 40line 39 didn't jump to line 40, because the condition on line 39 was never true
40 raise ValueError(f'Model ID ({self.coefficients.model["model_id"]}) is not covering year ({date.year})')
42 time = fractional_year_from_date(date) # date('Y') + date('z')/365
43 alt = alt_ft * FEET_TO_KILOMETER
45 model = self.coefficients.model
46 coeffs = self.coefficients
47 dt = time - model['epoch']
48 glat = lat_dd
49 glon = lon_dd
50 rlat = math.radians(glat)
51 rlon = math.radians(glon)
52 srlon = math.sin(rlon)
53 srlat = math.sin(rlat)
54 crlon = math.cos(rlon)
55 crlat = math.cos(rlat)
56 srlat_2 = srlat * srlat
57 crlat_2 = crlat * crlat
58 coeffs.sp[1] = srlon
59 coeffs.cp[1] = crlon
61 # watch these:
62 r = ca = sa = 1.0
63 st = ct = 0.0
65 # Convert from geodetic to spherical coordinates
66 q = math.sqrt(coeffs.a2 - coeffs.c2 * srlat_2)
67 q1 = alt * q
68 q2 = ((q1 + coeffs.a2) / (q1 + coeffs.b2)) * ((q1 + coeffs.a2) / (q1 + coeffs.b2))
69 ct = srlat / math.sqrt(q2 * crlat_2 + srlat_2)
70 st = math.sqrt(1.0 - (ct * ct))
71 r2 = (alt * alt) + 2.0 * q1 + (coeffs.a4 - coeffs.c4 * srlat_2) / (q * q)
72 r = math.sqrt(r2)
73 d = math.sqrt(coeffs.a2 * crlat_2 + coeffs.b2 * srlat_2)
74 ca = (alt + d) / r
75 sa = coeffs.c2 * crlat * srlat / (r * d)
77 for m in range(2, coeffs.maxord + 1):
78 coeffs.sp[m] = coeffs.sp[1] * coeffs.cp[m - 1] + coeffs.cp[1] * coeffs.sp[m - 1]
79 coeffs.cp[m] = coeffs.cp[1] * coeffs.cp[m - 1] - coeffs.sp[1] * coeffs.sp[m - 1]
81 aor = coeffs.re / r
82 ar = aor * aor
83 br = bt = bp = bpp = 0.0
84 for n in range(1, coeffs.maxord + 1):
85 ar *= aor
87 m = 0
88 D4 = n + m + 1 # (n+m+D3)/D3 with D3 = 1
89 while D4 > 0: # for (m=0,D3=1,D4=(n+m+D3)/D3;D4>0;D4--,m+=D3):
90 # Compute unnormalized associated Legendre polynomials and derivatives from recursion formula
91 if n == m:
92 coeffs.p[m][n] = st * coeffs.p[m - 1][n - 1]
93 coeffs.dp[m][n] = st * coeffs.dp[m - 1][n - 1] + ct * coeffs.p[m - 1][n - 1]
94 elif n == 1 and m == 0:
95 coeffs.p[m][n] = ct * coeffs.p[m][n - 1]
96 coeffs.dp[m][n] = ct * coeffs.dp[m][n - 1] - st * coeffs.p[m][n - 1]
97 elif n > 1 and n != m: 97 ↛ 107line 97 didn't jump to line 107, because the condition on line 97 was never false
98 if m > n - 2:
99 coeffs.p[m][n - 2] = 0
100 coeffs.dp[m][n - 2] = 0.0
101 coeffs.p[m][n] = ct * coeffs.p[m][n - 1] - coeffs.k[m][n] * coeffs.p[m][n - 2]
102 coeffs.dp[m][n] = (
103 ct * coeffs.dp[m][n - 1] - st * coeffs.p[m][n - 1] - coeffs.k[m][n] * coeffs.dp[m][n - 2]
104 )
106 # Time adjust the Gauss coefficients
107 coeffs.tc[m][n] = coeffs.c[m][n] + dt * coeffs.cd[m][n]
108 if m != 0:
109 coeffs.tc[n][m - 1] = coeffs.c[n][m - 1] + dt * coeffs.cd[n][m - 1]
111 # Accumulate terms of the spherical harmonic expansions
112 par = ar * coeffs.p[m][n]
114 temp1 = coeffs.tc[m][n] * coeffs.cp[m]
115 temp2 = coeffs.tc[m][n] * coeffs.sp[m]
116 if n != 0: 116 ↛ 120line 116 didn't jump to line 120, because the condition on line 116 was never false
117 temp1 += coeffs.tc[n][m - 1] * coeffs.sp[m]
118 temp2 -= coeffs.tc[n][m - 1] * coeffs.cp[m]
120 bt -= ar * temp1 * coeffs.dp[m][n]
121 bp += coeffs.fm[m] * temp2 * par
122 br += coeffs.fn[n] * temp1 * par
124 # Special case: North/South geographic poles
125 if st == 0.0 and m == 1: 125 ↛ 126line 125 didn't jump to line 126
126 coeffs.pp[n] = (
127 coeffs.pp[n - 1] if n == 1 else ct * coeffs.pp[n - 1] - coeffs.k[m][n] * coeffs.pp[n - 2]
128 )
129 parp = ar * coeffs.pp[n]
130 bpp += coeffs.fm[m] * temp2 * parp
132 D4 -= 1
133 m += 1
135 bp = bpp if st == 0.0 else bp / st
137 # Rotate magnetic vector components from spherical to geodetic coordinates
138 bx = -bt * ca - br * sa
139 by = bp
140 bz = bt * sa - br * ca
142 # Compute declination (dec), inclination (dip), and total intensity (ti)
143 bh = math.sqrt((bx * bx) + (by * by))
144 ti = math.sqrt((bh * bh) + (bz * bz))
145 dec = math.degrees(math.atan2(by, bx))
146 dip = math.degrees(math.atan2(bz, bh))
148 # Either compute magnetic grid variation if the current geodetic position is in the arctic or antarctic
149 # i.e. glat > +55 degrees or glat < -55 degrees
150 # or, set magnetic grid variation to marker value
151 gv = DEFAULT_MAG_VAR
152 if math.fabs(glat) >= ABS_LAT_DD_ANY_ARCITC:
153 if glat > 0.0:
154 gv = dec - glon if glon >= 0.0 else dec + math.fabs(glon)
155 if glat < 0.0:
156 gv = dec + glon if glon >= 0.0 else dec - math.fabs(glon)
157 while gv > +180.0:
158 gv = gv - 360.0
159 while gv < -180.0: 159 ↛ 160line 159 didn't jump to line 160, because the condition on line 159 was never true
160 gv = gv + 360.0
162 return MagneticVector(
163 dec=dec, dip=dip, ti=ti, bh=bh, bx=bx, by=by, bz=bz, lat=lat_dd, lon=lon_dd, alt=alt_ft, time=time
164 )
167def calc(options: argparse.Namespace) -> int:
168 """Drive the calculation."""
169 quiet = options.quiet
170 verbose = options.verbose
171 if quiet:
172 logging.getLogger().setLevel(logging.ERROR)
173 elif verbose:
174 logging.getLogger().setLevel(logging.DEBUG)
175 lat_dd = options.lat_dd
176 lon_dd = options.lon_dd
177 alt_ft = options.alt_ft
178 date = options.date_in
179 start_time = dti.datetime.now(tz=dti.timezone.utc)
180 log.info(f'starting calculation at {date} for {alt_ft=}, {lat_dd=}, and {lon_dd=}')
181 model = Model(date.year)
182 log.info(f'- fetched model with id {model.coefficients.model["model_id"]} ')
183 data = model.at(lat_dd=lat_dd, lon_dd=lon_dd, alt_ft=alt_ft, date=date)
184 log.info(f'- result for {data.time=}, {data.alt=}, {data.lat=}, and {data.lon=}:')
185 log.info(f' Declination D = {data.dec}, inclination I = {data.dip} in degrees')
186 log.info(f' North X = {data.bx}, East Y = {data.by}, Down Z = {data.bz} component in nT')
187 log.info(f' Horizontal H = {data.bh}, total T = {data.ti} intensity in nT')
188 end_time = dti.datetime.now(tz=dti.timezone.utc)
189 log.info(f'calculation complete after {(end_time - start_time).total_seconds()} seconds')
190 print(f'INP(time={data.time}, alt={data.alt}, lat={data.lat}, long={data.lon})')
191 print(f'OUT(D={data.dec}, I={data.dip}, X={data.bx}, Y={data.by}, Z={data.bz}, H={data.bh}, T={data.ti})')
192 return 0