#!/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(", \\")
