mistakes in calculating shearInc

Hi all,

I have problems in well understanding the code inside 'ScGeom.cpp' when calculating shearInc,
the results given by my python code is quite different from what given by yade simulations.
Does anyone know what's wrong in my python code? Here it is:

import numpy as np
i=O.interactions[237,259] ## choose a random contact for instance
n=np.hstack(i.geom.normal) ## transform a vector into an array to simply the programming process
p=i.geom.penetrationDepth
ra=O.bodies[i.id1].shape.radius
rb=O.bodies[i.id2].shape.radius
vel1=np.hstack(O.bodies[i.id1].state.vel)
angVel1=np.hstack(O.bodies[i.id1].state.angVel)
vel2=np.hstack(O.bodies[i.id2].state.vel)
angVel2=np.hstack(O.bodies[i.id2].state.angVel)

## c++ code given by ScGeom.cpp
#alpha = (radius1+radius2)/(radius1+radius2-penetrationDepth);
#relativeVelocity = (rbp2->vel-rbp1->vel)*alpha + rbp2->angVel.cross(-radius2*normal) - rbp1->angVel.cross(radius1*normal);
#relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
#shearInc = relativeVelocity*scene->dt;

alpha=(ra+rb)/(ra+rb-p)
relVel=(vel2-vel1)*alpha+np.cross(rb*n,angVel2)-np.cross(ra*n,angVel1) ##getIncidentVel
relVel=relVel-np.dot(relVel,n)*n
shearInc=relVel*O.dt

shearInc ### =array([ -4.46043319e-07, -8.65634198e-07, -1.11027027e-07])
i.geom.shearInc ## =Vector3(3.662833531906604e-07,-1.6118379736082308e-06,2.7062622606841303e-09)
O.step()
i.geom.shearInc ## =Vector3(3.0053927015226938e-07,-1.5987564318912085e-06,-8.2929936709194397e-09)

Thanks in advance!
Best regards.
Lingran