New Multi-BC drag model implementation:

- DragModel can take a List of BCpoints (for different Velocities) instead of a single BC, in which case it will interpolate the BC and the DragModel.drag_table will reflect the final CD (i.e., trajectory_calc will set BC=1).
- When given a DragTable in List-of-dicts form (as in drag_tables.py), DragModel converts it to more efficient List of DragDataPoints.
This code still needs refinement, and also still needs to be implemented in the Cython extensions.

py_ballisticcalc/__init__.py

@@ -11,7 +11,6 @@
 from .backend import *
 from .drag_tables import *
 from .settings import *
-from .multiple_bc import *
 from .interface import *
 from .trajectory_data import *
 from .conditions import *

py_ballisticcalc/drag_model.py

@@ -1,88 +1,110 @@
-"""definitions for a bullet's drag model calculation"""
+"""Drag model of projectile"""
 
 import typing
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 
 import math
+import numpy
 
 from .settings import Settings as Set
-from .unit import Weight, Distance
-from .drag_tables import DragTablesSet
+from .unit import Weight, Distance, Velocity
+#from .drag_tables import DragTablesSet
 
-__all__ = ('DragModel', 'make_data_points')
+__all__ = ('DragModel', 'BCpoint')
 
+cSpeedOfSoundMetric = 340.0  # Speed of sound in standard atmosphere, in m/s
 
 @dataclass
 class DragDataPoint:
-    CD: float    # Drag coefficient
     Mach: float  # Velocity in Mach units
+    CD: float    # Drag coefficient
 
     def __iter__(self):
-        yield self.CD
         yield self.Mach
+        yield self.CD
 
-    def __repr__(self):
-        return f"DragDataPoint(CD={self.CD}, Mach={self.Mach})"
-
+@dataclass(order=True)
+class BCpoint:
+    """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
+    V: Velocity = field(default_factory=lambda: Set.Units.velocity, compare=False)
+
+    def __post_init__(self):
+        # If Mach not defined then convert V using standard atmosphere
+        if self.Mach < 0:
+            self.Mach = (self.V >> Velocity.MPS) / cSpeedOfSoundMetric
+        if self.BC <= 0:
+            raise ValueError('BC must be positive')
 
 class DragModel:
     """
     :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: List of {Mach, Cd} pairs defining the standard drag model
+            where weight is in pounds, diameter is in inches, and
+            i is the bullet's form factor relative to the selected drag model.
+        Or List[BCpoint], and BC will be interpolated and applied to the .drag_table
+            (in which case self.BC = 1)
+    :param drag_table: If passed as List of {Mach, CD} dictionaries, this
+            will be converted to a List of DragDataPoints for efficiency.
     :param weight: Bullet weight in grains
     :param diameter: Bullet diameter in inches
     :param length: Bullet length in inches
     NOTE: .weight, .diameter, .length are only relevant for computing spin drift
     """
-    def __init__(self, BC: float,
+    def __init__(self, BC,
                  drag_table: typing.Iterable,
                  weight: [float, Weight]=0,
                  diameter: [float, Distance]=0,
                  length: [float, Distance]=0):
         table_len = len(drag_table)
         error = ''
         if table_len <= 0:
-            error = 'Custom drag table must be longer than 0'
-        elif BC <= 0:
+            error = 'Drag table must be longer than 0'
+        elif type(BC) is float and (BC <= 0):
             error = 'Ballistic coefficient must be greater than zero'
         if error:
             raise ValueError(error)
 
-        if drag_table in DragTablesSet:
-            self.BC = BC
-        elif table_len > 0:
+        if isinstance(drag_table[0], DragDataPoint):
+            self.drag_table = drag_table
+        else:  # Convert from list of dicts to list of DragDataPoints
+            self.drag_table = make_data_points(drag_table)
+
+        # BC is a list, so generate new drag table by interpolating the BCpoints
+        if hasattr(BC, '__getitem__'):
+            BC = sorted(BC)  # Make sure we're sorted for np.interp
+            self.BCinterp = numpy.interp([x['Mach'] for x in drag_table],
+                                         [x.Mach for x in BC],
+                                         [x.BC for x in BC])
+            for i in range(len(self.drag_table)):
+                self.drag_table[i].CD = self.drag_table[i].CD / self.BCinterp[i]
             self.BC = 1.0
         else:
-            raise ValueError('Wrong drag data')
+            self.BC = BC
 
         self.length = Set.Units.length(length)
         self.weight = Set.Units.weight(weight)
         self.diameter = Set.Units.diameter(diameter)
         if weight != 0 and diameter != 0:
             self.sectional_density = self._get_sectional_density()
             self.form_factor = self._get_form_factor(self.BC)
-        self.drag_table = drag_table
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return f"DragModel(BC={self.BC}, wgt={self.weight}, dia={self.diameter}, len={self.length})"
 
-    def _get_form_factor(self, bc: float):
+    def _get_form_factor(self, bc: float) -> float:
         return self.sectional_density / bc
 
     def _get_sectional_density(self) -> float:
         w = self.weight >> Weight.Grain
         d = self.diameter >> Distance.Inch
         return sectional_density(w, d)
 
-    @staticmethod
-    def from_mbc(mbc: 'MultiBC'):
-        return DragModel(1, mbc.cdm, mbc.weight, mbc.diameter)
-
 
 def make_data_points(drag_table: typing.Iterable) -> list:
-    return [DragDataPoint(point['CD'], point['Mach']) for point in drag_table]
+    "Convert drag table from list of dictionaries to list of DragDataPoints"
+    return [DragDataPoint(point['Mach'], point['CD']) for point in drag_table]
 
 
 def sectional_density(weight: float, diameter: float):

py_ballisticcalc/trajectory_calc.py

@@ -410,28 +410,28 @@ def calculate_ogw(bullet_weight: float, velocity: float) -> float:
 
 
 def calculate_curve(data_points) -> list[CurvePoint]:
-    """
-    :param DragTable: data_points
+    """Piecewise quadratic interpolation of drag curve
+    :param DragTable: List[{Mach, CD}] data_points in ascending Mach order
     :return: List[CurvePoints] to interpolate drag coefficient
     """
     # rate, x1, x2, x3, y1, y2, y3, a, b, c
     # curve = []
     # curve_point
     # num_points, len_data_points, len_data_range
 
