From 274461609b5d1f5bb272be9be95d0544ce4df31b Mon Sep 17 00:00:00 2001
From: Ian Jauslin <jauslin@ias.edu>
Date: Wed, 3 Oct 2018 17:58:38 +0000
Subject: As presented at Harvard University on 2018-10-03

---
 figs/atoms.fig/chiral.py | 102 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 figs/atoms.fig/chiral.py

(limited to 'figs/atoms.fig/chiral.py')

diff --git a/figs/atoms.fig/chiral.py b/figs/atoms.fig/chiral.py
new file mode 100644
index 0000000..eed6df8
--- /dev/null
+++ b/figs/atoms.fig/chiral.py
@@ -0,0 +1,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(", \\")
-- 
cgit v1.2.3-70-g09d2