1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
#!/usr/bin/env python3
from math import *
import random
import sys
# size of space
L=30
# number of lines
H=5
# heigh of lines
height=5
# number of rods per line
N=16
# aspect ratio
a=10
# spread in theta angle
spread_t=pi/60
# in phi
spread_p=pi/60
# in height
spread_h=1/120
# check whether two rods overlap
def check_overlap(rod1,rod2):
# relative placement
relative_pos=unrotate(subtract(rod2[0],rod1[0]), rod1[1])
if(abs(relative_pos[0])<2 and abs(relative_pos[1])<2 and abs(relative_pos[2])<2):
return(True)
# relative angle
relative_ang=cart_to_spherical(unrotate(spherical_to_cart(rod2[1]), rod1[1]))
# exclusion volume
# rotate other rod
relative_pos=unrotate(relative_pos, [0,relative_ang[1]])
#if(abs(relative_pos[1])<2 and abs(relative_pos[0])-2<abs(sin(relative_ang[0]))*a and abs(relative_pos[2])-a-2<abs(cos(relative_ang[0])*a)):
if(abs(relative_pos[1])<2 and abs(relative_pos[0])-2<abs(sin(relative_ang[0]))*a and abs(sin(relative_ang[0])*relative_pos[2]-cos(relative_ang[0])*relative_pos[0])-2<abs(sin(relative_ang[0])*a)):
return(True)
return(False)
# def subtract vectors
def subtract(x,y):
return([x[0]-y[0],x[1]-y[1],x[2]-y[2]])
# rotate vector
def unrotate(x,w):
ret=[x[0],x[1],x[2]]
# rotate phi
tmp=cos(w[1])*ret[0]+sin(w[1])*ret[1]
ret[1]=-sin(w[1])*ret[0]+cos(w[1])*ret[1]
ret[0]=tmp
# rotate theta
tmp=cos(w[0])*ret[0]-sin(w[0])*ret[2]
ret[2]=sin(w[0])*ret[0]+cos(w[0])*ret[2]
ret[0]=tmp
return(ret)
# convert coordinates
def spherical_to_cart(w):
return([cos(w[1])*sin(w[0]),sin(w[1])*sin(w[0]),cos(w[0])])
def cart_to_spherical(x):
w=[0,0]
w[0]=acos(x[2])
if(sin(w[0]==0)):
return([w[0],0])
c=x[0]/sin(w[0])
s=x[1]/sin(w[0])
# to avoid truncation errors
if(abs(c)>1 and abs(c)<1.0001):
if(c>0):
return([w[0],0])
else:
return([w[0],pi])
if(s>=0):
return([w[0],acos(c)])
return([w[0],2*pi-acos(c)])
# configuration
config=[]
# add rods
for h in range(H):
config_h=[]
while len(config_h)<N:
# random position and angles
x=[random.uniform(0,L), random.uniform(0,L), random.gauss(h*height,spread_h)]
w=[abs(random.gauss(pi/2,spread_t)), random.gauss(pi/2*(1-h/(H-1)),spread_p)]
# chek it does not interfere with other rods
fine=True
for rod in config+config_h:
if(check_overlap(rod,[x,w])):
fine=False
break
if fine:
config_h.append([x,w])
config=config+config_h
for i in range(len(config)):
rod=config[i]
print(str(rod[0][0])+"+("+str(cos(rod[1][1])*cos(rod[1][0]))+")*cos(u)*sin(v)+("+str(-sin(rod[1][1]))+")*sin(u)*sin(v)+("+str(cos(rod[1][1])*sin(rod[1][0])*a)+")*cos(v)", end=", ")
print(str(rod[0][1])+"+("+str(sin(rod[1][1])*cos(rod[1][0]))+")*cos(u)*sin(v)+("+str(cos(rod[1][1]))+")*sin(u)*sin(v)+("+str(sin(rod[1][1])*sin(rod[1][0])*a)+")*cos(v)", end=", ")
print(str(rod[0][2])+"+("+str(-sin(rod[1][0]))+")*cos(u)*sin(v)+("+str(cos(rod[1][0])*a)+")*cos(v)", end=" ")
print("with pm3d", end="")
if i<len(config)-1:
print(", \\")
|
