The code above is naturally implemented as a Python function. This function can take the most important physical parameters of the problem as input, along with information about the file with road shapes. We allow for defining road shapes either through a file on a web site or a local file.
def solve(url=None, m=60, b=80, k=60, v=5):
"""
Solve model for verticle vehicle vibrations.
========= ==============================================
variable description
========= ==============================================
url either URL of file with excitation force data,
or name of a local file
m mass of system
b friction parameter
k spring parameter
v (constant) velocity of vehicle
Return data (list) holding input and output data
[x, t, [h,a,u], [h,a,u], ...]
========= ==============================================
"""
# Download file (if url is not the name of a local file)
if url.startswith('http://') or url.startswith('file://'):
import urllib
filename = os.path.basename(url) # strip off path
urllib.urlretrieve(url, filename)
else:
# Check if url is the name of a local file
filename = url
if not os.path.isfile(filename):
print url, 'must be a URL or a filename'
sys.exit(1) # abort program
# else: ok
h_data = np.loadtxt(filename) # read numpy array from file
x = h_data[0,:] # 1st column: x coordinates
h_data = h_data[1:,:] # other columns: h shapes
t = x/v # time corresponding to x
dt = t[1] - t[0]
def f(u):
return k*u
data = [x, t] # key input and output data (arrays)
for i in range(h_data.shape[0]):
h = h_data[i,:] # extract a column
a = acceleration(h, x, v)
u = forced_vibrations(t=t, I=0.2, m=m, b=b, f=f, F=-m*a)
data.append([h, a, u])
return data
Note that function arguments can be given default values (known as keyword arguments in Python). Python has a lot of operating system functionality, such as checking if a file, directory, or link exist, creating or removing files and directories, running stand-alone applications, etc.
Since the roads have a quite noise shape, the force \( F=-ma \) looks
very noisy, while the response \( u(t) \) to this excitation is
significantly less noisy,
see Figure 5 for an example.
It may be useful to compute the root mean square value of the various
realizations of \( u(t) \) and add this array to the data list of input
and output data in the problem:
def rms(data):
t = data[1]
u_rms = np.zeros(t.size) # for accumulating the rms value
for h, a, u in data[2:]: # loop over results
u_rms += u**2
u_rms = np.sqrt(u_rms/u_rms.size)
data.append(u_rms)
return data
Observe here the for loop where three variables (h, a, and u)
are set equal to the three arrays in each element of the sublist
data[2:]. We only need the u array for the computation. An alternative
loop would be to have one loop variable for each list element and
get u as the 3rd element of each list element:
for res in data[2:]:
u = res[2] # res is [h,a,u]
...
After calling
road_url = 'http://hplbit.bitbucket.org/data/bumpy/bumpy.dat.gz'
data = solve(url=road_url,
m=60, b=200, k=60, v=6)
data = rms(data)
the data array contains single arrays and triplets of arrays,
[x, t, [h,a,u], [h,a,u], ..., [h,a,u], u_rms]
This list, or any Python object, can be stored on file for later retrieval of the results, using the pickling functionality in Python:
import cPickle
outfile = open('bumpy.res', 'w')
cPickle.dump(data, outfile)
outfile.close()