Skip to content

Commit

Permalink
Windage/spin-drift cleanup.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dbookstaber
dbookstaber committed Mar 13, 2024
1 parent 15cde78 commit 1dfa51a
Show file tree
Hide file tree
Showing 3 changed files with 133 additions and 124 deletions.
140 changes: 70 additions & 70 deletions Example.ipynb
View file

Large diffs are not rendered by default.

60 changes: 33 additions & 27 deletions py_ballisticcalc/trajectory_calc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ def _init_trajectory(self, shot_info: Shot):
self.muzzle_velocity = shot_info.ammo.get_velocity_for_temp(shot_info.atmo.temperature) >> Velocity.FPS
else:
self.muzzle_velocity = shot_info.ammo.mv >> Velocity.FPS
self.stability_coefficient = self.calc_stability_coefficient(shot_info.atmo)

def _zero_angle(self, shot_info: Shot, distance: Distance):
self._init_trajectory(shot_info)
Expand Down Expand Up @@ -183,15 +184,9 @@ def _trajectory(self, shot_info: Shot, maximum_range: float, step: float,
"""
ranges = []
ranges_length = int(maximum_range / step) + 1

time = 0
previous_mach = .0
drag = 0
stability_coefficient = 1.0
twist_coefficient = 0
if self.twist and self.length and self.diameter:
stability_coefficient = self.stability_coefficient(shot_info.atmo)
twist_coefficient = 1 if self.twist > 0 else -1

len_winds = len(shot_info.winds)
current_wind = 0
Expand Down Expand Up @@ -262,13 +257,9 @@ def _trajectory(self, shot_info: Shot, maximum_range: float, step: float,

# Record TrajectoryData row
if _flag & filter_flags:
windage = range_vector.z
if self.twist != 0: # Add spin drift to windage
windage += (1.25 * (stability_coefficient + 1.2)
* math.pow(time, 1.83) * twist_coefficient) / 12
ranges.append(create_trajectory_row(
time, range_vector, velocity_vector,
velocity, mach, windage, self.look_angle,
velocity, mach, self.spin_drift(time), self.look_angle,
density_factor, drag, self.weight, _flag.value
))
if current_item == ranges_length:
Expand Down Expand Up @@ -297,7 +288,7 @@ def _trajectory(self, shot_info: Shot, maximum_range: float, step: float,
if not filter_flags:
ranges.append(create_trajectory_row(
time, range_vector, velocity_vector,
velocity, mach, 0, self.look_angle,
velocity, mach, self.spin_drift(time), self.look_angle,
density_factor, drag, self.weight, _flag.value))
return ranges

Expand All @@ -313,20 +304,34 @@ def drag_by_mach(self, mach: float) -> float:
cd = calculate_by_curve(self._table_data, self._curve, mach)
return cd * 2.08551e-04 / self._bc

def stability_coefficient(self, atmo: Atmo) -> float:
twist_rate = math.fabs(self.twist) / self.diameter
length = self.length / self.diameter
# Miller stability formula
sd = 30 * self.weight / (
math.pow(twist_rate, 2) * math.pow(self.diameter, 3) * length * (1 + math.pow(length, 2))
)
# Velocity correction factor
fv = math.pow(self.muzzle_velocity / 2800, 1.0 / 3.0)
# Atmospheric correction
ft = atmo.temperature >> Temperature.Fahrenheit
pt = atmo.pressure >> Pressure.InHg
ftp = ((ft + 460) / (59 + 460)) * (29.92 / pt)
return sd * fv * ftp
def spin_drift(self, time) -> float:
"""Litz spin-drift approximation
:param time: Time of flight
:return: windage due to spin drift, in feet
"""
if self.twist != 0:
sign = 1 if self.twist > 0 else -1
return sign * (1.25 * (self.stability_coefficient + 1.2)
* math.pow(time, 1.83) ) / 12
return 0

def calc_stability_coefficient(self, atmo: Atmo) -> float:
"""Miller stability coefficient"""
if self.twist and self.length and self.diameter:
twist_rate = math.fabs(self.twist) / self.diameter
length = self.length / self.diameter
# Miller stability formula
sd = 30 * self.weight / (
math.pow(twist_rate, 2) * math.pow(self.diameter, 3) * length * (1 + math.pow(length, 2))
)
# Velocity correction factor
fv = math.pow(self.muzzle_velocity / 2800, 1.0 / 3.0)
# Atmospheric correction
ft = atmo.temperature >> Temperature.Fahrenheit
pt = atmo.pressure >> Pressure.InHg
ftp = ((ft + 460) / (59 + 460)) * (29.92 / pt)
return sd * fv * ftp
return 0


def wind_to_vector(wind: Wind) -> Vector:
Expand All @@ -346,8 +351,9 @@ def wind_to_vector(wind: Wind) -> Vector:


def create_trajectory_row(time: float, range_vector: Vector, velocity_vector: Vector,
velocity: float, mach: float, windage: float, look_angle: float,
velocity: float, mach: float, spin_drift: float, look_angle: float,
density_factor: float, drag: float, weight: float, flag: int) -> TrajectoryData:
windage = range_vector.z + spin_drift
drop_adjustment = get_correction(range_vector.x, range_vector.y)
windage_adjustment = get_correction(range_vector.x, windage)
trajectory_angle = math.atan(velocity_vector.y / velocity_vector.x)
Expand Down
57 changes: 30 additions & 27 deletions py_ballisticcalc_exts/py_ballisticcalc_exts/trajectory_calc.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,7 @@ cdef class TrajectoryCalc:
self.muzzle_velocity = shot_info.ammo.get_velocity_for_temp(shot_info.atmo.temperature) >> Velocity.FPS
else:
self.muzzle_velocity = shot_info.ammo.mv >> Velocity.FPS
self.stability_coefficient = self.calc_stability_coefficient(shot_info.atmo)

cdef _zero_angle(TrajectoryCalc self, object shot_info, object distance):
cdef:
Expand Down Expand Up @@ -196,15 +197,13 @@ cdef class TrajectoryCalc:
cdef _trajectory(TrajectoryCalc self, object shot_info,
double maximum_range, double step, CTrajFlag filter_flags):
cdef:
double density_factor, mach, velocity, windage, delta_time
object _flag, seen_zero # CTrajFlag
double density_factor, mach, velocity, delta_time
list ranges = []
int ranges_length = int(maximum_range / step) + 1

double time = .0
double previous_mach = .0
double drag = .0
double stability_coefficient = 1.0
int twist_coefficient = 0

int len_winds = len(shot_info.winds)
int current_item, current_wind
Expand All @@ -214,12 +213,6 @@ cdef class TrajectoryCalc:
Vector velocity_vector, velocity_adjusted
Vector range_vector, delta_range_vector, wind_vector

object _flag, seen_zero # CTrajFlag

if self.twist and self.length and self.diameter:
stability_coefficient = self.stability_coefficient(shot_info.atmo)
twist_coefficient = 1 if self.twist > 0 else -1

if len_winds < 1:
wind_vector = Vector(.0, .0, .0)
else:
Expand Down Expand Up @@ -285,13 +278,9 @@ cdef class TrajectoryCalc:

# Record TrajectoryData row
if _flag & filter_flags:
windage = range_vector.z
if self.twist != 0: # Add spin drift to windage
windage += (1.25 * (stability_coefficient + 1.2)
* pow(time, 1.83) * twist_coefficient) / 12
ranges.append(create_trajectory_row(
time, range_vector, velocity_vector,
velocity, mach, windage, self.look_angle,
velocity, mach, self.spin_drift(time), self.look_angle,
density_factor, drag, self.weight, _flag
))
if current_item == ranges_length:
Expand Down Expand Up @@ -320,7 +309,7 @@ cdef class TrajectoryCalc:
if not filter_flags:
ranges.append(create_trajectory_row(
time, range_vector, velocity_vector,
velocity, mach, 0, self.look_angle,
velocity, mach, self.spin_drift(time), self.look_angle,
density_factor, drag, self.weight, _flag))
return ranges

Expand All @@ -335,16 +324,29 @@ cdef class TrajectoryCalc:
cdef double cd = calculate_by_curve(self._table_data, self._curve, mach)
return cd * 2.08551e-04 / self._bc

cdef double stability_coefficient(self, object atmo):
cdef:
double twist_rate = fabs(self.twist) / self.diameter
double length = self.length / self.diameter
double sd = 30 * self.weight / (pow(twist_rate, 2) * pow(self.diameter, 3) * length * (1 + pow(length, 2)))
double fv = pow(self.muzzle_velocity / 2800, 1.0 / 3.0)
double ft = atmo.temperature >> Temperature.Fahrenheit
double pt = atmo.pressure >> Pressure.InHg
double ftp = ((ft + 460) / (59 + 460)) * (29.92 / pt)
return sd * fv * ftp
cdef double spin_drift(self, double time):
"""Litz spin-drift approximation
:param time: Time of flight
:return: windage due to spin drift, in feet
"""
if self.twist != 0:
sign = 1 if self.twist > 0 else -1
return sign * (1.25 * (self.stability_coefficient + 1.2) * pow(time, 1.83) ) / 12
return 0

cdef double calc_stability_coefficient(self, object atmo):
"""Miller stability coefficient"""
if self.twist and self.length and self.diameter:
cdef:
double twist_rate = fabs(self.twist) / self.diameter
double length = self.length / self.diameter
double sd = 30 * self.weight / (pow(twist_rate, 2) * pow(self.diameter, 3) * length * (1 + pow(length, 2)))
double fv = pow(self.muzzle_velocity / 2800, 1.0 / 3.0)
double ft = atmo.temperature >> Temperature.Fahrenheit
double pt = atmo.pressure >> Pressure.InHg
double ftp = ((ft + 460) / (59 + 460)) * (29.92 / pt)
return sd * fv * ftp
return 0

cdef Vector wind_to_vector(object wind):
cdef:
Expand All @@ -353,9 +355,10 @@ cdef Vector wind_to_vector(object wind):
return Vector(range_component, 0., cross_component)

cdef create_trajectory_row(double time, Vector range_vector, Vector velocity_vector,
double velocity, double mach, double windage, double look_angle,
double velocity, double mach, double spin_drift, double look_angle,
double density_factor, double drag, double weight, object flag):
cdef:
double windage = range_vector.z + spin_drift
double drop_adjustment = get_correction(range_vector.x, range_vector.y)
double windage_adjustment = get_correction(range_vector.x, windage)
double trajectory_angle = atan(velocity_vector.y / velocity_vector.x)
Expand Down

0 comments on commit 1dfa51a

Please sign in to comment.