In [1]:
import numpy as np
import pandas as pd
from lets_plot import *

LetsPlot.setup_html()
In [2]:
def bottle(theta, phi):
    x = -2/15*np.cos(theta)*(3*np.cos(phi) - 
                         30*np.sin(theta) + 
                         90*np.cos(theta)**4*np.sin(theta) - 
                         60*np.cos(theta)**6*np.sin(theta) + 
                         5*np.cos(theta)*np.cos(phi)*np.sin(theta))
    y = -1/15*np.sin(theta)*(3*np.cos(phi) - 
                         3*np.cos(theta)**2*np.cos(phi) - 
                         48*np.cos(theta)**4*np.cos(phi) + 
                         48*np.cos(theta)**6*np.cos(phi) - 60*np.sin(theta) + 
                         5*np.cos(theta)*np.cos(phi)*np.sin(theta) -
                         5*np.cos(theta)**3*np.cos(phi)*np.sin(theta)-
                         80*np.cos(theta)**5*np.cos(phi)*np.sin(theta)+
                         80*np.cos(theta)**7*np.cos(phi)*np.sin(theta))
    z = 2/15*(3 + 5*np.cos(theta)*np.sin(theta))*np.sin(phi)
    return x, y, z
In [3]:
n = 100
In [4]:
theta = np.linspace(0, np.pi, n)
phi = np.linspace(0, 2*np.pi, n)
In [5]:
thetas, phis = np.meshgrid(theta, phi)
theta, phi = thetas.reshape(-1), phis.reshape(-1)
In [6]:
x, y, z = bottle(theta, phi)
X = np.vstack((x, y, z)).T
In [7]:
def rotation_matrices(theta_x, theta_y, theta_z):
    Rx = [[1, 0, 0],
      [0, np.cos(theta_x), -np.sin(theta_x)],
      [0, np.sin(theta_x), np.cos(theta_x)]]
    Ry = [[np.cos(theta_y), 0, np.sin(theta_y)], 
      [0, 1, 0], 
      [-np.sin(theta_y), 0, np.cos(theta_y)]]
    Rz = [[np.cos(theta_z), -np.sin(theta_z), 0], 
      [np.sin(theta_z), np.cos(theta_z), 0], 
      [0, 0, 1]]
    return Rx, Ry, Rz
In [8]:
Rx, Ry, Rz = rotation_matrices(np.pi/6, np.pi/6, np.pi/6)
In [9]:
X_rotated = ((X.dot(Rx)).dot(Ry)).dot(Rz)
In [10]:
data = pd.DataFrame(X_rotated, columns=['x', 'y', 'z'])
data['cross-sectional'] = np.tile(np.arange(n), n)
data['longitudinal'] = np.repeat(np.arange(n), n)
In [11]:
p = ggplot(data) \
        + geom_path(aes(x='x', y='y', group='longitudinal'), alpha=0.5, color='gray', sampling='none')
p += theme(axis_text='blank', axis_ticks='blank', axis_line='blank', axis_title='blank')
p
Out[11]:
In [12]:
p += geom_path(aes(x='x', y='y', group='cross-sectional'), alpha=0.2, color='gray', sampling='none')
p
Out[12]: