Created
July 25, 2025 18:32
-
-
Save andriitishchenko/b7bfadee2eab355248fc28b7f418f050 to your computer and use it in GitHub Desktop.
Arrow Flight Calculations
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# https://colab.research.google.com/drive/1iUvsAFJ1hCD3qeSczQ7So6UYxWzrtiX5?usp=sharing | |
# === Arrow Ballistics Simulation === | |
# This script simulates the flight of an arrow, first without fletching (vanes), then with fletching. | |
# It optimizes the launch angle to hit a target at 70m, logs all key parameters at each step, and plots the trajectory with fletching. | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from scipy.optimize import minimize_scalar | |
# === PHYSICAL CONSTANTS === | |
GRAINS_TO_KG = 0.00006479891 # 1 grain = 0.00006479891 kg | |
ARROW_MASS = 320 * GRAINS_TO_KG # Arrow mass [kg] | |
DIAMETER = 0.0052 # Shaft diameter [m] | |
AREA = np.pi * (DIAMETER / 2) ** 2 # Shaft cross-section area [m^2] | |
# === FLETCHING PARAMETERS === | |
FLETCH_COUNT = 3 # Number of vanes (fletchings) | |
FLETCH_LENGTH = 0.025 # Vane length [m] | |
FLETCH_HEIGHT = 0.010 # Vane height [m] | |
FLETCH_THICKNESS = 0.0001 # Vane thickness [m] | |
# Parabolic area formula: (2/3) * base * height | |
FLETCH_AREA = (2/3) * FLETCH_LENGTH * FLETCH_HEIGHT # Area of one vane [m^2] | |
TOTAL_FLETCH_AREA = FLETCH_COUNT * FLETCH_AREA # Total fletching area [m^2] | |
# Angle of attack (cant of fletching relative to flow) | |
FLETCH_ALPHA_DEG = 3.0 # Fletching cant angle [deg] | |
FLETCH_ALPHA_RAD = np.radians(FLETCH_ALPHA_DEG) # [rad] | |
# Edge-on area: length * thickness | |
FLETCH_EDGE_AREA = FLETCH_LENGTH * FLETCH_THICKNESS | |
TOTAL_FLETCH_EDGE_AREA = FLETCH_COUNT * FLETCH_EDGE_AREA | |
# Base area: thickness * height (end of vane) | |
FLETCH_BASE_AREA = FLETCH_THICKNESS * FLETCH_HEIGHT | |
TOTAL_FLETCH_BASE_AREA = FLETCH_COUNT * FLETCH_BASE_AREA | |
# Drag coefficients (empirical, see literature NASA Glenn: Drag Coefficient Database): | |
FLETCH_CD = 1.5 # 0.8 # Main fletching drag (perpendicular to flow) # 1.2 – 1.8 | |
FLETCH_EDGE_CD = 0.15 # 0.2 # Edge-on drag (thin plate) # 0.05 – 0.15 | |
FLETCH_BASE_CD = 0.4 # 0.05 # Base drag (end of vane) #0.3 – 0.5 | |
# === ENVIRONMENT PARAMETERS === | |
RHO = 1.225 # Air density [kg/m^3] | |
MU = 1.81e-5 # Dynamic viscosity of air [kg/(m·s)] | |
G = 9.81 # Gravity [m/s^2] | |
# === LAUNCH CONDITIONS === | |
V0 = 50.0 # Initial speed [m/s] | |
H0 = 1.5 # Launch height [m] | |
TARGET_X = 70.0 # Target distance [m] | |
TARGET_Y = 1.3 # Target height [m] | |
TARGET_RADIUS = 0.7 # Target hit radius [m] (1.3 m center, 70 cm radius) | |
# === DRAG FORCE FUNCTIONS === | |
# Shaft only (no fletching) | |
def drag_force_shaft(v): | |
""" | |
Calculate drag force for arrow shaft (no fletching). | |
v: speed [m/s] | |
Returns: drag force [N] | |
Uses updated Cd ranges based on recent arrow aerodynamic studies: | |
- Laminar regime (Re < 1.2e4): Cd ≈ 1.5 | |
- Transition (1.2e4 ≤ Re < 2.0e4): linear transition to Cd ≈ 2.6 | |
- Turbulent regime (Re ≥ 2.0e4): Cd ≈ 2.6 (may vary with shape, angle) | |
Sources: | |
- Aerodynamic properties of an archery arrow, ResearchGate (2012) | |
https://www.researchgate.net/publication/235908721_Aerodynamic_properties_of_an_archery_arrow | |
- Archery flight dynamics blog | |
https://archerycoachesguildblog2.wordpress.com/technical_library/arrows-flight-dynamics/ | |
""" | |
Re = (RHO * v * DIAMETER) / MU | |
Cd_laminar = 1.5 | |
Cd_turbulent = 2.6 | |
Re_laminar_start = 1.2e4 | |
Re_transition_end = 2.0e4 | |
if Re < Re_laminar_start: | |
cd = Cd_laminar | |
elif Re_laminar_start <= Re < Re_transition_end: | |
# Linear interpolation between laminar and turbulent Cd | |
cd = Cd_laminar + (Re - Re_laminar_start) / (Re_transition_end - Re_laminar_start) * (Cd_turbulent - Cd_laminar) | |
else: | |
cd = Cd_turbulent | |
drag = 0.5 * cd * RHO * AREA * v**2 | |
return drag | |
# Fletching drag calculation (modular, all aerodynamic effects) | |
def fletching_drag(v, alpha_deg, params, omega=None): | |
""" | |
Calculate total fletching drag force for given velocity, fletching cant angle, and arrow spin. | |
Considers: | |
- Main drag (perpendicular area, depends on angle) | |
- Edge-on drag (thin plate) | |
- Base drag (end of vane) | |
- Gyroscopic stabilization (reduces side drag) | |
- Spin-dependent Cd correction | |
Parameters: | |
v: arrow velocity [m/s] | |
alpha_deg: fletching cant angle [deg] | |
params: dict with keys: | |
'count', 'area', 'edge_area', 'base_area', | |
'cd_main', 'cd_edge', 'cd_base', 'rho' | |
omega: arrow spin rate [rad/s] (optional) | |
Returns: | |
Total fletching drag force [N] | |
""" | |
alpha_rad = np.radians(alpha_deg) | |
# Main fletching drag (perpendicular to flow, depends on angle) | |
effective_area = params['count'] * params['area'] * np.sin(alpha_rad) | |
# --- Spin dynamics --- | |
# Estimate spin rate if not provided (simple model: proportional to alpha and velocity) | |
if omega is None: | |
# Empirical: 1 deg cant ~ 1 turn per 2m, so omega ~ v * alpha_rad / (2 * DIAMETER) | |
omega = v * alpha_rad / (2 * DIAMETER) | |
# Gyroscopic stabilization factor: higher omega reduces side drag | |
# Model: stabilization_factor = 1 / (1 + omega/omega_ref), omega_ref ~ 200 rad/s | |
omega_ref = 200.0 | |
stabilization_factor = 1.0 / (1.0 + omega / omega_ref) | |
# Spin-dependent Cd correction: reduce Cd_main by up to 20% at high spin | |
cd_main_eff = params['cd_main'] * (1.0 - 0.2 * min(omega / (2*omega_ref), 1.0)) | |
# Main fletching drag (perpendicular to flow, depends on angle and stabilization) | |
F_main = 0.5 * cd_main_eff * params['rho'] * effective_area * v**2 * stabilization_factor | |
# Edge-on drag (thin plate, less affected by spin) | |
F_edge = 0.5 * params['cd_edge'] * params['rho'] * params['count'] * params['edge_area'] * v**2 | |
# Base drag (end of vane, less affected by spin) | |
F_base = 0.5 * params['cd_base'] * params['rho'] * params['count'] * params['base_area'] * v**2 | |
# For logging/debugging, you can print omega, stabilization_factor, cd_main_eff if needed | |
return F_main + F_edge + F_base | |
# Shaft + fletching (all components) | |
def drag_force_full(v): | |
F_shaft = drag_force_shaft(v) | |
# Fletching drag (all effects, with spin) | |
fletch_params = { | |
'count': FLETCH_COUNT, | |
'area': FLETCH_AREA, | |
'edge_area': FLETCH_EDGE_AREA, | |
'base_area': FLETCH_BASE_AREA, | |
'cd_main': FLETCH_CD, | |
'cd_edge': FLETCH_EDGE_CD, | |
'cd_base': FLETCH_BASE_CD, | |
'rho': RHO | |
} | |
# Estimate spin rate for current velocity and fletching angle | |
alpha_deg = FLETCH_ALPHA_DEG | |
omega = v * np.radians(alpha_deg) / (2 * DIAMETER) | |
F_fletch = fletching_drag(v, alpha_deg, fletch_params, omega=omega) | |
return F_shaft + F_fletch | |
# === SIMULATION FUNCTION === | |
def simulate(angle_deg, drag_func, dt=0.05, max_time=10, log_speeds=False): | |
""" | |
Simulate arrow flight for given launch angle and drag function. | |
Logs all key parameters at each step if log_speeds=True. | |
Returns: x_hit, v_hit, angle_deg, xs, ys, vs | |
""" | |
angle = np.radians(angle_deg) | |
vx, vy = V0 * np.cos(angle), V0 * np.sin(angle) | |
x, y = 0.0, H0 | |
xs, ys, vs = [x], [y], [V0] | |
t = 0 | |
if log_speeds: | |
print(f"\n--- Simulation log for angle {angle_deg:.2f}° ---") | |
print("t [s] | x [m] | y [m] | v [m/s] | Re | F_drag [N] | Slow speed [%]") | |
while y >= 0 and t < max_time: | |
v = np.sqrt(vx**2 + vy**2) | |
Re = (RHO * v * DIAMETER) / MU | |
Fd = drag_func(v) | |
ax = -Fd * vx / (v * ARROW_MASS) | |
ay = -G - (Fd * vy) / (v * ARROW_MASS) | |
vx += ax * dt | |
vy += ay * dt | |
x += vx * dt | |
y += vy * dt | |
t += dt | |
slow_speed = (V0 - vs[-1]) / V0 * 100 | |
xs.append(x) | |
ys.append(y) | |
vs.append(np.sqrt(vx**2 + vy**2)) | |
if log_speeds: | |
print(f"{t:.3f} | {x:.2f} | {y:.2f} | {vs[-1]:.2f} | {Re:.0f} | {Fd:.3f} | {slow_speed:.1f}%") | |
# Основное условие: пересечение TARGET_Y | |
if ys[-2] > TARGET_Y > y : | |
ratio = (TARGET_Y - y) / (ys[-2] - y) | |
x_cross = x + ratio * (xs[-2] - x) | |
v_hit = vs[-1] + ratio * (vs[-2] - vs[-1]) | |
return x_cross, v_hit, angle_deg, xs, ys, vs | |
return None, None, angle_deg, xs, ys, vs | |
# === MAIN SIMULATION LOGIC === | |
# 1. Optimize angle without fletching | |
def objective_shaft(angle_deg): | |
result = simulate(angle_deg, drag_force_shaft) | |
x_hit, _, _, xs, ys, _ = result | |
# Проверка попадания: если simulate вернул x_hit, значит попадание | |
if x_hit is not None: | |
# Попадание только если x_hit == TARGET_X (с точностью до шага) | |
return abs(x_hit - TARGET_X) | |
return 1e6 | |
result_shaft = minimize_scalar(objective_shaft, bounds=(0, 20), method='bounded') | |
opt_angle_shaft = result_shaft.x | |
x_hit_shaft, v_final_shaft, _, _, _, vs_shaft = simulate(opt_angle_shaft, drag_force_shaft, log_speeds=True) | |
print(f"\n[NO FLETCHING] Launch angle: {opt_angle_shaft:.3f}°") | |
if v_final_shaft: | |
print(f"Final speed at {TARGET_X}m: {v_final_shaft:.2f} m/s") | |
else: | |
print("Arrow did not reach target") | |
# 2. Simulate with fletching, using found angle | |
x_hit_fletch, v_final_fletch, _, traj_x_fletch, traj_y_fletch, vs_fletch = simulate(opt_angle_shaft, drag_force_full, log_speeds=True) | |
print(f"\n[WITH FLETCHING] Launch angle: {opt_angle_shaft:.3f}°") | |
if v_final_fletch: | |
print(f"Final speed at {TARGET_X}m: {v_final_fletch:.2f} m/s") | |
else: | |
print("Arrow did not reach target") | |
# Simulate bare shaft trajectory for plotting | |
_, _, _, traj_x_bare, traj_y_bare, _ = simulate(opt_angle_shaft, drag_force_shaft) | |
# Plot both trajectories | |
plt.figure(figsize=(10, 4)) | |
plt.plot(traj_x_bare, traj_y_bare, label='Bare shaft', color='blue') | |
plt.plot(traj_x_fletch, traj_y_fletch, label='Fletched arrow', color='green') | |
plt.axhline(y=TARGET_Y, color='r', linestyle='--', label='Target center (1.3 m)') | |
plt.axvline(x=TARGET_X, color='k', linestyle='--', label=f"{TARGET_X} m") | |
plt.xlabel('Distance (m)') | |
plt.ylabel('Height (m)') | |
plt.title(f"Trajectories - Angle: {opt_angle_shaft:.2f}°") | |
plt.legend() | |
plt.grid(True) | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment