Improvements to zero-finding algorithm

[dbookstaber]

1. Determine if there is a zero angle before searching. Is there a heuristic that can determine whether zero is reachable? E.g., max range is typically achieved about 40 degrees.

2. MaxIterations for the proportional iterative search should probably be no more than 10. Especially if, instead of starting from elevation of zero, we take a better initial guess. For example, we could use the vacuum angle to zero as the starting point:

def calculate_drag_free_launch_angles_in_degrees(
    distance: float, height: float, velocity: float, g: float = EARTH_GRAVITY_CONSTANT
) -> Optional[Tuple[float, float]]:
    """
    Calculate both possible launch angles for hitting a target.

    Args:
        distance (float): Horizontal distance to target in meters
        height (float): Vertical height difference to target in meters
        velocity (float): Initial launch velocity in m/s

    Returns:
        Optional[Tuple[float, float]]: Tuple of (low angle, high angle) in degrees,
                                     or None if target is unreachable
    """
    # Check if target is reachable with given velocity
    velocity_squared = velocity**2
    discriminant = velocity_squared**2 - g * (
        g * distance**2 + 2 * height * velocity_squared
    )

    if discriminant < 0:
        return None

    # Calculate the two possible angles using quadratic formula
    term1 = velocity_squared
    term2 = math.sqrt(discriminant)
    term3 = g * distance

    angle1 = math.atan((term1 - term2) / term3)
    angle2 = math.atan((term1 + term2) / term3)

    # Return angles in ascending order (low, high)
    return (math.degrees(min(angle1, angle2)), math.degrees(max(angle1, angle2)))

3. When first iterative search fails, an "integrative" correction can work. From Serhiy:

    summary_error = 0
    summary_error_corrections = 0
    history = []
    while True:
        while zero_finding_error > _cZeroFindingAccuracy and iterations_count < _cMaxIterations:
            # Check height of trajectory at the zero distance (using current self.barrel_elevation)
            try:
               t = self._integrate(shot_info, zero_distance, zero_distance, TrajFlag.NONE)[0]
            except RangeError as e:
               t = e.incomplete_trajectory[0]
            height = t.height >> Distance.Foot
            height_diff = height - height_at_zero
            zero_finding_error = math.fabs(height_diff)
            summary_error+=height_diff
            if zero_finding_error > _cZeroFindingAccuracy:
                # Adjust barrel elevation to close height at zero distance
                self.barrel_elevation -= (height_diff) / zero_distance
                history.append((self.barrel_elevation, height_diff))
            else:  # last barrel_elevation hit zero!
                return Angular.Radian(self.barrel_elevation)
            iterations_count += 1
        if math.fabs(summary_error)>0 and summary_error_corrections<3:
            print(f'Got to integral correction {summary_error=}')
            print(f'Iterations: {len(history)=} minimum deviation by height {min(history, key=lambda x: abs(x[1]))=} {history=}')
            self.barrel_elevation-=(summary_error/iterations_count)/(zero_distance)
            summary_error = 0
            iterations_count=0
            summary_error_corrections += 1
        else:
            if zero_finding_error > _cZeroFindingAccuracy:
                print(f'The search has not converged  Iteration Count {len(history)} '
                      f'minimum deviation by height {min(history, key=lambda x: abs(x[1]))=} {history=}')
                # ZeroFindingError contains an instance of last barrel elevation; so caller can check how close zero is
                raise ZeroFindingError(zero_finding_error, iterations_count, Angular.Radian(self.barrel_elevation))

4. When that fails, a bracket search (e.g., https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brenth.html).

[serhiy-yevtushenko]

Here is the version, which was able to pass the test at maximal range:

The crucial details was to set

cMaxIterations: Final[int] = 60

in py_ballisticcalc/trajectory_calc/init.py

def zero_angle(self, shot_info: Shot, distance: Distance) -> Angular:
        """Iterative algorithm to find barrel elevation needed for a particular zero
        :param shot_info: Shot parameters
        :param distance: Zero distance
        :return: Barrel elevation to hit height zero at zero distance
        """
        self._init_trajectory(shot_info)

        _cZeroFindingAccuracy = self._config.cZeroFindingAccuracy
        _cMaxIterations = self._config.cMaxIterations

        distance_feet = distance >> Distance.Foot  # no need convert it twice
        zero_distance = math.cos(self.look_angle) * distance_feet
        height_at_zero = math.sin(self.look_angle) * distance_feet

        iterations_count = 0
        zero_finding_error = _cZeroFindingAccuracy * 2
        # x = horizontal distance down range, y = drop, z = windage
        while zero_finding_error > _cZeroFindingAccuracy and iterations_count < _cMaxIterations:
            # Check height of trajectory at the zero distance (using current self.barrel_elevation)
            try:
                t = self._integrate(shot_info, zero_distance, zero_distance, TrajFlag.NONE)[0]
                height = t.height >> Distance.Foot
            except RangeError as e:
                last_distance_foot = e.last_distance >> Distance.Foot
                proportion = ((last_distance_foot) / zero_distance)
                height = (e.incomplete_trajectory[-1].height>>Distance.Foot)/proportion

            zero_finding_error = math.fabs(height - height_at_zero)
 
            if zero_finding_error > _cZeroFindingAccuracy:
                # Adjust barrel elevation to close height at zero distance
                self.barrel_elevation -= (height - height_at_zero) / zero_distance
            else:  # last barrel_elevation hit zero!
                break
            iterations_count += 1
        if zero_finding_error > _cZeroFindingAccuracy:
            # ZeroFindingError contains an instance of last barrel elevation; so caller can check how close zero is
            raise ZeroFindingError(zero_finding_error, iterations_count, Angular.Radian(self.barrel_elevation))
        return Angular.Radian(self.barrel_elevation)

