closest Wind using bisect

README.md

@@ -1,6 +1,37 @@
 # BallisticCalculator
 LGPL library for small arms ballistic calculations based on point-mass (3 DoF) plus spin drift.
 
+[![license]][LGPL-3]
+[![pypi]][PyPiUrl]
+[![downloads]][pepy]
+[![downloads/month]][pepy]
+[![versions]][sources]
+[![Made in Ukraine]][SWUBadge]
+
+
+[sources]:
+https://github.com/o-murphy/py-ballisticcalc
+[license]:
+https://img.shields.io/github/license/o-murphy/py-ballisticcalc?style=flat-square
+[LGPL-3]:
+https://opensource.org/licenses/LGPL-3.0-only
+[pypi]:
+https://img.shields.io/pypi/v/py-ballisticcalc?style=flat-square&logo=pypi
+[PyPiUrl]:
+https://pypi.org/project/py-ballisticcalc/
+[downloads]:
+https://img.shields.io/pepy/dt/py-ballisticcalc?style=flat-square
+[downloads/month]:
+https://static.pepy.tech/personalized-badge/py-ballisticcalc?style=flat-square&period=month&units=abbreviation&left_color=grey&right_color=blue&left_text=downloads%2Fmonth
+[pepy]:
+https://pepy.tech/project/py-ballisticcalc
+[versions]:
+https://img.shields.io/pypi/pyversions/py-ballisticcalc?style=flat-square
+[Made in Ukraine]:
+https://img.shields.io/badge/made_in-Ukraine-ffd700.svg?labelColor=0057b7&style=flat-square
+[SWUBadge]:
+https://stand-with-ukraine.pp.ua
+
 ### Table of contents
 * **[Installation](#installation)**
   * [Latest stable](#latest-stable-release-from-pypi)

py_ballisticcalc.exts/py_ballisticcalc_exts/trajectory_calc.pyx

@@ -138,10 +138,13 @@ cdef class _TrajectoryDataFilter:
 cdef class _WindSock:
     cdef:
         list[Wind_t] winds
-        int current
-        double next_range
-        CVector _last_vector_cache
-        int _length
+        list[double] distances
+        int previous_wind_index
+        double previous_range
+        # int current
+        # double next_range
+        # CVector _last_vector_cache
+        # int _length
 
     def __cinit__(_WindSock self, object winds):
         self.winds = [
@@ -152,37 +155,98 @@ cdef class _WindSock:
                 w.MAX_DISTANCE_FEET
             ) for w in winds
         ]
-        self.current = 0
-        self.next_range = cMaxWindDistanceFeet
-        self._last_vector_cache = CVector(0.0, 0.0, 0.0)
-        self._length = len(self.winds)
-
-        # Initialize cache correctly
-        self.update_cache()
-
-    cdef CVector current_vector(_WindSock self):
-        return self._last_vector_cache
-
-    cdef void update_cache(_WindSock self):
-        cdef Wind_t cur_wind
-        if self.current < self._length:
-            cur_wind = self.winds[self.current]
-            self._last_vector_cache = wind_to_c_vector(&cur_wind)
-            self.next_range = cur_wind.until_distance
+        self.distances = [wind.until_distance._feet for wind in winds]
+        self.previous_wind_index = -1
+        self.previous_range = -1.0
+
+    cdef bint should_change_wind(self, double next_range):
+        cdef int new_index
+        # Skip processing if the range hasn't changed since the last call
+        if self.previous_range == next_range:
+            return False
+
+        # Find the closest vector for the new range
+        new_index = self.find_closest_index(next_range)
+
+        # Check if the wind region has changed
+        if new_index == self.previous_wind_index:
+            # Update the cached range but not the vector
+            self.previous_range = next_range
+            return False
+
+        # Update cache with the new index and range
+        self.previous_wind_index = new_index
+        self.previous_range = next_range
+        return True
+
+    cdef CVector vector_for_range(self, double next_range):
+        """
+        Returns the closest or interpolated wind vector for the given range if the wind region changes.
+        Otherwise, returns None if the wind vector does not need updating.
+        """
+        return self.find_closest_vector(next_range)
+
+    cdef int find_closest_index(self, double target_range):
+        cdef int start_index
+        # Start searching from the cached index to optimize
+        if self.previous_wind_index != -1:
+            start_index = max(0, min(self.previous_wind_index, len(self.distances) - 1))
         else:
-            self._last_vector_cache = CVector(0.0, 0.0, 0.0)
-            self.next_range = cMaxWindDistanceFeet
+            start_index = 0
 
-    cdef CVector vector_for_range(_WindSock self, double next_range):
-        if next_range >= self.next_range:
-        # if next_range + 1e-6 >= self.next_range:
-            self.current += 1
-            if self.current >= self._length:
-                self._last_vector_cache = CVector(0.0, 0.0, 0.0)
-                self.next_range = cMaxWindDistanceFeet
-            else:
-                self.update_cache()  # This will trigger cache updates.
-        return self._last_vector_cache
+        # Find the appropriate position using bisect_left starting from the cached index
+        return bisect_left(self.distances, target_range, lo=start_index)
+
+
+    cdef CVector find_closest_vector(self, double target_range):
+        """
+        Finds and returns an interpolated wind vector for the given range (target_range) and the wind index:
+        - If the range is within the distance of two winds, interpolate between their vectors.
+        - If the range surpasses the maximum wind distance, return a zero vector.
+        """
+
+        cdef:
+            int pos
+            Wind_t lower_wind, upper_wind, wind_0
+            double lower_distance, upper_distance, factor
+            CVector lower_wind_vector, upper_wind_vector, interpolated_vector
+            CVector adjusted_lower_wind_vector, adjusted_upper_wind_vector
+
+        if not self.winds:
+            return CVector(0.0, 0.0, 0.0)  # If there are no winds, return a zero vector and invalid index.
+
+        # Check if the target range exceeds the maximum wind distance
+        if target_range > self.distances[-1]:
+            return CVector(0.0, 0.0, 0.0)
+
+        pos = self.find_closest_index(target_range)
+
+        if pos == 0:
+            # Target is smaller than the smallest distance, return the first wind vector
+            wind_0 = self.winds[0]
+            return wind_to_c_vector(&wind_0)
+        else:
+            # Get the two closest winds for interpolation
+            lower_wind = self.winds[pos - 1]
+            upper_wind = self.winds[pos]
+
+            lower_distance = self.distances[pos - 1]
+            upper_distance = self.distances[pos]
+
+            if target_range < lower_distance:
+                # Return the lower wind if the target is exactly at or before it
+                return wind_to_c_vector(&lower_wind)
+
+            # Calculate the interpolation factor (0 <= factor <= 1)
+            factor = (target_range - lower_distance) / (upper_distance - lower_distance)
+
+            # Interpolate between the two vectors
+            lower_wind_vector = wind_to_c_vector(&lower_wind)
+            upper_wind_vector = wind_to_c_vector(&upper_wind)
+            adjusted_lower_wind_vector = mul_c(&lower_wind_vector, 1 - factor)
+            adjusted_upper_wind_vector = mul_c(&upper_wind_vector, factor)
+            interpolated_vector =  add(&adjusted_lower_wind_vector, &interpolated_vector)
+            return interpolated_vector
 
 
 cdef class TrajectoryCalc:
@@ -328,7 +392,8 @@ cdef class TrajectoryCalc:
             double calc_step = self.__shot.calc_step
 
             # region Initialize wind-related variables to first wind reading (if any)
-            CVector wind_vector = self.ws.current_vector()
+            # CVector wind_vector = self.ws.current_vector()
+            CVector wind_vector = self.ws.vector_for_range(0)
             # endregion
 
             _TrajectoryDataFilter data_filter
@@ -367,7 +432,9 @@ cdef class TrajectoryCalc:
             data_filter.current_flag = CTrajFlag.NONE
 
             # Update wind reading at current point in trajectory
-            if range_vector.x >= self.ws.next_range:  # require check before call to improve performance
+            # if range_vector.x >= self.ws.next_range:  # require check before call to improve performance
+            #     wind_vector = self.ws.vector_for_range(range_vector.x)
+            if self.ws.should_change_wind(range_vector.x):
                 wind_vector = self.ws.vector_for_range(range_vector.x)
 
             # overwrite density_factor and mach by pointer

py_ballisticcalc/trajectory_calc/_trajectory_calc.py

@@ -4,6 +4,7 @@
 
 import math
 import warnings
+from bisect import bisect_left
 
 from typing_extensions import NamedTuple, Union, List, Tuple
 
@@ -140,52 +141,94 @@ def check_zero_crossing(self, range_vector: Vector):
                     self.current_flag |= TrajFlag.ZERO_DOWN
                     self.seen_zero |= TrajFlag.ZERO_DOWN
 
-
 class _WindSock:
     """
-    Currently this class assumes that requests for wind readings will only be made in order of increasing range.
-    This assumption is violated if the projectile is blown or otherwise moves backwards.
+    Finds and returns the wind vector for the given range (target_range) based on thresholds.
+    Caches the previous wind index to optimize repeated calls for similar ranges.
     """
     winds: tuple['Wind', ...]
-    current: int
-    next_range: float
+    distances: list[float]
+    previous_wind_index: int  # Cache to store the last used wind index
+    previous_range: float
 
     def __init__(self, winds: Union[Tuple["Wind", ...], None]):
         self.winds: Tuple["Wind", ...] = winds or tuple()
-        self.current: int = 0
-        self.next_range: float = Wind.MAX_DISTANCE_FEET
-        self._last_vector_cache: Union["Vector", None] = None
-        self._length = len(self.winds)
-
-        # Initialize cache correctly
-        self.update_cache()
-
-    def current_vector(self) -> "Vector":
-        """Returns the current cached wind vector."""
-        if not self._last_vector_cache:
-            raise RuntimeError("No cached wind vector")
-        return self._last_vector_cache
-
-    def update_cache(self) -> None:
-        """Updates the cache only if needed or if forced during initialization."""
-        if self.current < self._length:
-            cur_wind = self.winds[self.current]
-            self._last_vector_cache = cur_wind.vector
-            self.next_range = cur_wind.until_distance >> Distance.Foot
+        self.distances = [wind.until_distance >> Distance.Foot for wind in self.winds]
+        self.previous_wind_index: int = -1  # Initialize with an invalid index
+        self.previous_range: float = -1.0  # Initialize with an invalid range
+
+    def should_change_wind(self, next_range) -> bool:
+        # Skip processing if the range hasn't changed since the last call
+        if self.previous_range == next_range:
+            return False
+
+        # Find the closest vector for the new range
+        new_index = self.find_closest_index(next_range)
+
+        # Check if the wind region has changed
+        if new_index == self.previous_wind_index:
+            # Update the cached range but not the vector
+            self.previous_range = next_range
+            return False
+
+        # Update cache with the new index and range
+        self.previous_wind_index = new_index
+        self.previous_range = next_range
+        return True
+
+    def vector_for_range(self, next_range: float) -> Union["Vector", None]:
+        """
+        Returns the closest or interpolated wind vector for the given range if the wind region changes.
+        Otherwise, returns None if the wind vector does not need updating.
+        """
+        return self.find_closest_vector(next_range)
+
+    def find_closest_index(self, target_range: float):
+        # Start searching from the cached index to optimize
+        if self.previous_wind_index != -1:
+            start_index = max(0, min(self.previous_wind_index, len(self.distances) - 1))
         else:
-            self._last_vector_cache = Vector(0.0, 0.0, 0.0)
-            self.next_range = Wind.MAX_DISTANCE_FEET
+            start_index = 0
+
+        # Find the appropriate position using bisect_left starting from the cached index
+        return bisect_left(self.distances, target_range, lo=start_index)
+
+    def find_closest_vector(self, target_range: float) -> "Vector":
+        """
+        Finds and returns an interpolated wind vector for the given range (target_range) and the wind index:
+        - If the range is within the distance of two winds, interpolate between their vectors.
+        - If the range surpasses the maximum wind distance, return a zero vector.
+        """
+        if not self.winds:
+            return Vector(0.0, 0.0, 0.0)  # If there are no winds, return a zero vector and invalid index.
+
+        # Check if the target range exceeds the maximum wind distance
+        if target_range > self.distances[-1]:
+            return Vector(0.0, 0.0, 0.0)
+
+        pos = self.find_closest_index(target_range)
+
+        if pos == 0:
+            # Target is smaller than the smallest distance, return the first wind vector
+            return self.winds[0].vector
+        else:
+            # Get the two closest winds for interpolation
+            lower_wind = self.winds[pos - 1]
+            upper_wind = self.winds[pos]
+
+            lower_distance = self.distances[pos - 1]
+            upper_distance = self.distances[pos]
+
+            if target_range < lower_distance:
+                # Return the lower wind if the target is exactly at or before it
+                return lower_wind.vector
+
+            # Calculate the interpolation factor (0 <= factor <= 1)
+            factor = (target_range - lower_distance) / (upper_distance - lower_distance)
 
-    def vector_for_range(self, next_range: float) -> "Vector":
-        """Updates the wind vector if `next_range` surpasses `self.next_range`."""
-        if next_range >= self.next_range:
-            self.current += 1
-            if self.current >= self._length:
-                self._last_vector_cache = Vector(0.0, 0.0, 0.0)
-                self.next_range = Wind.MAX_DISTANCE_FEET
-            else:
-                self.update_cache()  # This will trigger cache updates.
-        return self.current_vector()
+            # Interpolate between the two vectors
+            interpolated_vector = lower_wind.vector * (1 - factor) + upper_wind.vector * factor
+            return interpolated_vector
 
 
 # pylint: disable=too-many-instance-attributes
@@ -311,7 +354,7 @@ def _integrate(self, shot_info: Shot, maximum_range: float, step: float,
 
         # region Initialize wind-related variables to first wind reading (if any)
         wind_sock = _WindSock(shot_info.winds)
-        wind_vector = wind_sock.current_vector()
+        wind_vector = wind_sock.vector_for_range(0)
         # endregion
 
         # region Initialize velocity and position of projectile
@@ -339,7 +382,7 @@ def _integrate(self, shot_info: Shot, maximum_range: float, step: float,
             data_filter.clear_current_flag()
 
             # Update wind reading at current point in trajectory
-            if range_vector.x >= wind_sock.next_range:  # require check before call to improve performance
+            if wind_sock.should_change_wind(range_vector.x):
                 wind_vector = wind_sock.vector_for_range(range_vector.x)
 
             # Update air density at current point in trajectory

tests/test_computer.py

@@ -2,10 +2,13 @@
 
 import unittest
 import copy
+import logging
 from py_ballisticcalc import (
     DragModel, Ammo, Weapon, Calculator, Shot, Wind, Atmo, TableG7, RangeError,
 )
 from py_ballisticcalc.unit import *
+from py_ballisticcalc.logger import logger
+logger.setLevel(logging.DEBUG)
 
 
 class TestComputer(unittest.TestCase):