103 lines
3.2 KiB
Python
103 lines
3.2 KiB
Python
#!/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(", \\")
|