DragModel

def make_data_points(drag_table: DragTableDataType) -> List[DragDataPoint]:
    """Convert drag table from list of dictionaries to list of DragDataPoints.

    Handles both DragDataPoint objects and dictionaries with 'Mach' and 'CD' keys.
    Validates input format and provides clear error messages for invalid data.

    Args:
        drag_table: Either list of DragDataPoint objects or list of dictionaries
                    with 'Mach' and 'CD' keys

    Returns:
        List of DragDataPoint objects ready for use in ballistic calculations

    Raises:
        TypeError: If drag_table items are not DragDataPoint objects or valid
                   dictionaries with required keys
    """
    try:
        return [
            point if isinstance(point, DragDataPoint) else DragDataPoint(point["Mach"], point["CD"])
            for point in drag_table
        ]
    except (KeyError, TypeError) as exc:
        raise TypeError(
            "All items in drag_table must be of type DragDataPoint or dict with 'Mach' and 'CD' keys"
        ) from exc

def sectional_density(weight: float, diameter: float) -> float:
    """Calculate sectional density of a projectile.

    Args:
        weight: Projectile weight in grains
        diameter: Projectile diameter in inches

    Returns:
        Sectional density in lb/in² (pounds per square inch)

    Note:
        Formula: SD = weight / (diameter² * 7000)
        where 7000 converts grains to pounds (7000 grains = 1 pound)
    """
    return weight / math.pow(diameter, 2) / 7000

def DragModelMultiBC(
    bc_points: List[BCPoint],
    drag_table: DragTableDataType,
    weight: Union[float, Weight] = 0,
    diameter: Union[float, Distance] = 0,
    length: Union[float, Distance] = 0,
) -> DragModel:
    """Create a drag model with multiple ballistic coefficients.

    Constructs a DragModel using multiple BC measurements at different velocities,
    interpolating between them to create a more accurate drag function.

    Args:
        bc_points: List of BCPoint objects with BC measurements at specific velocities
        drag_table: Standard drag table (G1, G7, etc.) or custom drag data
        weight: Projectile weight in grains (default: 0)
        diameter: Projectile diameter in inches (default: 0)
        length: Projectile length in inches (default: 0)

    Returns:
        DragModel with interpolated drag coefficients based on multiple BCs

    Example:
        ```python
        from py_ballisticcalc.drag_tables import TableG7
        DragModelMultiBC([BCPoint(.21, V=Velocity.FPS(1500)), BCPoint(.22, V=Velocity.FPS(2500))],
                         drag_table=TableG7)
        ```

    Note:
        If weight and diameter are provided, BC is set to sectional density.
        Otherwise, BC=1 and the drag_table contains final drag terms.
        BC points are automatically sorted by Mach number for interpolation.
    """
    weight = PreferredUnits.weight(weight)
    diameter = PreferredUnits.diameter(diameter)
    if weight > 0 and diameter > 0:
        bc = sectional_density(weight >> Weight.Grain, diameter >> Distance.Inch)
    else:
        bc = 1.0

    drag_table = make_data_points(drag_table)  # Convert from list of dicts to list of DragDataPoints

    bc_points.sort(key=lambda p: p.Mach)  # Make sure bc_points are sorted for linear interpolation
    bc_interp = linear_interpolation(
        [x.Mach for x in drag_table], [x.Mach for x in bc_points], [x.BC / bc for x in bc_points]
    )

    for i, point in enumerate(drag_table):
        point.CD = point.CD / bc_interp[i]
    return DragModel(bc, drag_table, weight, diameter, length)

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]]:
    """Perform piecewise linear interpolation.

    Interpolates y-values for given x-values using linear interpolation between known data points.
    Handles extrapolation by returning boundary values for x-values outside the range of xp.

    Args:
        x: List of points for which we want interpolated values
        xp: List of existing x-coordinates (must be sorted in ascending order)
        yp: List of corresponding y-values for existing points

    Returns:
        List of interpolated y-values corresponding to input x-values

    Raises:
        AssertionError: If `xp` and `yp` lists have different lengths

    Note:
        - For x-values below `min(xp)`, returns `yp[0]`
        - For x-values above `max(xp)`, returns `yp[-1]`
        - Uses binary search for efficient interval location
    """
    assert len(xp) == len(yp), "xp and yp lists must have same length"
    # Validate xp strictly increasing to prevent zero-division and undefined intervals
    for i in range(1, len(xp)):
        if xp[i] <= xp[i - 1]:
            raise ValueError("xp must be strictly increasing with no duplicates")

    y = []
    for xi in x:
        if xi <= xp[0]:
            y.append(yp[0])
        elif xi >= xp[-1]:
            y.append(yp[-1])
        else:
            # Binary search to find interval containing xi
            left, right = 0, len(xp) - 1
            while left < right:
                mid = (left + right) // 2
                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
                if xi < xp[mid]:
                    right = mid
                else:
                    left = mid + 1
            if left == right:
                y.append(yp[left])
    return y