import numpy as np
from .Ordinate import Ordinate
from .Absolute import Absolute
from .Angle import to_dms
from .MeasurementType import MeasurementType as mt
from pydantic import BaseModel
class Calculate(BaseModel):
"""
Class object for performing calculations.
Contains the following:
angle: average angle across a measurement type
residual: average residual across a measurement type
hs: Multiplier for inclination claculations. +1 if measurment was taken in northern hemisphere, -1 if measurement was taken in the southern hemishpere.
ordinate: Variometer data. Ordinate object(contains a datapoint for H, E, Z, and F)
ud: Multiplier for inclination calculations. +1 if instrument is oriented upward, -1 if instrument if oriented downward.
shift: Degree shift in inclination measurements.
"""
angle: float = None
residual: float = None
ordinate: Ordinate = None
f: float = None
ud: int = None
pm: int = None
shift: int = None
def calculate(reading):
"""
Calculate/recalculate absolute from a Reading object's
ordinates, measurements, and metadata.
Outputs a list of absolutes containing baseline, absolute,
and element name. Also reutrns the scale value.
"""
# gather oridinates, measuremenets, and metadata objects from reading
metadata = reading.metadata
ordinates = reading.ordinates
ordinate_index = reading.ordinate_index()
measurements = reading.measurements
measurement_index = reading.measurement_index()
# define measurement types used to calculate inclination
inclination_types = [mt.NORTH_DOWN, mt.NORTH_UP, mt.SOUTH_DOWN, mt.SOUTH_UP]
# get ordinate values across h, e, z, and f for inclination measurement types
inclination_ordinates = [
o for o in ordinates if o.measurement_type in inclination_types
]
# get average ordinate values across h, e, z, and f
mean = average_ordinate(inclination_ordinates, None)
# calculate inclination
inclination, f = calculate_I(
measurement_index, inclination_ordinates, ordinate_index, mean, metadata,
)
# calculate absolutes
h_abs, z_abs = calculate_absolutes(f, inclination)
# calculate baselines for H and Z
h_b, z_b = calculate_baselines(h_abs, z_abs, mean, reading.pier_correction)
# calculate scale value
scale_ordinates = ordinate_index[mt.NORTH_DOWN_SCALE]
scale_measurements = measurement_index[mt.NORTH_DOWN_SCALE]
scale = calculate_scale(f, scale_ordinates, scale_measurements, inclination,)
# calculate declination absolute and baseline
d_b, d_abs = calculate_D(
ordinate_index, measurements, measurement_index, metadata["mark_azimuth"], h_b,
)
# return results as a set of Absolute objects along with the calculated scale value
resultH = Absolute(element="H", baseline=h_b, absolute=h_abs)
resultD = Absolute(element="D", baseline=d_b, absolute=d_abs)
resultZ = Absolute(element="Z", baseline=z_b, absolute=z_abs)
result = [resultH, resultD, resultZ]
return result, scale
def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
"""
Calculate inclination angles from measurements, ordinates,
average ordinates from every measurement, and metadata.
Returns inclination angle and calculated average f
"""
# get first inclination angle, assumed to be the southdown angle
Iprime = average_angle(measurements, mt.SOUTH_DOWN)
if Iprime >= 100:
Iprime -= 200
# get multiplier for hempisphere the observatory is located in
# 1 if observatory is in northern hemisphere
# -1 if observatory is in southern hemisphere
hs = metadata["hemisphere"]
# gather calculation objects for each measurement type
southdown = Calculate(
shift=-200,
ud=1,
pm=1,
angle=average_angle(measurements, mt.SOUTH_DOWN),
residual=average_residual(measurements, mt.SOUTH_DOWN),
ordinate=average_ordinate(ordinates_index, mt.SOUTH_DOWN),
)
southup = Calculate(
shift=200,
ud=-1,
pm=-1,
angle=average_angle(measurements, mt.SOUTH_UP),
residual=average_residual(measurements, mt.SOUTH_UP),
ordinate=average_ordinate(ordinates_index, mt.SOUTH_UP),
)
northup = Calculate(
shift=0,
ud=-1,
pm=1,
angle=average_angle(measurements, mt.NORTH_UP),
residual=average_residual(measurements, mt.NORTH_UP),
ordinate=average_ordinate(ordinates_index, mt.NORTH_UP),
)
northdown = Calculate(
shift=400,
ud=1,
pm=-1,
angle=average_angle(measurements, mt.NORTH_DOWN),
residual=average_residual(measurements, mt.NORTH_DOWN),
ordinate=average_ordinate(ordinates_index, mt.NORTH_DOWN),
)
# set inclination value for looping = Iprime
inclination = Iprime
# add one to inclination value to enter the loop
Inclination = inclination + 1
# loop condition
while abs(Inclination - inclination) > 0.0001:
# set temporary inclination variable to hold previous step's inclination
Inclination = inclination
# calculate f for inclination measurement types
southdown.f = calculate_f(southdown.ordinate, mean, inclination)
northdown.f = calculate_f(northdown.ordinate, mean, inclination)
southup.f = calculate_f(southup.ordinate, mean, inclination)
northup.f = calculate_f(northup.ordinate, mean, inclination)
# gather measurements into array(necessary because f were recalculated)
measurements = [southup, southdown, northup, northdown]
# average f values for inclination measurement types
f = np.average([i.f for i in measurements])
# calculation inclination for each inclination measurement type and average
inclination = np.average(
[calculate_measurement_inclination(m, hs) for m in measurements]
)
# loop exits once the difference in average inclination between steps is lower than 0.0001
return inclination, f
def calculate_D(ordinates_index, measurements, measurements_index, azimuth, h_b):
"""
Calculate declination absolute and declination baseline from
ordinates, measurements, measurement_index(dictionary), azimuth and H baseline
Returns absolute and baseline for declination.
"""
# specify mark measurement types
mark_types = [
mt.FIRST_MARK_DOWN,
mt.FIRST_MARK_UP,
mt.SECOND_MARK_DOWN,
mt.SECOND_MARK_UP,
]
# average mark angles
# note that angles are being converted to geon
average_mark = np.average(
[
convert_to_geon(m.angle)
for m in measurements
if m.measurement_type in mark_types
]
)
# add 100 if mark up is greater than mark down
# subtract 100 otherwise
if (
measurements_index["FirstMarkUp"][0].angle
< measurements_index["FirstMarkDown"][0].angle
):
average_mark += 100
else:
average_mark -= 100
# gather calculation objects for each declination measurement type
# note that the pm(plus minus) multiplier has been repurposed.
# West facing measurements have a multiplier of -1
# East facing measurements have a multipllier of 1
westdown = Calculate(
angle=average_angle(measurements_index, mt.WEST_DOWN),
residual=average_residual(measurements_index, mt.WEST_DOWN),
ordinate=average_ordinate(ordinates_index, mt.WEST_DOWN),
pm=-1,
)
westup = Calculate(
angle=average_angle(measurements_index, mt.WEST_UP),
residual=average_residual(measurements_index, mt.WEST_UP),
ordinate=average_ordinate(ordinates_index, mt.WEST_UP),
pm=-1,
)
eastdown = Calculate(
angle=average_angle(measurements_index, mt.EAST_DOWN),
residual=average_residual(measurements_index, mt.EAST_DOWN),
ordinate=average_ordinate(ordinates_index, mt.EAST_DOWN),
pm=1,
)
eastup = Calculate(
angle=average_angle(measurements_index, mt.EAST_UP),
residual=average_residual(measurements_index, mt.EAST_UP),
ordinate=average_ordinate(ordinates_index, mt.EAST_UP),
pm=1,
)
# convert azimuth to geon
azimuth = (int(azimuth / 100) + (azimuth % 100) / 60) / 0.9
# gather declination measurements into array
measurements = [westdown, westup, eastdown, eastup]
# average meridian terms calculated from each declination measurements
meridian = np.average([calculate_meridian_term(i, h_b) for i in measurements])
# add subtract average mark angle from average meridian angle and add
# azimuth(in geon) to get the declination baseline
D_b = (meridian - average_mark) + azimuth
# convert decimal baseline into dms baseline
d_b = round(D_b * 54, 2)
d_b_dms = int(d_b / 60) * 100 + ((d_b / 60) % 1) * 60
# gather first declination measurements' H ans E ordinates
wd_E_1 = ordinates_index[mt.WEST_DOWN][0].e
wd_H_1 = ordinates_index[mt.WEST_DOWN][0].h
# calculate Declination baseline
d_abs = D_b + np.arctan(wd_E_1 / (h_b + wd_H_1)) * (200 / np.pi)
d_abs = round(d_abs * 54, 1)
# convert decimal absolute into dms absolute
d_abs_dms = int(d_abs / 60) * 100 + ((d_abs / 60) % 1) * 60
return d_b_dms, d_abs_dms
def calculate_absolutes(f, inclination):
"""
Calculate absolutes for H, Z and F from computed
average f value(from inclination computations) and
calculated inclination angle.
Returns baselines for H and Z
"""
# convert inclination to radians
i = (np.pi / 200) * (inclination)
h_abs = f * np.cos(i)
z_abs = f * np.sin(i)
return h_abs, z_abs
def calculate_baselines(h_abs, z_abs, mean, pier_correction):
"""
Calculate baselines with H and Z absolutes,
average ordinates across all measurements,
and pier correction(metadata).
Returns H and Z baselines
"""
correction = pier_correction / 5
h_b = round(np.sqrt(h_abs ** 2 - mean.e ** 2) - mean.h, 1) - correction
z_b = round(z_abs - mean.z, 1) - correction
return h_b, z_b
def calculate_scale(f, ordinates, measurements, inclination):
"""
Calculate scale value from calulated f(from inclination computations),
inclination, and the measurements/ordinates taken for scaling purposes.
Such measurements usually present themselves as a set of three North Down
measurements, where the final two measuremets(and ordinates) are used for scaling.
"""
first_ord = ordinates[0]
second_ord = ordinates[1]
first_measurement = measurements[0]
second_measurement = measurements[1]
field_change = (
(
-np.sin(inclination * np.pi / 200) * (second_ord.h - first_ord.h) / f
+ np.cos(inclination * np.pi / 200) * (second_ord.z - first_ord.z) / f
)
* 200
/ np.pi
)
field_change += 0.1852
residual_change = abs(second_measurement.residual - first_measurement.residual)
scale_value = (f * field_change / residual_change) * np.pi / 200
return scale_value
def average_angle(measurements, type):
"""
Compute average angle from a dictionary of
measurements and specified measurement type.
"""
return np.average([convert_to_geon(m.angle) for m in measurements[type]])
def average_residual(measurements, type):
"""
Compute average residual from a dictionary
of measurements and specified measurement type.
"""
return np.average([m.residual for m in measurements[type]])
def average_ordinate(ordinates, type):
"""
Compute average ordinate from a dictionary
of ordinates and specified measurement type.
"""
if type is not None:
ordinates = ordinates[type]
o = Ordinate(measurement_type=type)
avgs = np.average([[o.h, o.e, o.z, o.f] for o in ordinates], axis=0,)
o.h, o.e, o.z, o.f = avgs
return o
def calculate_f(ordinate, mean, inclination):
"""
calculate f for a measurement type using a measurement's
average ordinates, average ordinate across all measurements,
and calculated inclination.
"""
# get channel means form all ordinates
# calculate f using current step's inclination angle
f = (
mean.f
+ (ordinate.h - mean.h) * np.cos(inclination * np.pi / 200)
+ (ordinate.z - mean.z) * np.sin(inclination * np.pi / 200)
+ ((ordinate.e) ** 2 - (mean.e) ** 2) / (2 * mean.f)
)
return f
def calculate_measurement_inclination(calculation, hs):
"""
Calculate a measurement's inclination value using
Calculate items' elements.
"""
return calculation.shift + calculation.pm * (
+calculation.angle
+ calculation.ud
* (hs * np.arcsin(calculation.residual / calculation.f) * 200 / np.pi)
)
def calculate_meridian_term(calculation, h_b):
"""
Calculate meridian value from a measurement type
using a Calculate object and H's baseline value.
"""
A1 = np.arcsin(
calculation.residual
/ np.sqrt((calculation.ordinate.h + h_b) ** 2 + (calculation.ordinate.e) ** 2)
)
A2 = np.arctan(calculation.ordinate.e / (calculation.ordinate.h + h_b))
A1 = (200 / np.pi) * (A1)
A2 = (200 / np.pi) * (A2)
meridian_term = calculation.angle + (calculation.pm * A1) - A2
return meridian_term
def convert_to_geon(angle, include_seconds=True):
"""
Convert angles from measurements to geon
"""
degrees = int(angle)
minutes = int((angle % 1) * 100) / 60
if include_seconds:
seconds = ((angle * 100) % 1) / 36
else:
seconds = 0
dms = (degrees + minutes + seconds) / 0.9
return dms