#!/usr/bin/env python
# coding: utf-8

# ## Tubular surfaces ##

# A tubular surface (or tube surface) is generated by a 3D curve, called spine, and a moving circle of radius r, with center on the spine and  included in planes orthogonal to curve.
# 
# Tubular surfaces are associated to spines that are biregular, that is, they have a $C^2$ parameterization, $c:[a,b]\to \mathbb{R}^3$, with 
# velocity, $\dot{c}(t)$, and  acceleration, $\ddot{c}(t)$, that are  non-null and non-colinear vectors:
# $\dot{c}(t)\times \ddot{c}(t)\neq 0$.

# ### Tubular surface defined by a spine curve parameterized by arc length ###

# A tube of prescribed [curvature](https://en.wikipedia.org/wiki/Curvature#Curvature_of_space_curves) and [torsion](https://en.wikipedia.org/wiki/Torsion_of_a_curve) is defined by a spine parameterized by the arc length, i.e. by
# $c(s)$, with constant speed, $||\dot{c}(s)||=1$, and non-null acceleration, $\ddot{c}(s)\neq 0$, for all $s$.
# 
# The given curvature and torsion, $\kappa(s)$,  $\tau(s)$, define the Frenet-Serre equations:
# $$\begin{array}{lll}
#     \dot{e}_1(s)&=&\kappa(t)e_2(s)\\
#     \dot{e}_2(s)&=&-\kappa(s)e_1(s)+\tau(s)e_3(s)\\
#     \dot{e}_3(s)&=&-\tau(s)e_2(s),\\
#     \end{array} $$
#     
# where $e_1(s), e_2(s), e_3(s)$ are  respectively the unit vectors of tangent, principal normal and binormal along the curve.
# 
# Frenet-Serre equations completed with the equation $ \dot{c}(s)=e_1(s)$ define a system of ordinary differential equations, with  12 equations and 12 unknown functions. The last three
# coordinates of a solution represent the  discretized curve,  $c(s)$, starting from an initial point, with a prescribed Frenet frame at that point.

# We define below a tubular surface with highly oscillating curvature and constant torsion of the spine.

# In[1]:


import numpy as np
from scipy import integrate


# In[2]:


def curv(s):#curvature
    return   3*np.sin(s/10.)*np.sin(s/10.)

def tors(s):#torsion is constant
    return 0.35 

def Frenet_eqns(x, s):# right side vector field of the system of ODE
    return [ curv(s)*x[3],
             curv(s)*x[4], 
             curv(s)*x[5], 
            -curv(s)*x[0]+tors(s)*x[6], 
            -curv(s)*x[1]+tors(s)*x[7], 
            -curv(s)*x[2]+tors(s)*x[8],
            -tors(s)*x[3], 
            -tors(s)*x[4],
            -tors(s)*x[5],
             x[0], x[1], x[2]]


# Integrate the system, with an initial point consisting in the initial Frenet frame (of three orthonormal vectors)
# and the initial position of the curve,  $c(0)$:

# In[3]:


x_init=np.array([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0])
                                                     
s_final=150# [0, s_final] is the interval of integration
N=1000
s_div=np.linspace(0, s_final, N)  
X=integrate.odeint(Frenet_eqns, x_init, s_div)

normal=X[:, 3:6].T
binormal=X[:, 6:9].T
curve=X[:, 9:].T
xc, yc, zc=curve# lists of coordinates of the spine  points 


# Now we define a tubular surface that has as spine the above curve.

# 
# A  tubular surface having as spine  a curve, $c(s)$, parameterized by the arclength, is defined as follows:
# $r(s,u)=c(s)+\varepsilon(e_2(s)cos(u)+e_3(s)sin(u))$, $0<\varepsilon <<1$, $u\in[0, 2\pi]$.
# $\varepsilon$ is the radius of circles orthogonal to the spine.

# In[4]:


import plotly.plotly as py
from plotly.graph_objs import *


