yade-users team mailing list archive

[Jan Stránský <question691560@xxxxxxxxxxxxxxxxxxxxx>]

Question #691560 on Yade changed:
https://answers.launchpad.net/yade/+question/691560

    Status: Open => Answered

Jan Stránský proposed the following answer:
Hello,

thanks for the code. I have tried it and the results are IMHO nicely isotropic.
Evaluation:
I chose 200 random unit vectors and for these vectors I counted normals that are inclined from the vector maximally by given angle.
>From these counts I computed relative standard deviation, whose value is:
1.1 % for angle 30 deg
3.2 % for angle 15 deg
11 % for angle 5 deg
Plotting the values gives relatively nice semi-sphere.
Of course it does not give the same number for every direction, but I did not find any systematic problem ("there is little contact in the vertical direction, and the middle is recessed")
###
# Extract normals from file
with open("normals.txt") as f:
   lines = f.readlines()
lines = [l.split()[2:] for l in lines]
lines = [[float(v) for v in l] for l in lines]
normals = [Vector3(*l) for l in lines]
print("Number of normals:",len(normals))

# Function to test isotropy.
# For random direction counts normals "close" to this direction
def countSimilarNormals(limAngle=pi/6):
   # random unit vector with positive z coordinate
   from random import gauss
   v = [gauss(0, 1) for i in range(3)]
   mag = sqrt(sum(pow(x,2) for x in v))
   v = Vector3(*[x/mag for x in v])
   if v[2] < 0: v *= -1
   # counting "similar" normals
   ret = 0
   for n in normals:
      # compute angle
      dot = v.dot(n)
      dot = abs(dot)
      angle = acos(dot)
      # if it is within the limit, count it
      if angle <= limAngle:
         ret += 1
   return v,ret

# count similar normals for a few random directions
data = [countSimilarNormals(pi/6/2/3) for i in range(200)]
directions,counts = zip(*data)
# some statistics
cmax = max(counts)
import statistics
cmean = statistics.mean(counts)
cstddev = statistics.stdev(counts)
print("Isotropy fulfilled with",cstddev/cmean*100,"% relative standard deviation")

# plot
dvals = [d*c/cmax for (d,c) in zip(directions,counts)]
x,y,z = [[d[i] for d in dvals] for i in (0,1,2)]
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection="3d")
ax.auto_scale_xyz((-1,1),(-1,1),(0,1))
ax.scatter3D(x,y,z)
plt.show()
# sorry, I don't know how to set aspect ratio, so the result is displayed as semi-ellipsoid
###

> I plotted the spherical histogram using Matlab.
> Am I doing something wrong here?

please provide the code. As discussed, this stage may contain a mistake.

> I couldn't attach my histogram plot here

yes, thanks for respecting the rule, but you can add the data (i.e. a
few lines of e.g. angle-number pairs) or the code producing the plot or
the data.

cheers
Jan

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.