diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py index d912d58101507b9235b2255ab9a2b69aa98f339e..55a40d314c1ac7a795564e55a77aceaf4401ea47 100644 --- a/geomagio/residual/Calculation.py +++ b/geomagio/residual/Calculation.py @@ -59,6 +59,7 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading: h_baseline=absoluteH.baseline, measurements=reading.measurements, reference=reference, + shift=reading.absolutes[0].shift, ) d_computed = calculate_D_computed( @@ -145,6 +146,7 @@ def calculate_D_absolute( azimuth: float, h_baseline: float, reference: Measurement, + shift: float = 0, ) -> Tuple[Absolute, Diagnostics]: """Calculate D absolute. @@ -154,6 +156,7 @@ def calculate_D_absolute( azimuth: azimuth to mark in decimal degrees. h_baseline: calculated H baseline value. reference: reference measurement used to adjust absolute. + shift: usually 0 or +/-180 for Declination Returns ------- @@ -192,12 +195,26 @@ def calculate_D_absolute( for m in declination_measurements ] ) - shift = 0 - if azimuth > 180: - shift = -180 + # add subtract average mark angle from average meridian angle and add # azimuth to get the declination baseline - d_b = (meridian - average_mark.angle) + azimuth + shift + d_b = (meridian - average_mark.angle) + azimuth + + # given operational limitations, there is always an 180-degree ambiguity + # in our calculated d_b; therefore, we have to assume the magnetic field + # points north(ish), and force d_b to fall between -90 and +90 degrees + while d_b < -90: + d_b += 180 + while d_b > 90: + d_b -= 180 + + # the preceeding lead to a wrong answer if the magnetic field actually + # points south(ish); + d_b += shift + + # (for example, there are places in Antarctica where the magnetic field + # actually points south-ish) + # calculate absolute d_abs = d_b + np.degrees(np.arctan(reference.e / (reference.h + h_baseline)))