# Define a function that sets the plot layout:

# In[9]:


axis = dict(
showbackground=True, 
backgroundcolor="rgb(230, 230,230)",
gridcolor="rgb(255, 255, 255)",      
zerolinecolor="rgb(255, 255, 255)",  
          )


noaxis=dict(showbackground=False,
            showgrid=False,
            showline=False,
            showticklabels=False,
            ticks='',
            title='',
            zeroline=False
           )


def set_layout(title='', width=800, height=800, axis_type=axis, aspect=(1, 1, 1)):
    return Layout(
                 title=title,
                 autosize=False,
                 width=width,
                 height=height,
                 showlegend=False,
                 scene=Scene(xaxis=XAxis(axis_type),
                             yaxis=YAxis(axis_type),
                             zaxis=ZAxis(axis_type), 
                            
                             aspectratio=dict(x=aspect[0],
                                              y=aspect[1],
                                              z=aspect[2]
                                             )
                             ) 
                 )   


# The colorscale for the tubular surface:

# In[10]:


my_colorscale=[[0.0, 'rgb(46, 107, 142)'],
               [0.1, 'rgb(41, 121, 142)'],
               [0.2, 'rgb(36, 134, 141)'],
               [0.3, 'rgb(31, 147, 139)'],
               [0.4, 'rgb(30, 160, 135)'],
               [0.5, 'rgb(40, 174, 127)'],
               [0.6, 'rgb(59, 186, 117)'],
               [0.7, 'rgb(85, 198, 102)'],
               [0.8, 'rgb(116, 208, 84)'],
               [0.9, 'rgb(151, 216, 62)'],
               [1.0, 'rgb(189, 222, 38)']]


# Define a function that evaluates the tube parameterization,  $r(s,u)=(x, y, z)$, at the meshgrid `np.meshgrid(s_div, u)`:

# In[13]:


def create_tube(spine_points,  normal, binormal,
                epsilon=0.2, colorscale=my_colorscale, zmin=None, zmax=None):
    #returns an instance of the Plotly Surface, representing a tube 
    
    u=np.linspace(0, 2*np.pi, 100)
    x,y,z=[np.outer(spine_points[k,:], np.ones(u.shape))+
          epsilon*(np.outer(normal[k, :], np.cos(u))+np.outer(binormal[k,:], np.sin(u))) 
          for k in range(3)]
   
    if zmin is not None and zmax is not None:
        return Surface(x=x, y=y, z=z, zmin=zmin, zmax=zmax,
                   colorscale=colorscale, 
                   colorbar=dict(thickness=25, lenmode='fraction', len=0.75)) 
    else:
        return Surface(x=x, y=y, z=z, 
                   colorscale=colorscale, 
                   colorbar=dict(thickness=25, lenmode='fraction', len=0.75)) 


# The keywords `zmin`, `zmax` are set when we connect at least two tubular surfaces. They define the color bounds for
# the tubular structure.

# In[14]:


tube=create_tube(curve,  normal, binormal, epsilon=0.1)


# In[15]:


data1=Data([tube])
layout1=set_layout(title='Tubular surface', aspect=(1,1,1.05))
fig1 = Figure(data=data1, layout=layout1)

py.sign_in('empet', '')
py.iplot(fig1, filename='tubular-cst-torsion')


# ### Tubular surface with a spine curve of given  parameterization ###

# If a general  biregular  parameterization,  $c(t)$, of the spine  is given,
# then we have to do some analytical computations by hand, in order to get the
# directions $\dot{c}(t)$, $\ddot{c}(t)$, $\dot{c}(t)\times \ddot{c}(t)$, of the velocity (tangent), acceleration, and binormals along the curve.
# 
# Then we define  Python functions, `tangent`, `acceleration`,  `curve_normals`, that compute the unit vectors of these directions.
# Finally the unit vector of the principal normal is computed as $n(t)=b(t)\times tg(t)$, where $b(t), tg(t)$ are the unit vectors of binormals and tangents.
#  

