How to plot the acceleration of the boulder and or the grains?

Dear all,

I am trying to plot the velocity and acceleration:

1. How to add plot data for acceleration (for the boulder and or the grains)?

2. Based on my script, is that the way correct if I wanna plot velocity for the boulder like in my script?

##SCRIPT

#Sphere Cylinder pack
readParamsFromTable(maxLoad=1000,minLoad=100)
from yade.params.table import *
from yade import pack, plot
from pylab import rand

O.bodies.append(geom.facetCylinder(Vector3(1,1,25),4,50,wallMask=6))

#Define material of the grains
O.materials.append(FrictMat(young=1e7,poisson=.3,density=2000,frictionAngle=20))

sp=pack.SpherePack()
sp.makeCloud((-1.6,2.6,.1),(3,-1,20),rMean=.3,rRelFuzz=.5,num=950,)
sp.toSimulation(color=(0.6+0.15*rand(),0.5+0.15*rand(),0.15+0.15*rand()))

O.bodies.append(wall((0,0,0),axis=2))

##Engines and Constitutive Law
O.engines=[ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]), #move for sphere, facet, wall
  InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()], #interaction between them
   [Ip2_FrictMat_FrictMat_FrictPhys()], # ip2 list (just one list!)
   [Law2_ScGeom_FrictPhys_CundallStrack()]), # law2 list
   NewtonIntegrator(damping=.05,gravity=[0,0,-9.81]),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),]

O.dt=.5*PWaveTimeStep()

def checkUnbalanced():
 if O.iter<1000:return # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if unbalancedForce()>1:return # the rest will be run only if unbalanced is < 1 (stabilized packing)
 O.bodies.append([utils.sphere(center=(1,1,50),radius=3.7)])
 O.materials.append(FrictMat(young=100e7,poisson=.2,density=2700))
 global boulder # without this line, the BOULDER variable would only exist inside this function
 boulder=O.bodies[-1] # the last particles is the boulder
 boulder.state.vel=(0,0,-.1) # start plotting the data now, it was not interesting before
 O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=500)]
 checker.command='unloadBoulder()'

def unloadBoulder():
 if abs(O.forces.f(boulder.id)[2])>maxLoad:
  boulder.state.vel*=-1 # prescribing a velocity will therefore make it move at constant velocity (downwards)
  checker.command='stopUnloading()'

def stopUnloading():
 if abs(O.forces.f(boulder.id)[2])<minLoad:
  O.run()

def addPlotData():
  b=O.bodies[-1]
  plot.addData(position=b.state.pos[2], velocity=b.state.vel.norm(), itera=O.iter, time=O.time)

plot.plots={'time':('velocity',), 'itera':('position',)}
plot.plot()

yade.qt.Controller();

O.saveTmp()

############################################

Thank you for your any suggestions.

Cheers!