diff --git a/README.md b/README.md index 585c6d1..97141aa 100644 --- a/README.md +++ b/README.md @@ -56,7 +56,7 @@ from py_ballisticcalc.unit import * ## Simple Zero ```python -# Establish 100-yard zero for a standard .308, G7 BC=0.22, muzzle velocity 2600fps +# Establish 100-yard zero for a standard .308, G7 bc=0.22, muzzle velocity 2600fps zero = Shot(weapon=Weapon(sight_height=2), ammo=Ammo(DragModel(0.22, TableG7), mv=Velocity.FPS(2600))) calc = Calculator() zero_distance = Distance.Yard(100) diff --git a/py_ballisticcalc/conditions.py b/py_ballisticcalc/conditions.py index 5204ee7..8fd7526 100644 --- a/py_ballisticcalc/conditions.py +++ b/py_ballisticcalc/conditions.py @@ -3,14 +3,14 @@ import math from dataclasses import dataclass, field +from .munition import Weapon, Ammo from .settings import Settings as Set from .unit import Distance, Velocity, Temperature, Pressure, TypedUnits, Angular -from .munition import Weapon, Ammo __all__ = ('Atmo', 'Wind', 'Shot') -cStandardHumidity: float = 0.0 # Relative Humidity -cPressureExponent: float = 5.255876 # =g*M/R*L +cStandardHumidity: float = 0.0 # Relative Humidity +cPressureExponent: float = 5.255876 # =g*M/R*L cA0: float = 1.24871 cA1: float = 0.0988438 cA2: float = 0.00152907 @@ -18,20 +18,21 @@ cA4: float = 4.21329e-07 cA5: float = 3.342e-04 # ISA, metric units: (https://www.engineeringtoolbox.com/international-standard-atmosphere-d_985.html) -cDegreesCtoK: float = 273.15 # °K = °C + 273.15 -cStandardTemperatureC: float = 15.0 # °C -cLapseRateMetric: float = -6.5e-03 # Lapse Rate, °C/m +cDegreesCtoK: float = 273.15 # °K = °C + 273.15 +cStandardTemperatureC: float = 15.0 # °C +cLapseRateMetric: float = -6.5e-03 # Lapse Rate, °C/m cStandardPressureMetric: float = 1013.25 # hPa -cSpeedOfSoundMetric: float = 331.3 # Mach1 in m/s = cSpeedOfSound * sqrt(°K) -cStandardDensityMetric: float = 1.2250 # kg/m^3 -cDensityImperialToMetric: float = 16.0185 # lb/ft^3 to kg/m^3 +cSpeedOfSoundMetric: float = 331.3 # Mach1 in m/s = cSpeedOfSound * sqrt(°K) +cStandardDensityMetric: float = 1.2250 # kg/m^3 +cDensityImperialToMetric: float = 16.0185 # lb/ft^3 to kg/m^3 # ICAO standard atmosphere: -cDegreesFtoR: float = 459.67 # °R = °F + 459.67 -cStandardTemperatureF: float = 59.0 # °F +cDegreesFtoR: float = 459.67 # °R = °F + 459.67 +cStandardTemperatureF: float = 59.0 # °F cLapseRateImperial: float = -3.56616e-03 # Lapse rate, °F/ft -cStandardPressure: float = 29.92 # InHg -cSpeedOfSoundImperial: float = 49.0223 # Mach1 in fps = cSpeedOfSound * sqrt(°R) -cStandardDensity: float = 0.076474 # lb/ft^3 +cStandardPressure: float = 29.92 # InHg +cSpeedOfSoundImperial: float = 49.0223 # Mach1 in fps = cSpeedOfSound * sqrt(°R) +cStandardDensity: float = 0.076474 # lb/ft^3 + @dataclass class Atmo(TypedUnits): # pylint: disable=too-many-instance-attributes @@ -39,7 +40,7 @@ class Atmo(TypedUnits): # pylint: disable=too-many-instance-attributes altitude: [float, Distance] = field(default_factory=lambda: Set.Units.distance) pressure: [float, Pressure] = field(default_factory=lambda: Set.Units.pressure) temperature: [float, Temperature] = field(default_factory=lambda: Set.Units.temperature) - humidity: float = 0.0 # Relative humidity [0% to 100%] + humidity: float = 0.0 # Relative humidity [0% to 100%] density_ratio: float = field(init=False) # Density / cStandardDensity mach: Velocity = field(init=False) # Mach 1 in reference atmosphere _mach1: float = field(init=False) # Mach 1 in reference atmosphere in fps @@ -58,7 +59,7 @@ def __post_init__(self): if not self.temperature: self.temperature = Atmo.standard_temperature(self.altitude) if not self.pressure: - self.pressure = Atmo.standard_pressure(self.altitude) + self.pressure = Atmo.standard_pressure(self.altitude) self._t0 = self.temperature >> Temperature.Fahrenheit self._p0 = self.pressure >> Pressure.InHg @@ -72,29 +73,29 @@ def __post_init__(self): def standard_temperature(altitude: Distance) -> Temperature: """ICAO standard temperature for altitude""" return Temperature.Fahrenheit(cStandardTemperatureF - + (altitude >> Distance.Foot) * cLapseRateImperial) + + (altitude >> Distance.Foot) * cLapseRateImperial) @staticmethod def standard_pressure(altitude: Distance) -> Pressure: """ICAO standard pressure for altitude""" return Pressure.InHg(0.02953 - * math.pow(3.73145 - (2.56555e-05) * (altitude >> Distance.Foot), - cPressureExponent) - ) + * math.pow(3.73145 - 2.56555e-05 * (altitude >> Distance.Foot), + cPressureExponent) + ) # # Metric formula # Pressure.hPa(cStandardPressureMetric # * math.pow(1 - cLapseRateMetric * (altitude >> Distance.Meter) / (cStandardTemperatureC + cDegreesCtoK), # cPressureExponent)) @staticmethod - def standard(altitude: [float, Distance] = 0, temperature: Temperature=None): + def standard(altitude: [float, Distance] = 0, temperature: Temperature = None): """Creates standard ICAO atmosphere at given altitude. If temperature not specified uses standard temperature. """ return Atmo.icao(altitude, temperature) @staticmethod - def icao(altitude: [float, Distance] = 0, temperature: Temperature=None): + def icao(altitude: [float, Distance] = 0, temperature: Temperature = None): """Creates standard ICAO atmosphere at given altitude. If temperature not specified uses standard temperature. """ @@ -118,7 +119,7 @@ def machF(fahrenheit: float) -> float: @staticmethod def machC(celsius: float) -> float: """:return: Mach 1 in m/s for Celsius temperature""" - return math.sqrt(1+celsius / cDegreesCtoK) * cSpeedOfSoundMetric + return math.sqrt(1 + celsius / cDegreesCtoK) * cSpeedOfSoundMetric @staticmethod def air_density(t: Temperature, p: Pressure, humidity: float) -> float: @@ -128,11 +129,11 @@ def air_density(t: Temperature, p: Pressure, humidity: float) -> float: tC = t >> Temperature.Celsius pM = (p >> Pressure.hPa) * 100 # Pressure in Pascals # Tetens approximation to saturation vapor pressure: - psat = 6.1078*math.pow(10, 17.27 * tC / (tC + 237.3)) - pv = humidity*psat # Pressure of water vapor in Pascals - pd = pM - pv # Partial pressure of dry air in Pascals + psat = 6.1078 * math.pow(10, 17.27 * tC / (tC + 237.3)) + pv = humidity * psat # Pressure of water vapor in Pascals + pd = pM - pv # Partial pressure of dry air in Pascals # Density in metric units kg/m^3 - density = (pd*0.0289652 + pv*0.018016)/(8.31446 *(tC + cDegreesCtoK)) + density = (pd * 0.0289652 + pv * 0.018016) / (8.31446 * (tC + cDegreesCtoK)) return density / cDensityImperialToMetric @property @@ -181,7 +182,7 @@ def get_density_factor_and_mach_for_altitude(self, altitude: float): mach = self._mach1 else: # https://en.wikipedia.org/wiki/Density_of_air#Exponential_approximation - density_ratio = math.exp(-altitude/34112.0) + density_ratio = math.exp(-altitude / 34112.0) t = self.temperature_at_altitude(altitude) mach = Atmo.machF(t) return density_ratio, mach @@ -229,6 +230,7 @@ class Shot(TypedUnits): ammo: Ammo = field(default=None) atmo: Atmo = field(default=None) winds: list[Wind] = field(default=None) + # NOTE: Calculator assumes that winds are sorted by Wind.until_distance (ascending) @property @@ -237,14 +239,14 @@ def barrel_elevation(self) -> Angular: return Angular.Radian((self.look_angle >> Angular.Radian) + math.cos(self.cant_angle >> Angular.Radian) * ((self.weapon.zero_elevation >> Angular.Radian) - + (self.relative_angle >> Angular.Radian))) + + (self.relative_angle >> Angular.Radian))) @property def barrel_azimuth(self) -> Angular: """Horizontal angle of barrel relative to sight line""" return Angular.Radian(math.sin(self.cant_angle >> Angular.Radian) * ((self.weapon.zero_elevation >> Angular.Radian) - + (self.relative_angle >> Angular.Radian))) + + (self.relative_angle >> Angular.Radian))) def __post_init__(self): if not self.look_angle: diff --git a/py_ballisticcalc/drag_model.py b/py_ballisticcalc/drag_model.py index f6bfb74..aa35fd6 100644 --- a/py_ballisticcalc/drag_model.py +++ b/py_ballisticcalc/drag_model.py @@ -1,26 +1,27 @@ """Drag model of projectile""" import math -import typing from dataclasses import dataclass, field +from typing import Mapping, Union from .settings import Settings as Set from .unit import Weight, Distance, Velocity -#from .drag_tables import DragTablesSet __all__ = ('DragModel', 'DragDataPoint', 'BCpoint', 'DragModelMultiBC') cSpeedOfSoundMetric = 340.0 # Speed of sound in standard atmosphere, in m/s + @dataclass class DragDataPoint: """Drag coefficient at Mach number""" Mach: float # Velocity in Mach units - CD: float # Drag coefficient + CD: float # Drag coefficient + @dataclass(order=True) class BCpoint: - """For multi-BC drag models, designed to sort by Mach ascending""" + """For multi-bc drag models, designed to sort by Mach ascending""" BC: float = field(compare=False) # Ballistic Coefficient at the given Mach number Mach: float = field(default=-1, compare=True) # Velocity in Mach units # Velocity only referenced if Mach number not supplied @@ -31,11 +32,12 @@ def __post_init__(self): if self.Mach < 0: self.Mach = (self.V >> Velocity.MPS) / cSpeedOfSoundMetric if self.BC <= 0: - raise ValueError('BC must be positive') + raise ValueError('bc must be positive') + class DragModel: """ - :param BC: Ballistic Coefficient of bullet = weight / diameter^2 / i, + :param bc: Ballistic Coefficient of bullet = weight / diameter^2 / i, where weight is in pounds, diameter is in inches, and i is the bullet's form factor relative to the selected drag model. :param drag_table: If passed as List of {Mach, CD} dictionaries, this @@ -45,23 +47,24 @@ class DragModel: :param length: Bullet length in inches NOTE: .weight, .diameter, .length are only relevant for computing spin drift """ - def __init__(self, BC: float, - drag_table: typing.Iterable, - weight: [float, Weight]=0, - diameter: [float, Distance]=0, - length: [float, Distance]=0): + + def __init__(self, bc: float, + drag_table: list[Mapping[float, float]], + weight: [float, Weight] = 0, + diameter: [float, Distance] = 0, + length: [float, Distance] = 0): error = '' if len(drag_table) <= 0: error = 'Received empty drag table' - elif BC <= 0: + elif bc <= 0: error = 'Ballistic coefficient must be positive' if error: raise ValueError(error) - self.drag_table = drag_table if isinstance(drag_table[0], DragDataPoint)\ - else make_data_points(drag_table) # Convert from list of dicts to list of DragDataPoints + self.drag_table = drag_table if isinstance(drag_table[0], DragDataPoint) \ + else make_data_points(drag_table) # Convert from list of dicts to list of DragDataPoints - self.BC = BC + self.BC = bc self.length = Set.Units.length(length) self.weight = Set.Units.weight(weight) self.diameter = Set.Units.diameter(diameter) @@ -70,7 +73,7 @@ def __init__(self, BC: float, self.form_factor = self._get_form_factor(self.BC) def __repr__(self) -> str: - return f"DragModel(BC={self.BC}, wgt={self.weight}, dia={self.diameter}, len={self.length})" + return f"DragModel(bc={self.BC}, wgt={self.weight}, dia={self.diameter}, len={self.length})" def _get_form_factor(self, bc: float) -> float: return self.sectional_density / bc @@ -81,7 +84,7 @@ def _get_sectional_density(self) -> float: return sectional_density(w, d) -def make_data_points(drag_table: typing.Iterable) -> list[BCpoint]: +def make_data_points(drag_table: list[Mapping[float, float]]) -> list[DragDataPoint]: """Convert drag table from list of dictionaries to list of DragDataPoints""" return [DragDataPoint(point['Mach'], point['CD']) for point in drag_table] @@ -95,14 +98,17 @@ def sectional_density(weight: float, diameter: float) -> float: return weight / math.pow(diameter, 2) / 7000 -def DragModelMultiBC(BCpoints: list[BCpoint], drag_table: typing.Iterable, - weight: [float, Weight]=0, - diameter: [float, Distance]=0, - length: [float, Distance]=0) -> DragModel: - """Compute a drag model based on multiple BCs. - If weight and diameter are provided then we set BC=sectional density. - Otherwise we set BC=1 and the drag_table contains final drag terms. - :param BC: +def DragModelMultiBC(bc_points: list[BCpoint], + drag_table: list[Mapping[float, float]], + weight: [float, Weight] = 0, + diameter: [float, Distance] = 0, + length: [float, Distance] = 0) -> DragModel: + """ + Compute a drag model based on multiple BCs. + If weight and diameter are provided then we set bc=sectional density. + Otherwise, we set bc=1 and the drag_table contains final drag terms. + :param bc_points: + :param drag_table: list of dicts containing drag table data :param weight: Bullet weight in grains :param diameter: Bullet diameter in inches :param length: Bullet length in inches @@ -114,19 +120,21 @@ def DragModelMultiBC(BCpoints: list[BCpoint], drag_table: typing.Iterable, else: BC = 1.0 - drag_table = drag_table if isinstance(drag_table[0], DragDataPoint)\ - else make_data_points(drag_table) # Convert from list of dicts to list of DragDataPoints + drag_table = drag_table if isinstance(drag_table[0], DragDataPoint) \ + else make_data_points(drag_table) # Convert from list of dicts to list of DragDataPoints - BCpoints = sorted(BCpoints) # Make sure BCpoints are sorted for linear interpolation + bc_points = sorted(bc_points) # Make sure bc_points are sorted for linear interpolation BCinterp = linear_interpolation([x.Mach for x in drag_table], - [x.Mach for x in BCpoints], - [x.BC / BC for x in BCpoints]) + [x.Mach for x in bc_points], + [x.BC / BC for x in bc_points]) for i in range(len(drag_table)): drag_table[i].CD = drag_table[i].CD / BCinterp[i] return DragModel(BC, drag_table, weight, diameter, length) -def linear_interpolation(x: list[float], xp: list[float], yp: list[float]) -> list[float]: +def linear_interpolation(x: Union[list[float], tuple[float]], + xp: Union[list[float], tuple[float]], + yp: Union[list[float], tuple[float]]) -> Union[list[float], tuple[float]]: """Piecewise linear interpolation :param x: List of points for which we want interpolated values :param xp: List of existing points (x coordinate), *sorted in ascending order* @@ -135,7 +143,7 @@ def linear_interpolation(x: list[float], xp: list[float], yp: list[float]) -> li """ if len(xp) != len(yp): raise ValueError("xp and yp lists must have same length") - + y = [] for xi in x: @@ -148,7 +156,7 @@ def linear_interpolation(x: list[float], xp: list[float], yp: list[float]) -> li left, right = 0, len(xp) - 1 while left < right: mid = (left + right) // 2 - if xp[mid] <= xi < xp[mid+1]: + if xp[mid] <= xi < xp[mid + 1]: slope = (yp[mid + 1] - yp[mid]) / (xp[mid + 1] - xp[mid]) y.append(yp[mid] + slope * (xi - xp[mid])) # Interpolated value for xi break diff --git a/py_ballisticcalc/drag_tables.py b/py_ballisticcalc/drag_tables.py index da0f297..fb6115c 100644 --- a/py_ballisticcalc/drag_tables.py +++ b/py_ballisticcalc/drag_tables.py @@ -1,5 +1,6 @@ """Templates of the common used drag tables""" + TableG1 = [ {'Mach': 0.00, 'CD': 0.2629}, {'Mach': 0.05, 'CD': 0.2558}, @@ -666,5 +667,3 @@ {'Mach': 3.95, 'CD': 0.9295}, {'Mach': 4.00, 'CD': 0.9280}, ] - -DragTablesSet = [value for key, value in globals().items() if key.startswith("Table")] diff --git a/py_ballisticcalc/trajectory_calc.py b/py_ballisticcalc/trajectory_calc.py index cf34cfa..d5e93eb 100644 --- a/py_ballisticcalc/trajectory_calc.py +++ b/py_ballisticcalc/trajectory_calc.py @@ -312,7 +312,7 @@ def drag_by_mach(self, mach: float) -> float: cStandardDensity of Air = 0.076474 lb/ft^3 S is cross-section = d^2 pi/4, where d is bullet diameter in inches m is bullet mass in pounds - BC contains m/d^2 in units lb/in^2, which we multiply by 144 to convert to lb/ft^2 + bc contains m/d^2 in units lb/in^2, which we multiply by 144 to convert to lb/ft^2 Thus: The magic constant found here = StandardDensity * pi / (4 * 2 * 144) :return: Drag coefficient at the given mach number """ diff --git a/py_ballisticcalc_exts/py_ballisticcalc_exts/trajectory_calc.pyx b/py_ballisticcalc_exts/py_ballisticcalc_exts/trajectory_calc.pyx index 98f15eb..c018443 100644 --- a/py_ballisticcalc_exts/py_ballisticcalc_exts/trajectory_calc.pyx +++ b/py_ballisticcalc_exts/py_ballisticcalc_exts/trajectory_calc.pyx @@ -320,7 +320,7 @@ cdef class TrajectoryCalc: cStandardDensity of Air = 0.076474 lb/ft^3 S is cross-section = d^2 pi/4, where d is bullet diameter in inches m is bullet mass in pounds - BC contains m/d^2 in units lb/in^2, which we multiply by 144 to convert to lb/ft^2 + bc contains m/d^2 in units lb/in^2, which we multiply by 144 to convert to lb/ft^2 Thus: The magic constant found here = StandardDensity * pi / (4 * 2 * 144) """ cdef double cd = calculate_by_curve(self._table_data, self._curve, mach) diff --git a/tests/test_computer.py b/tests/test_computer.py index de06f65..d3d9a7c 100644 --- a/tests/test_computer.py +++ b/tests/test_computer.py @@ -142,7 +142,7 @@ def test_pressure(self): #region Ammo def test_ammo_drag(self): - """Increasing ballistic coefficient (BC) should decrease drop""" + """Increasing ballistic coefficient (bc) should decrease drop""" tdm = DragModel(self.dm.BC+0.5, self.dm.drag_table, self.dm.weight, self.dm.diameter, self.dm.length) slick = Ammo(tdm, self.ammo.mv) shot = Shot(weapon=self.weapon, ammo=slick, atmo=self.atmosphere) diff --git a/tests/test_mbc.py b/tests/test_mbc.py index 9cc41ad..1849a6f 100644 --- a/tests/test_mbc.py +++ b/tests/test_mbc.py @@ -1,4 +1,4 @@ -"Unit tests of multiple-BC drag models" +"Unit tests of multiple-bc drag models" import unittest from py_ballisticcalc import * @@ -18,7 +18,7 @@ def setUp(self) -> None: self.baseline_trajectory = self.calc.fire(shot=self.baseline_shot, trajectory_range=self.range, trajectory_step=self.step).trajectory def test_mbc1(self): - "We should get the same trajectory whether we give single BC or use multi-BC with single value" + "We should get the same trajectory whether we give single bc or use multi-bc with single value" dm_multi = DragModelMultiBC([BCpoint(.22, V=Velocity.FPS(2500)), BCpoint(.22, V=Velocity.FPS(1500)), BCpoint(BC=.22, Mach=3)], TableG7) multi_shot = Shot(weapon=self.weapon, ammo=Ammo(dm_multi, self.ammo.mv)) multi_trajectory = self.calc.fire(shot=multi_shot, trajectory_range=self.range, trajectory_step=self.step).trajectory @@ -26,7 +26,7 @@ def test_mbc1(self): self.assertEqual(multi_trajectory[i].formatted(), self.baseline_trajectory[i].formatted()) def test_mbc2(self): - "Setting different BC above muzzle velocity should have no effect" + "Setting different bc above muzzle velocity should have no effect" dm_multi = DragModelMultiBC([BCpoint(.22, V=Velocity.FPS(2700)), BCpoint(.5, V=Velocity.FPS(3500))], TableG7) multi_shot = Shot(weapon=self.weapon, ammo=Ammo(dm_multi, self.ammo.mv)) multi_trajectory = self.calc.fire(shot=multi_shot, trajectory_range=self.range, trajectory_step=self.step).trajectory @@ -34,8 +34,8 @@ def test_mbc2(self): self.assertEqual(multi_trajectory[i].formatted(), self.baseline_trajectory[i].formatted()) def test_mbc3(self): - "Setting higher BC should result in higher downrange velocities" - # So here we'll boost the BC for velocities lower than the baseline's velocity at 200 yards + "Setting higher bc should result in higher downrange velocities" + # So here we'll boost the bc for velocities lower than the baseline's velocity at 200 yards dm_multi = DragModelMultiBC([BCpoint(.5, V=self.baseline_trajectory[3].velocity), BCpoint(.22, V=self.baseline_trajectory[2].velocity)], TableG7) multi_shot = Shot(weapon=self.weapon, ammo=Ammo(dm_multi, self.ammo.mv)) multi_trajectory = self.calc.fire(shot=multi_shot, trajectory_range=self.range, trajectory_step=self.step).trajectory