# The tube parameterization, $$r(t,u)=c(t)+\varepsilon(n(t)\cos(u)+b(t)\sin(u)), t\in[tm, tM], u\in[0,2\pi],$$
# is evaluated at a meshgrid.

# We illustrate a tubular structure, called [Hopf link](https://en.wikipedia.org/wiki/Hopf_link), defined by two tubes, having the spines parameterized by:
#     $$c(t)=(\pm a+\cos(t), \sin(t), \pm b\sin(t)), t\in[0, 2\pi]$$
# The first spine corresponds to $a=0.5, b=0.2$, and the second one, to $a=-0.5, b=-0.2$.

# In[16]:


from numpy import sin, cos, pi 


# In[17]:


def spine_param( a, b, tm, tM, nr): 
    #spine parameterization c:[tm, tM]-->R^3 
    # a, b are parameters on which  the spine parameterization depends
    # nr is the number of points to be evaluated on spine
    t=np.linspace(tm, tM, nr )# nr is the number of points to ve evaluated on spine
    return t,  a+cos(t), sin(t), b*sin(t)


# In[18]:


def tangent( a, b, t):
    # returns the unit tangent vectors along the spine curve
    v=np.vstack((-sin(t), cos(t), b*cos(t)))
    return v/np.vstack((np.linalg.norm(v, axis=0),)*3)

def acceleration( a, b, t):
    # returns the unit acceleration  vectors along the spine
    v=np.array([ -cos(t), -sin(t), -b*sin(t)])
    return v/np.vstack((np.linalg.norm(v, axis=0),)*3)


# In[19]:


def curve_normals(a, b): 
    # computes and returns  the point coordinates on spine, and the unit normal vectors
    
    t, xc, yc, zc=spine_param(a,b, 0.0, 2*pi, 100)
    tang=tangent(a,b, t) 
    binormal=np.cross(tang, acceleration(a, b, t), axis=0)
    binormal=binormal/np.vstack((np.linalg.norm(binormal, axis=0),)*3)
    normal=np.cross(binormal, tang, axis=0)
    return np.vstack((xc, yc, zc)), normal, binormal


# In[20]:


epsilon=0.025 # the radius of each tube
zm=[]# list of min z-values on both tubes
zM=[]# list of max z-values on both tubes
spine1, normal1, binormal1=curve_normals(0.5, 0.2)
zm.append(min(spine1[2,:]))
zM.append(max(spine1[2,:]))
spine2, normal2, binormal2=curve_normals(-0.5, -0.2)
zm.append(min(spine2[2,:]))
zM.append(max(spine2[2,:]))

zmin=min(zm)
zmax=max(zM)

tube1=create_tube(spine1, normal1, binormal1, epsilon=epsilon, zmin=zmin, zmax=zmax)
tube2=create_tube(spine2, normal2, binormal2, epsilon=epsilon, zmin=zmin, zmax=zmax)
layout2=set_layout(title='Hopf link',  aspect=(1, 0.75, 0.35))


# In[21]:


data2=Data([tube1,tube2])
fig2 = Figure(data=data2, layout=layout2)


# In[22]:


py.sign_in('empet', '')
py.iplot(fig2, filename='Hopf-link')


# If we take all combinations of signs for the parameters, a,  b, we get an interesting configuration of tubes
# communicating with each other:

# In[23]:


from IPython.display import HTML
HTML('<iframe src=https://plot.ly/~empet/13930/comunicating-rings/ width=900 height=700></iframe>')


# ###  Canal (Channels) surfaces ###

