diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py
index dbce845d062a71634dea7b06e53e73a81a13d3c3..72708859b8382b38829db4c34d3cc926417ac739 100644
--- a/geomagio/residual/Calculation.py
+++ b/geomagio/residual/Calculation.py
@@ -37,24 +37,44 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading:
hemisphere=reading.hemisphere, measurements=reading.measurements
)
corrected_f = f + reading.pier_correction
- # calculate absolutes
- absoluteH, absoluteZ = calculate_HZ_absolutes(
+
+ # calculate absolutes,
+ (absoluteH, absoluteZ) = calculate_HZ_absolutes(
corrected_f=corrected_f,
inclination=inclination,
mean=i_mean,
reference=reference,
)
+
+ horizontal_component = calculate_horizontal_component(
+ inclination=inclination, corrected_f=corrected_f
+ )
+
+ vertical_component = calculate_vertical_component(
+ inclination=inclination, corrected_f=corrected_f
+ )
+
absoluteD, meridian = calculate_D_absolute(
azimuth=reading.azimuth,
h_baseline=absoluteH.baseline,
measurements=reading.measurements,
reference=reference,
)
- # populate diagnostics object with averaged measurements
- diagnostics = Diagnostics(
- inclination=inclination,
- meridian=meridian,
+
+ d_computed = calculate_D_computed(
+ measurements=reading.measurements,
+ h_baseline=absoluteH.baseline,
)
+
+ r_computed = calculate_R_computed(
+ measurements=reading.measurements,
+ h_baseline=absoluteH.baseline,
+ )
+
+ magnetic_south_meridian = calculate_magnetic_south_meridian(
+ measurements=reading.measurements,
+ )
+
# calculate scale
scale_value = None
scale_measurements = reading[mt.NORTH_DOWN_SCALE]
@@ -64,14 +84,59 @@ def calculate(reading: Reading, adjust_reference: bool = True) -> Reading:
inclination=inclination,
measurements=scale_measurements,
)
+
+ mean_mark = calculate_mean_mark(
+ measurements=reading.measurements,
+ )
+
+ magnetic_azimuth_of_mark = calculate_magnetic_azimuth_of_mark(
+ mean_mark=mean_mark, magnetic_south_meridian=magnetic_south_meridian
+ )
+
+ ordinate_f = f
+
+ ordinate_d = calculate_ordinate_d(
+ absoluteD_absolute=absoluteD.absolute,
+ absoluteD_baseline=absoluteD.baseline,
+ )
+
+ ordinate_h = calculate_ordinate_h(
+ absoluteH_absolute=absoluteH.absolute,
+ absoluteH_baseline=absoluteH.baseline,
+ ordinate_d=ordinate_d,
+ )
+
+ ordinate_z = calculate_ordinate_z(
+ absoluteZ_absolute=absoluteZ.absolute,
+ absoluteZ_baseline=absoluteZ.baseline,
+ )
+
+ # populate diagnostics object with averaged measurements
+ diagnostics = Diagnostics(
+ corrected_f=corrected_f,
+ d_computed=d_computed,
+ horizontal_component=horizontal_component,
+ inclination=inclination,
+ magnetic_azimuth_of_mark=magnetic_azimuth_of_mark,
+ magnetic_south_meridian=magnetic_south_meridian,
+ mean_mark=mean_mark,
+ ordinate_d=ordinate_d,
+ ordinate_f=ordinate_f,
+ ordinate_h=ordinate_h,
+ ordinate_z=ordinate_z,
+ r_computed=r_computed,
+ vertical_component=vertical_component,
+ )
+
# create new reading object
calculated = Reading(
absolutes=[absoluteD, absoluteH, absoluteZ],
- scale_value=scale_value,
diagnostics=diagnostics,
+ scale_value=scale_value,
# copy other attributes
- **reading.dict(exclude={"absolutes", "scale_value", "diagnostics"}),
+ **reading.dict(exclude={"absolutes", "diagnostics", "scale_value"}),
)
+
return calculated
@@ -154,8 +219,8 @@ def calculate_HZ_absolutes(
corrected_f: float,
mean: AverageMeasurement,
reference: Measurement,
-) -> Tuple[Absolute, Absolute]:
- """Calculate H and Z absolutes.
+) -> Tuple[Absolute, Absolute, float, float]:
+ """Calculate H and Z absolutes
Parameters
----------
@@ -292,3 +357,235 @@ def calculate_scale_value(
residual_change = m2.residual - m1.residual
scale_value = corrected_f * field_change / np.abs(residual_change)
return scale_value
+
+
+def calculate_magnetic_south_meridian(measurements: List[Measurement]) -> float:
+ """Calculate magnetic_south_meridian.
+
+ Parameters
+ ----------
+ measurements: measurements, should have type WEST_DOWN, EAST_DOWN, WEST_UP, EAST_UP.
+
+ Returns
+ -------
+ Calculated magnetic_south_meridian.
+ """
+ # declination measurements
+ declination_measurements = [
+ average_measurement(measurements, [t]) for t in DECLINATION_TYPES
+ ]
+
+ magnetic_south_meridian = np.average([m.angle for m in declination_measurements])
+ return magnetic_south_meridian
+
+
+def calculate_D_computed(
+ measurements: List[Measurement],
+ h_baseline: float,
+) -> float:
+ """Calculate d_computed.
+
+ Parameters
+ ----------
+ measurements: measurements,
+ h_baseline
+
+ Returns
+ -------
+ Calculated d_computed.
+ """
+ # declination measurements
+ declination_measurements = [
+ average_measurement(measurements, [t]) for t in DECLINATION_TYPES
+ ]
+
+ d_computed = np.average(
+ [
+ -1 * np.degrees(np.arctan(m.e / (m.h + h_baseline)))
+ for m in declination_measurements
+ ]
+ )
+
+ return d_computed
+
+
+def calculate_R_computed(
+ measurements: List[Measurement],
+ h_baseline: float,
+) -> float:
+ """Calculate r_computed.
+
+ Parameters
+ ----------
+ measurements: measurements,
+ h_baseline
+
+ Returns
+ -------
+ Calculated r_computed.
+ """
+ # declination measurements
+ declination_measurements = [
+ average_measurement(measurements, [t]) for t in DECLINATION_TYPES
+ ]
+
+ r_computed = np.average(
+ [
+ np.degrees(
+ m.measurement_type.meridian
+ * (np.arcsin(m.residual / np.sqrt((m.h + h_baseline) ** 2 + m.e**2)))
+ )
+ for m in declination_measurements
+ ]
+ )
+
+ return r_computed
+
+
+def calculate_mean_mark(
+ measurements: List[Measurement],
+) -> float:
+ """Calculate mean_mark.
+
+ Parameters
+ ----------
+ measurements: measurements,
+
+
+ Returns
+ -------
+ Calculated mean_mark angle.
+ """
+
+ # average mark
+ average_mark = average_measurement(measurements, MARK_TYPES)
+
+ return average_mark.angle
+
+
+def calculate_magnetic_azimuth_of_mark(
+ mean_mark: float,
+ magnetic_south_meridian: float,
+) -> float:
+ """Calculate magnetic_azimuth_of_mark.
+
+ Parameters
+ ----------
+ mean_mark: float,
+ magnetic_south_meridian:float,
+
+
+ Returns
+ -------
+ Calculated magnetic_azimuth_of_mark
+ """
+
+ magnetic_azimuth_of_mark = mean_mark - magnetic_south_meridian
+
+ return magnetic_azimuth_of_mark
+
+
+def calculate_ordinate_d(
+ absoluteD_absolute: float,
+ absoluteD_baseline: float,
+) -> float:
+ """Calculate ordinate_d.
+
+ Parameters
+ ----------
+ absoluteD_absolute: float,
+ absoluteD_baseline:float,
+
+
+ Returns
+ -------
+ Calculated ordinate_d
+ """
+
+ ordinate_d = absoluteD_absolute - absoluteD_baseline
+
+ return ordinate_d
+
+
+def calculate_ordinate_z(
+ absoluteZ_absolute: float,
+ absoluteZ_baseline: float,
+) -> float:
+ """Calculate ordinate_z.
+
+ Parameters
+ ----------
+ absoluteZ_absolute: float,
+ absoluteZ_baseline:float,
+
+
+ Returns
+ -------
+ Calculated ordinate_z
+ """
+
+ ordinate_z = absoluteZ_absolute - absoluteZ_baseline
+
+ return ordinate_z
+
+
+def calculate_ordinate_h(
+ absoluteH_absolute: float,
+ absoluteH_baseline: float,
+ ordinate_d: float,
+) -> float:
+ """Calculate ordinate_z.
+
+ Parameters
+ ----------
+ absoluteH_absolute: float,
+ absoluteH_baseline:float,
+ ordinate_d:float,
+
+ Returns
+ -------
+ Calculated ordinate_h
+ """
+ ordinate_e = absoluteH_absolute * np.sin(np.radians(ordinate_d))
+
+ ordinate_h = np.sqrt(absoluteH_absolute**2 - ordinate_e**2) - absoluteH_baseline
+
+ return ordinate_h
+
+
+def calculate_horizontal_component(inclination: float, corrected_f: float) -> float:
+ """Calculate horizontal_component.
+
+ Parameters
+ ----------
+ inclination: calculated inclination.
+ corrected_f: calculated f with pier correction.
+
+ Returns
+ -------
+ Calculated horizontal_component
+ """
+
+ inclination_radians = np.radians(inclination)
+ horizontal_component = corrected_f * np.cos(inclination_radians)
+
+ return horizontal_component
+
+
+def calculate_vertical_component(inclination: float, corrected_f: float) -> float:
+ """Calculate vertical_component.
+
+ Parameters
+ ----------
+ inclination: calculated inclination.
+ corrected_f: calculated f with pier correction.
+
+ Returns
+ -------
+ Calculated vertical_component
+ """
+
+ inclination_radians = np.radians(inclination)
+ vertical_component = corrected_f * np.sin(inclination_radians)
+
+ return vertical_component
diff --git a/geomagio/residual/Diagnostics.py b/geomagio/residual/Diagnostics.py
index 42f8572a7372b9c03d8b90406db5d9235509ebd3..7c75de63106022b30bead8ffea2653549d916c3a 100644
--- a/geomagio/residual/Diagnostics.py
+++ b/geomagio/residual/Diagnostics.py
@@ -8,9 +8,32 @@ class Diagnostics(BaseModel):
Attributes
----------
- inclination: Average of inclination measurements
- meridian: Calculated meridian value
+ corrected_f: Optional[float] = None
+ d_computed: Optional[float] = None
+ horizontal_component: Optional[float] = None
+ inclination: Optional[float] = None
+ magnetic_azimuth_of_mark: Optional[float] = None
+ magnetic_south_meridian: Optional[float] = None
+ mean_mark: Optional[float] = None
+ meridian: Optional[float] = None
+ ordinate_d: Optional[float] = None
+ ordinate_f: Optional[float] = None
+ ordinate_h: Optional[float] = None
+ ordinate_z: Optional[float] = None
+ r_computed: Optional[float] = None
+ vertical_component: Optional[float] = None
"""
- inclination: float
- meridian: Optional[float] = None
+ corrected_f: Optional[float] = None
+ d_computed: Optional[float] = None
+ horizontal_component: Optional[float] = None
+ inclination: Optional[float] = None
+ magnetic_azimuth_of_mark: Optional[float] = None
+ magnetic_south_meridian: Optional[float] = None
+ mean_mark: Optional[float] = None
+ ordinate_d: Optional[float] = None
+ ordinate_f: Optional[float] = None
+ ordinate_h: Optional[float] = None
+ ordinate_z: Optional[float] = None
+ r_computed: Optional[float] = None
+ vertical_component: Optional[float] = None
diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py
index 101ba4b4b7ed3c7e36024e6c8234b03cb1aea55c..752ffc4b0c4d778072584b2de468785444f61f52 100644
--- a/geomagio/residual/Reading.py
+++ b/geomagio/residual/Reading.py
@@ -33,12 +33,12 @@ class Reading(BaseModel):
absolutes: List[Absolute] = []
azimuth: float = 0
+ diagnostics: Diagnostics = None
hemisphere: Literal[-1, 1] = 1
measurements: List[Measurement] = []
metadata: Dict = {}
pier_correction: float = 0
scale_value: float = None
- diagnostics: Diagnostics = None
def __getitem__(self, measurement_type: MeasurementType):
"""Provide access to measurements by type.