# 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. Upwards it is lower than zero (-) and downwards 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