Skip to main content

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