calculation of a trajectory for example
#
# Trajectory calculation of a bullet
# You can run this script with GNU Octave ( https://www.gnu.org/software/octave/ )
# Simple create a File "trajectory.m" and copy this code in it.
# Open a Terminal and run this command:
#
# octave trajectory.m
#
#
# diameter and weight of the bullet
d = 0.008 # diameter in m
m = 0.007 # weight in kg
BC = 0.4 # ballistic coefficient
d_inch = d * 39.3701 # diameter in inch
m_pound = m * 2.20462 # mass in pound
# If the bullet is equal to the G1 standard bullet, the BC is 1.0 and the formfactor is 1.0
# The standard bullet weight 1 pound and with a diameter of 1.0 inch
# http://en.wikipedia.org/wiki/Ballistic_coefficient
formfactor = m_pound / (d_inch ^ 2.0 * BC)
# The calculation is split in time steps.
# the smaller the time steps the more accuracy is the result
timestep = 0.01 # time step in sec
# release direction
vx = 200.0 # the horizontal component of the velocity in m/s
vy = 10.0 # the vertical component of the velocity in m/s,
# this is zero for a horizontal shoot. Downwards it is lower than zero (-) and upwards it is bigger than zero (+)
# Start values for the trajectory
sx = 0.0 # horizontal component of the trajectory in m
sy = 0.0 # vertical component of the trajectory in m
rho = 1.2041 # density of the air in kg/m³
g = 9.81 # acceleration of the earth in m/s²
t = 0.0 # flight time begin with 0
A = d^2.0 * pi / 4.0 # cross section of the bullet in m²
mach_1 = 340.3; # velocity of mach = 1 in in m/s
# the mach values with the coresponding cw values of the G1 drag table
# from www.jbmballistics.com/ballistics/downloads/downloads.shtml
mach_values = [
0.00, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.50, 3.00, 4.00, 5.00
]
cw_values = [
0.26, 0.23, 0.22, 0.21, 0.20, 0.20, 0.21, 0.25, 0.34, 0.48, 0.58, 0.63, 0.65, 0.66, 0.65, 0.64, 0.63, 0.62, 0.60, 0.59, 0.53, 0.51, 0.50, 0.49
]
# run while the range is smaller than 350m.
# each pass of the loop one time step is calculated.
while sx < 350.0
printf("------------------------------------------------------------\n")
# the total amount of the velocity.
v = sqrt(vx^2.0 + vy^2.0)
# fetch the cw value from the G1 table:
mach = v / mach_1
cw = spline(mach_values,cw_values,mach)
# force of the air drag in N
# http://de.wikipedia.org/wiki/Luftwiderstand
F = 0.5 * rho * A * v^2.0 * cw * formfactor
# split the drag force in a horizontal and a vertical component.
Fx = F * vx / v
Fy = F * vy / v
# horizontal and vertical component of the retardment (negativ acceleration) in m/s²
# The vertical component is add with the earth accelaration.
ax = -Fx / m
ay = -Fy / m - g
# diff of the velocity
delta_vx = ax * timestep
delta_vy = ay * timestep
# distance in the time step
delta_sx = vx * timestep + delta_vx * timestep / 2.0
delta_sy = vy * timestep + delta_vy * timestep / 2.0
# add the distance in the time step to the whole distance
sx = sx + delta_sx
sy = sy + delta_sy
# velocities after the time step
vx = vx + delta_vx
vy = vy + delta_vy
# add the timestep to the total amount of time
t = t + timestep
endwhile