-    rate = (data_points[1]['CD'] - data_points[0]['CD']) \
-           / (data_points[1]['Mach'] - data_points[0]['Mach'])
-    curve = [CurvePoint(0, rate, data_points[0]['CD'] - data_points[0]['Mach'] * rate)]
+    rate = (data_points[1].CD - data_points[0].CD) \
+           / (data_points[1].Mach - data_points[0].Mach)
+    curve = [CurvePoint(0, rate, data_points[0].CD - data_points[0].Mach * rate)]
     len_data_points = int(len(data_points))
     len_data_range = len_data_points - 1
 
     for i in range(1, len_data_range):
-        x1 = data_points[i - 1]['Mach']
-        x2 = data_points[i]['Mach']
-        x3 = data_points[i + 1]['Mach']
-        y1 = data_points[i - 1]['CD']
-        y2 = data_points[i]['CD']
-        y3 = data_points[i + 1]['CD']
+        x1 = data_points[i - 1].Mach
+        x2 = data_points[i].Mach
+        x3 = data_points[i + 1].Mach
+        y1 = data_points[i - 1].CD
+        y2 = data_points[i].CD
+        y3 = data_points[i + 1].CD
         a = ((y3 - y1) * (x2 - x1) - (y2 - y1) * (x3 - x1)) / (
                 (x3 * x3 - x1 * x1) * (x2 - x1) - (x2 * x2 - x1 * x1) * (x3 - x1))
         b = (y2 - y1 - a * (x2 * x2 - x1 * x1)) / (x2 - x1)
@@ -440,10 +440,10 @@ def calculate_curve(data_points) -> list[CurvePoint]:
         curve.append(curve_point)
 
     num_points = len_data_points
-    rate = (data_points[num_points - 1]['CD'] - data_points[num_points - 2]['CD']) / \
-           (data_points[num_points - 1]['Mach'] - data_points[num_points - 2]['Mach'])
+    rate = (data_points[num_points - 1].CD - data_points[num_points - 2].CD) / \
+           (data_points[num_points - 1].Mach - data_points[num_points - 2].Mach)
     curve_point = CurvePoint(
-        0, rate, data_points[num_points - 1]['CD'] - data_points[num_points - 2]['Mach'] * rate
+        0, rate, data_points[num_points - 1].CD - data_points[num_points - 2].Mach * rate
     )
     curve.append(curve_point)
     return curve
@@ -461,12 +461,12 @@ def calculate_by_curve(data: list, curve: list, mach: float) -> float:
 
     while mhi - mlo > 1:
         mid = int(math.floor(mhi + mlo) / 2.0)
-        if data[mid]['Mach'] < mach:
+        if data[mid].Mach < mach:
             mlo = mid
         else:
             mhi = mid
 
-    if data[mhi]['Mach'] - mach > mach - data[mlo]['Mach']:
+    if data[mhi].Mach - mach > mach - data[mlo].Mach:
         m = mlo
     else:
         m = mhi