The drawback of solution that it takes almost a minute to find the weapon zero angle at maximum range an my laptop, but nevertheless, this version works as well at maximal range

[serhiy-yevtushenko]

I was using the following test case for checking:

import math

import pytest

from py_ballisticcalc import Distance, InterfaceConfigDict, Calculator, DragModel, TableG1, Weight, Velocity, Ammo, \
    Weapon, Shot
from py_ballisticcalc.helpers import find_index_of_apex_point

DISTANCES_FOR_CHECKING = list(range(100, 1000, 100)) + list(range(1000, 3000, 1000))+list(range(3000, 4000, 100))+\
                         list(range(4000, 7000, 500))+ list(range(6600, 7100, 100))+[7126.05]

def create_shot():
    drag_model = DragModel(
        bc=0.759,
        drag_table=TableG1,
        weight=Weight.Gram(108),
        diameter=Distance.Millimeter(23),
        length=Distance.Millimeter(108.2),
    )
    ammo = Ammo(drag_model, Velocity.MPS(930))
    gun = Weapon()
    shot = Shot(
        weapon=gun,
        ammo=ammo,
    )
    return shot

@pytest.fixture()
def zero_min_velocity_calc():
    config = InterfaceConfigDict(
        cMinimumVelocity=0,
    )
    return Calculator(_config=config)


@pytest.mark.parametrize("distance", DISTANCES_FOR_CHECKING)
def test_set_weapon_zero(distance, zero_min_velocity_calc):
    shot = create_shot()
    zero_min_velocity_calc.set_weapon_zero(shot, Distance.Meter(distance))
    print(f"{math.degrees(shot.barrel_elevation)=}")
    hit_result=zero_min_velocity_calc.fire(shot, Distance.Meter(distance), extra_data=True)
    print(f'{hit_result[-1].distance>>Distance.Meter=} {hit_result[-1].time=} {hit_result[-1].velocity>>Velocity.MPS=}')
    index_of_apex_point = find_index_of_apex_point(hit_result)
    apex_point = hit_result.trajectory[index_of_apex_point]
    print(f'{apex_point.height>>Distance.Meter=} {apex_point.distance>>Distance.Meter=} {apex_point.time=} {apex_point.velocity>>Velocity.MPS=} ')
    assert abs((hit_result[-1].distance>>Distance.Meter)-distance)<=1.0

[dbookstaber]

Added as optional test zero_finding.py in o-murphy#168

[dbookstaber]

I got significant improvements by changing the iterative correction term to incorporate the data we already compute for the trajectory angle: instead of dividing the height error by distance alone, we divide by $d (1 + \tan(\theta_0) \tan(\theta_d))$. This converges in a fraction of the iterations required with the previous formula. We can also detect that we are likely near the range limit where this doesn't work by looking at the product of tangents and switching to root-finders in those scenarios. See latest commit to the Zeroing branch, and this study that helped get there.

[dbookstaber]

Regarding the max range estimate: An upper-bound on the reachable region is given by the Envelope of Vacuum Trajectories $$Y_E = \frac{V_0^2}{2g} - \frac{g X^2}{2V_0^2}$$. (This equation from McCoy 3.28)

[serhiy-yevtushenko]

This upper bound overestimates real bound in case of drag by about 4 times

[dbookstaber]

From ZeroStudy.ipynb we can see that the v2.2 iterative zeroing method is efficient and reliable out to about 80% of max range, and that for targets at look-angles about 15-30 degrees it works well almost to the range limit. See first diagram from that notebook here:

In these images, "Iterations" is the number of trajectories that have to be calculated to reach the firing solution for the target point.

Beyond that the guaranteed-convergence methods (which start by checking the maximum range and then use a root-finder) should be preferred. The iterations to solve using these methods are shown here:

Perhaps we can nail down a heuristic that automatically selects the correct method based on a measurable characteristic of a test trajectory?

[dbookstaber]

Further detail on this: The expensive part of the current .find_zero_angle() is actually the call to .find_max_range():

The root-finding is then relatively fast: