from math import sin,cos,pi,sqrt
from pylab import plot,show,axis,savefig


# positions and orientations of the segments
L,r=1.,2.
alpha1,betha1=pi/2,0
alpha2,betha2=0,0
dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)


def poly1(a,b,c):
 if a==0:
  if b<>0:
   return -c/(2*b)
  else:
   return 'N'
 D=b**2-a*c
 if D<0:
  return 'N'
 else:
  l=(-b-sqrt(D))/a
  if l<0:
   return 'N'
  else:
   return l

def poly2(a,b,c):
 if a==0:
  return 'N'
 D=b**2-a*c
 if D<0:
  return 'N'
 else:
  l=(-b+sqrt(D))/a
  if l<0:
   return 'N'
  else:
   return l



def voroll1(L,r,alpha1,betha1,alpha2,betha2,theta,phi):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2

 a=sd1**2-sd2**2
 b=dy2*r*sd2-r*sy
 c=(1-dy2**2)*r**2
 l=poly1(a,b,c)
 if l=='N':
  return 'N'
 else:
   t1,t2=sd1*l,sd2*l-dy2*r
   if (t1>-L and t1<L and t2>-L and t2<L):
    return l
   else:
    return 'N'

def voroll2(L,r,alpha1,betha1,alpha2,betha2,theta,phi):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2

 a=sd1**2-sd2**2
 b=dy2*r*sd2-r*sy
 c=(1-dy2**2)*r**2
 l=poly2(a,b,c)
 if l=='N':
  return 'N'
 else:
   t1,t2=sd1*l,sd2*l-dy2*r
   if (t1>-L and t1<L and t2>-L and t2<L):
    return l
   else:
    return 'N'


def vorolp1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,t2):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2
   
 a=sd1**2
 b=-sd2*t2-r*sy
 c=t2**2+2*dy2*t2*r+r**2
 if a==0:
  return 'N'
 l=poly1(a,b,c)
 if l=='N':
  return l
 else:
  DD2=(sd2*l-dy2*r)*t2
  t1,t2=sd1*l,sd2*l-dy2*r
  if (t1>-L and t1<L) and (t2<-L or t2>L) and DD2>0:
   return l
  else:
   return 'N'

def vorolp2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,t2):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2
   
 a=sd1**2
 b=-sd2*t2-r*sy
 c=t2**2+2*dy2*t2*r+r**2
 l=poly2(a,b,c)
 if l=='N':
  return l
 else:
  DD2=(sd2*l-dy2*r)*t2
  t1,t2=sd1*l,sd2*l-dy2*r
  if (t1>-L and t1<L) and (t2<-L or t2>L) and DD2>0:
   return l
  else:
   return 'N'


def voropl1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,t1):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2
   
 a=sd2**2
 b=-(sd1*t1+dy2*r*sd2-r*sy)
 c=t1**2-(1-dy2**2)*r**2
 l=poly1(a,b,c)
 if l=='N':
  return l
 else:
  DD1=sd1*t1
  t1,t2=sd1*l,sd2*l-dy2*r
  if (t1<-L or t1>L) and (t2>-L and t2<L) and DD1>0:
   return l
  else:
   return 'N'

def voropl2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,t1):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2
   
 a=sd2**2
 b=-(sd1*t1+dy2*r*sd2-r*sy)
 c=t1**2-(1-dy2**2)*r**2
 l=poly2(a,b,c)
 if l=='N':
  return l
 else:
  DD1=sd1*t1
  t1,t2=sd1*l,sd2*l-dy2*r
  if (t1<-L or t1>L) and (t2>-L and t2<L) and DD1>0:
   return l
  else:
   return 'N'

def voropp(L,r,alpha1,betha1,alpha2,betha2,theta,phi,t1,t2):
 dx1,dy1,dz1=sin(alpha1)*cos(betha1),sin(alpha1)*sin(betha1),cos(alpha1)
 dx2,dy2,dz2=sin(alpha2)*cos(betha2),sin(alpha2)*sin(betha2),cos(alpha2)
 sx,sy,sz=sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)
 sd1,sd2=sx*dx1+sy*dy1+sz*dz1,sx*dx2+sy*dy2+sz*dz2

 k=2*(sd2*t2-sd1*t1+r*sy)
 if k==0:
  return 'N'
 else:
  l=(t2**2-t1**2+2*dy2*t2*r+r**2)/k
 if l<0:
  return 'N'
 else:
  DD1,DD2=sd1*t1,(sd2*l-dy2*r)*t2
  t1,t2=sd1*l,sd2*l-dy2*r
  if (t1<-L or t1>L) and (t2<-L or t2>L) and (DD1>0 and DD2>0):
   return l
  else:
   return 'N'


N=100
X,Y,Z=[],[],[]
for j in range(N):
 for i in range(N):
  theta=(2*pi*i/(N-1.))
  phi=(pi*j/(N-1.))
  l=voroll1(L,r,alpha1,betha1,alpha2,betha2,theta,phi)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voroll2(L,r,alpha1,betha1,alpha2,betha2,theta,phi)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=vorolp1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=vorolp1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=vorolp2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=vorolp2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropl1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropl1(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropl2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropl2(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropp(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropp(L,r,alpha1,betha1,alpha2,betha2,theta,phi,-L,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropp(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L,-L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))
  l=voropp(L,r,alpha1,betha1,alpha2,betha2,theta,phi,L,L)
  if l<>'N':
   X.append(l*sin(phi)*cos(theta))
   Y.append(l*sin(phi)*sin(theta))
   Z.append(l*cos(phi))

f=open('3Dlim.txt','w') #boundary to gnuplot (splot)
for i in range(len(X)):
 print >>f, X[i],Y[i],Z[i]
f.close()


X,Y,Z=[],[],[]
for i in range(11):
 t=(2*i-5)*L/10.
 X.append(dx1*t)
 Y.append(dy1*t)
 Z.append(dz1*t)
 X.append(dx2*t)
 Y.append(dy2*t+r)
 Z.append(dz2*t)


f=open('3Dline.txt','w') #segments to gnuplot (splot)
for i in range(len(X)):
 print >>f, X[i],Y[i],Z[i]
f.close()