# Tubular surfaces are particular surfaces in the class of  canal surfaces. A canal surface
# is again defined by a biregular spine, $c(t)$, but the circles
# ortogonal to spine have variable radii,  gigen by a $C^1$-function, $r(t)$, with $|r'(t)|<||\dot{c}(t)||$.
# 
# The parameterization of a canal surface is:
#     
# $$r(t,u)=c(t)-\displaystyle\frac{r(t)r'(t)}{||\dot{c}(t)||^2}\dot{c}(t)+
#     \displaystyle\frac{r(t)\sqrt{||\dot{c}(t) ||^2-r'(t)^2}}{||\dot{c}(t) ||}(n(t)\cos(u)+b(t)\sin(u))$$

# We plot the canal surface of spine, $c(t)=(10\cos(t),  10\sin(t), 0)$, and  radius function
# $r(t)=2+\cos(2t)$, $t\in[0,2\pi]$.

# In[24]:


def radius_deriv(t):
    return 2+cos(2*t), -2*sin(2*t)


# In[25]:


def create_canal(spine,  normal, binormal, term,
                 colorscale=my_colorscale, zmin=None, zmax=None):
    #returns an instance of the Plotly Surface, representing a  canal surface
    #term is the second term in the parameterization
    
    u=np.linspace(0, 2*np.pi, 100)
    x,y,z=[np.outer(spine[k,:]-term[k, :], np.ones(u.shape))+\
           np.outer(normal[k, :], np.cos(u))+np.outer(binorm[k,:], np.sin(u)) for k in range(3)]
     
        
    if zmin is not None and zmax is not None:
        return Surface(x=x, y=y, z=z, zmin=zmin, zmax=zmax,
                   colorscale=colorscale, 
                   colorbar=dict(thickness=25, lenmode='fraction', len=0.75)) 
    else:
        return Surface(x=x, y=y, z=z, 
                   colorscale=colorscale, 
                   colorbar=dict(thickness=25, lenmode='fraction', len=0.75)) 


# In[26]:


t=np.linspace(0, 3*pi/2, 50)
xc, yc, zc= 10*cos(t), 10*sin(t), np.zeros(t.shape)
spine=np.vstack((xc,yc, zc))
rt,rdt=radius_deriv(t)# rt is the variable radius r(t), and rdt its derivative


# In[27]:


tang=np.vstack((-10*sin(t), 10*cos(t),  np.zeros(t.shape))) #c'(t)
cdot_norm=np.vstack((np.linalg.norm(tang, axis=0),)*3)# ||c'(t)||
 
factor=rt*rdt/cdot_norm**2
term=factor*tang#term.shape=(3, t.shape[0])# second term in canal surface parameterization
R=rt*np.sqrt(cdot_norm**2-rdt**2)/cdot_norm  # R.shape (3, t.shape[0])  is the scalar factor in the third term

tangu= (tang/cdot_norm) #unit tangent vector
acceler=np.vstack((-10*cos(t), -10*sin(t),  np.zeros(t.shape)))
acceler= acceler/np.vstack((np.linalg.norm(acceler, axis=0),)*3)#unit acceleration vector
binorm=np.cross(tangu, acceler, axis=0)
binorm=binorm/np.vstack((np.linalg.norm(binorm, axis=0),)*3)#unit binormal vector
normal=np.cross(binorm, tangu, axis=0)# unit normal vector
binorm=R*binorm
normal=R*normal


# In[28]:


canal=create_canal(spine,  normal, binorm, term,  colorscale=my_colorscale)


# In[29]:


layout3=set_layout(title='Canal surface', axis_type=axis, aspect=(1, 1, 0.25))
data3=Data([canal])
fig3 = Figure(data=data3, layout=layout3) 

py.sign_in('empet', '')
py.iplot(fig3, filename='Canal-surf')


# Finally, we stress that in order to get a   tubular looking surface, we have to set the aspect ratio 
# of the plot that respects the real ratios between axes lengths. Otherwise the tube is deformed.

# In[31]:


from IPython.core.display import HTML
def  css_styling():
    styles = open("./custom.css", "r").read()
    return HTML(styles)
css_styling()