Skip to content

Commit

Permalink
typings update
Browse files Browse the repository at this point in the history
  • Loading branch information
o-murphy committed Apr 4, 2024
1 parent f2a2cb3 commit 4040f03
Show file tree
Hide file tree
Showing 8 changed files with 83 additions and 74 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
62 changes: 32 additions & 30 deletions py_ballisticcalc/conditions.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,43 +3,44 @@
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
cA3: float = -3.07031e-06
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
"""Atmospheric conditions and density calculations"""
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
Expand All @@ -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
Expand All @@ -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.
"""
Expand All @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand Down
74 changes: 41 additions & 33 deletions py_ballisticcalc/drag_model.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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]

Expand All @@ -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
Expand All @@ -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*
Expand All @@ -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:
Expand All @@ -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
Expand Down
3 changes: 1 addition & 2 deletions py_ballisticcalc/drag_tables.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
"""Templates of the common used drag tables"""


TableG1 = [
{'Mach': 0.00, 'CD': 0.2629},
{'Mach': 0.05, 'CD': 0.2558},
Expand Down Expand Up @@ -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")]
2 changes: 1 addition & 1 deletion py_ballisticcalc/trajectory_calc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion tests/test_computer.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 4040f03

Please sign in to comment.