Bicycle Simulator
Introduction
Below is the source code for a python script that models the motion of a bicycle using conservation of mechanical energy. You will need the Python interpreters to make this work, available for free at http://www.vpython.org. You can also download the Visual Python package at this website, but note this particular program does not include a graphical component.
This program can be used in an introductory physics course or laboratory to introduce physics majors to computer modeling. The program is based on the bicycle simulation discussed in Computational Physics, by Giordano and Nakanishi, chapter 2.
To run the program you need to create an input file using any text editor, such as Notepad. In the input file, you use keywords to specify the initial conditions for the simulation. Legal keywords and the types of arguments that go with them are discussed below.
The program creates a subdirectory in the same directory from which it's run. The name of the subdirectory is based on the date and current time when the simulation is executed, so that all subdirectories will have unique names. The subdirectory will contain a copy of the input file, a log file summarizing the results, and a CSV (comma separated variables) file containing the numerical output of the simulation. The CSV file can be opened in MS-Excel, and from there you can create graphical representations of the data.
Keywords
The input file uses keywords that specify the initial conditions of the simulation. Following is a list of legal keywords and the argument range for each.
- MAX_TIME, the maximum time of the simulation, in seconds (default value is 200 seconds)
- TIME_STEP, the time-step, or dt, in seconds (default is 0.1 seconds)
- GRAVITY, the acceleration due to gravity, in SI units (default value is 9.8 m/s2)
- EFFICIENCY, the mechanical efficiency of the bicycle (default value is 0.95)
- DRAG, the drag coefficient associated with air resistance (default value is 1.0)
- ROLLING, the drag coefficient associated with rolling resistance (default value is 0.003)
- AREA, effective frontal surface area of the bicycle plus rider, in square meters (default value is 0.33 m2)
- MASS, total mass of bicycle plus rider, in kilograms (default value is 70 kg)
- AIR_DENSITY, density of air, in SI units (default value is 1.247 kg/m3)
- POWER, the sustained input power provided by the rider, in watts (default is 150 W)
- WIND_SPEED, the wind speed, in SI units (default is 0)
- WIND_DIR, wind direction, measured relative to the foward direction of the bicycle, in degrees (default is 0)
- GRADE, the angle of the road, in radians (default is 0; use negative numbers for downhill angles)
Source Code
The source code appears between the bars below. Copy and paste the code into the Python IDLE environment, and hit F5 to run the code.
When running the code, you will be prompted to enter the name of the input file.
from __future__ import division
from visual import *
import time
import os
import shutil
import sys
"""
Filename: bicycle_power.py
Author: Rodney Dunning
Last edited: 18 June 2007
This program calculates the velocity at each
time step for a bicycle moving under the
influence of several resistive forces.
General algorithm:
1. Get the input file from the user.
2. Read the input file (sets initial conditions).
3. Create a subdirectory for the output files.
3. Integrate the velocity.
4. Write the output files (CSV data file, log file,
and a copy of the input file) to the subdirectory.
5. Close open files.
6. End program.
"""
class BicycleSimulation:
def __init__(self):
"""
Constructor.
"""
# SI units
self.t = 0.0 # elapsed time
self.t_max = 200.0 # max time
self.dt = 0.1 # time step
self.dt_2 = self.dt / 2 # half-step
self.x = 0 # bicycle's position
self.v = 0 # bicycle's speed
self.a = 0 # bicycle's acceleration
self.g = 9.81 # acceleration due to gravity
self.e = 0.95 # bicycle's mechanical efficiency
self.Cd = 1.0 # bicycle's drag coefficient
self.Cr = 0.003 # bicycle's rolling coefficient
self.area = 0.33 # effective frontal area
self.M = 70 # total mass (bicycle + rider)
self.rho = 1.247 # air density (at 10 deg C)
self.P = 150 # rider's power input
self.v_wind = 0 # wind speed
self.wind_dir = 0 # direction of the wind (radians)
self.grade = 0 # grade of the road (radians)
self.F = self.P / 2.68224 # constant force
# See Giordano, problem 2.5, pp. 24-25.
self.current_time_string = "now"
self.data_path = ".\/"
def get_input_file(self):
"""
Prompts the user for an input file.
Receives: nothing.
Precondition: The user specifies an input file that exists in
the same directory as the program.
Postcondition: The input file is opened in read mode.
"""
finished = 0
while not finished:
user_input = raw_input("Filename: ")
try:
self.input_file = open(name = user_input,
mode = 'r')
finished = 1
except IOError:
print "" # spacer
print "File not found."
print "Please try again."
print "" # spacer
def read_input_file(self):
"""
Parses an input file provided by the user.
Receives: nothing.
Precondition: An input file exists in the same
directory as the program that calls this method.
Postcondition: The simulation is fully initialized,
and control is passed back to the calling unit.
"""
input_file = self.input_file
input_file.seek(0)
comment = "#"
i = 0 # line number
for line in input_file.xreadlines():
nextline = line.upper().strip()
i += 1
if nextline.startswith(comment):
continue
words = nextline.split()
if len(words)==0:
continue
elif words[0] == "MAX_TIME":
self.t_max = float(words[1])
elif words[0] == "TIME_STEP":
self.dt = float(words[1])
self.dt_2 = self.dt / 2
elif words[0] == "GRAVITY":
self.g = float(words[1])
elif words[0] == "EFFICIENCY":
self.e = float(words[1])
elif words[0] == "DRAG":
self.Cd = float(words[1])
elif words[0] == "ROLLING":
self.Cr = float(words[1])
elif words[0] == "AREA":
self.area = float(words[1])
elif words[0] == "MASS":
self.M = float(words[1])
elif words[0] == "AIR_DENSITY":
self.rho = float(words[1])
elif words[0] == "POWER":
self.P = float(words[1])
self.F = self.P / 2.68224 # constant force
elif words[0] == "WIND_SPEED":
self.v_wind = float(words[1])
elif words[0] == "WIND_DIR":
self.wind_dir = float(words[1]) * acos(-1.0)/180
elif words[0] == "GRADE":
self.grade = float(words[1]) # assumes radians!
else:
print "INPUT FILE ERROR."
print "Unrecognized keyword: ", words[0]
print "Exiting program."
sys.exit()
def create_data_path(self):
"""
Creates a subdirectory where the output data files
are stored.
Receives: nothing.
Precondition: The system has been fully initialized.
Postcondition: A subdirectory is created in the same
directory as the main program file.
The subdirectory is named in part after the current
(wall) time, which insures that each simulation produces
a uniquely-name subdirectory for its data files.
"""
# Create the directory
current_time = time.localtime()
current_time_string = []
current_time_string.append(str(current_time[0]))
for i in xrange(1,6):
if current_time[i] < 10:
current_time_string.append("0" + str(current_time[i]))
else:
current_time_string.append(str(current_time[i]))
self.current_time_string = ""
for i in xrange(0,3):
self.current_time_string += current_time_string[i] + "_"
for i in xrange(3,6):
self.current_time_string += current_time_string[i]
dir_name = "Bicycle." +\
str(current_time_string[0]) + "_" +\
str(current_time_string[1]) + "_" +\
str(current_time_string[2]) + "_" +\
str(current_time_string[3]) + \
str(current_time_string[4]) +\
str(current_time_string[5])
os.mkdir(dir_name)
self.data_path = ".\/" + dir_name + ".\/" #+ file_name
def create_output_files(self):
"""
Creates the output files in the subdirectory.
Receives: nothing.
Precondition: The subdirectory has been created.
Postcondition: The output files are open in the subdirectory,
and a header line has been written to the data file. Also,
a copy of the input file is placed in the subdirectory.
"""
file_name = "data." + self.current_time_string + ".csv"
self.data_file = open(name = self.data_path + file_name, mode = 'w')
self.data_file.write("t(s), v(m/s), v(mph), a (m/s/s)\n")
file_name = "log." + self.current_time_string + ".txt"
self.log_file = open(name = self.data_path + file_name, mode = 'w')
shutil.copy(self.input_file.name, self.data_path)
def integrate(self):
inverse_mass = 1/self.M
self.a = self.F / self.M
while self.F * self.v <= self.P * self.e:
self.x += self.dt * self.v + self.dt**2 * self.a
self.v += self.dt * self.a
self.t += self.dt
self.write_data()
# main loop, using Euler's method
while self.t_max - self.t >= self.dt_2:
v_air = self.v - self.v_wind * cos(self.wind_dir)
rider_contribute = self.P * self.e / self.v
air_drag = 0.5 * self.Cd * self.rho * self.area \
* v_air**2
roll_grade = self.g * (self.Cr * cos(self.grade)\
+ sin(self.grade))
self.a = inverse_mass * (rider_contribute - air_drag) - roll_grade
v_mid = self.v + self.a * self.dt_2
v_air = v_mid - self.v_wind * cos(self.wind_dir)
rider_contribute = self.P * self.e / v_mid
air_drag = 0.5 * self.Cd * self.rho * self.area \
* v_air**2
self.a = inverse_mass * (rider_contribute - air_drag) - roll_grade
self.v += self.a * self.dt
self.x += v_mid * self.dt
self.t += self.dt
self.write_data()
self.log_file.write("End of integration loop.\n")
self.log_file.write("Time step: %3.3f s \n" %(self.dt))
self.log_file.write("Power: %3.3f W \n" %(self.P))
self.log_file.write("Total mass: %3.3f kg\n" %(self.M))
self.log_file.write("Area: %3.3f sq.m\n" %(self.area))
self.log_file.write("Wind speed: %3.3f m/s "
"(%3.3f mph)\n" %(self.v_wind,
self.v_wind*2.23693629))
self.log_file.write("Wind direction: %3.3f (rad)\n" %(self.wind_dir))
self.log_file.write("Grade: %3.3f (rad)\n" %(self.grade))
self.log_file.write("Speed at t_max = %3.2f s is %3.3f m/s "\
"(%3.3f mph)\n"%(self.t_max,
self.v,
self.v*2.23693629))
self.log_file.write("Position at t_max = %3.2f s is %3.3f m "\
"(%3.3f mi)\n"%(self.t_max,
self.x,
self.x/1609.3))
def write_data(self):
self.data_file.write("%3.10f,%3.10f,%3.10f,%3.10f\n"\
%(self.t,
self.v,
self.v*2.23693629,
self.a))
def close_files(self):
"""
Closes all open files.
"""
self.input_file.close()
self.log_file.close()
self.data_file.close()
my_bike = BicycleSimulation()
my_bike.get_input_file()
my_bike.read_input_file()
my_bike.create_data_path()
my_bike.create_output_files()
my_bike.integrate()
my_bike.close_files()
print "Program finished."