```python
surf_choice = 'torus'
```


```python
from sympy import init_printing; init_printing();
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
```


```python
from silkpy.symbolic.surface.surface import ParametricSurface
from sympy import symbols, sin, cos, pi, cot, Array, refine, Q
from silkpy.sympy_utility import dot
u, v = symbols('u, v', real=True)
```


```python
if surf_choice=='cylindrical':
    R = symbols('R', positive=True)
    s = ParametricSurface([u, v], [R*cos(u), R*sin(u), v])
elif surf_choice=='cone':
    w = symbols('omega', real=True)
    s = ParametricSurface([u, v], [v*cos(u), v*sin(u), v*cot(w)])
elif surf_choice=='Mobius':
    theta = symbols('theta', real=True)
    s = ParametricSurface([theta, v],
        Array([cos(theta),                sin(theta),                0                ]) + 
        Array([sin(theta/2) * cos(theta), sin(theta/2) * sin(theta), cos(theta/2)]) * v)
elif surf_choice=='torus':
    from sympy import Q, ask
    from sympy.assumptions import global_assumptions
    a, r = symbols('a, r', real=True, positive=True)
    global_assumptions.add(Q.positive(a + r*cos(u)))
    s = ParametricSurface([u, v], [ (a+r*cos(u)) * cos(v), (a+r*cos(u)) * sin(v), r*sin(u)])
```


```python
s.christoffel_symbol.tensor()
```

    /home/wenyin/pyvenv/vsym/lib/python3.8/site-packages/IPython/lib/latextools.py:126: MatplotlibDeprecationWarning: 
    The to_png function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use mathtext.math_to_image instead.
      mt.to_png(f, s, fontsize=12, dpi=dpi, color=color)
    /home/wenyin/pyvenv/vsym/lib/python3.8/site-packages/IPython/lib/latextools.py:126: MatplotlibDeprecationWarning: 
    The to_rgba function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use mathtext.math_to_image instead.
      mt.to_png(f, s, fontsize=12, dpi=dpi, color=color)
    /home/wenyin/pyvenv/vsym/lib/python3.8/site-packages/IPython/lib/latextools.py:126: MatplotlibDeprecationWarning: 
    The to_mask function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use mathtext.math_to_image instead.
      mt.to_png(f, s, fontsize=12, dpi=dpi, color=color)
    /home/wenyin/pyvenv/vsym/lib/python3.8/site-packages/IPython/lib/latextools.py:126: MatplotlibDeprecationWarning: 
    The MathtextBackendBitmap class was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use mathtext.math_to_image instead.
      mt.to_png(f, s, fontsize=12, dpi=dpi, color=color)





$\displaystyle \left[\begin{matrix}\left[\begin{matrix}0 & 0\\0 & \frac{\left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{r}\end{matrix}\right] & \left[\begin{matrix}0 & - \frac{r \left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{a^{2} + 2 a r \cos{\left(u \right)} + r^{2} \cos^{2}{\left(u \right)}}\\- \frac{r \left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{a^{2} + 2 a r \cos{\left(u \right)} + r^{2} \cos^{2}{\left(u \right)}} & 0\end{matrix}\right]\end{matrix}\right]$




```python

```


```python
s.metric_tensor.tensor() 
s.metric_tensor.change_config('uu').tensor()
s.christoffel_symbol.tensor()
r_u, r_v = s.exprs.diff(u), s.exprs.diff(v); r_u, r_v
a_, b_ = r_u, r_v
s.weingarten_matrix
```




$\displaystyle \left[\begin{matrix}r^{2} & 0\\0 & \left(a + r \cos{\left(u \right)}\right)^{2}\end{matrix}\right]$






$\displaystyle \left[\begin{matrix}\frac{1}{r^{2}} & 0\\0 & \frac{1}{\left(a + r \cos{\left(u \right)}\right)^{2}}\end{matrix}\right]$






$\displaystyle \left[\begin{matrix}\left[\begin{matrix}0 & 0\\0 & \frac{\left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{r}\end{matrix}\right] & \left[\begin{matrix}0 & - \frac{r \left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{a^{2} + 2 a r \cos{\left(u \right)} + r^{2} \cos^{2}{\left(u \right)}}\\- \frac{r \left(a + r \cos{\left(u \right)}\right) \sin{\left(u \right)}}{a^{2} + 2 a r \cos{\left(u \right)} + r^{2} \cos^{2}{\left(u \right)}} & 0\end{matrix}\right]\end{matrix}\right]$






$\displaystyle \left( \left[\begin{matrix}- r \sin{\left(u \right)} \cos{\left(v \right)} & - r \sin{\left(u \right)} \sin{\left(v \right)} & r \cos{\left(u \right)}\end{matrix}\right], \  \left[\begin{matrix}- \left(a + r \cos{\left(u \right)}\right) \sin{\left(v \right)} & \left(a + r \cos{\left(u \right)}\right) \cos{\left(v \right)} & 0\end{matrix}\right]\right)$






$\displaystyle \left[\begin{matrix}\frac{1}{r} & 0\\0 & \frac{\cos{\left(u \right)}}{a + r \cos{\left(u \right)}}\end{matrix}\right]$




```python
Wa = s.weingarten_transform(a_) 
Wb = s.weingarten_transform(b_)
dot(Wa, b_), dot(a_, Wb)
s.K_H
s.prin_curvature_and_vector
from silkpy.sympy_utility import dot
(_, vec1), (_, vec2) = s.prin_curvature_and_vector
dot(vec1, vec2) # The two principal curvature vectors are perpendicular to each other.
```


```python
InteractiveShell.ast_node_interactivity = "last"
```


```python
from sympy import sin, cos, pi
from silkpy.numeric.surface.geodesic import geodesic_ncurve
theta = pi / 24 # symbols('theta', real=True)
t_arr, (u_arr, v_arr) = geodesic_ncurve(
    s.subs({a:5, r:2}), [pi/4, pi/4], [cos(theta), sin(theta)])
```


```python
from sympy import sin, cos, pi
from silkpy.numeric.surface.geodesic import geodesic_polar_ncoordinate
rho_arr, theta_arr, u_grid, v_grid = geodesic_polar_ncoordinate(
    s.subs({a:5, r:2}), [pi/4, pi/4], rho1=2.4, nrho=12, ntheta=48)

x_grid, y_grid, z_grid = s.subs({a:5, r:2}).lambdified()(u_grid, v_grid)
```


```python
from silkpy.symbolic.surface.draw import draw_surface_plotly
import plotly.graph_objects as go

if surf_choice=='torus':
    fig = draw_surface_plotly(s.subs({a: 5, r:2}), domain=[(-float(pi), float(pi)), (-float(pi), float(pi))])
else:
    raise NotImplementedError()
    
fig.add_trace(go.Scatter3d(
    x=x_arr, y=y_arr, z=z_arr,
    mode='lines',
    line=dict(color=t_arr, width=2)
))

# import numpy as np
# for i in range(len(theta_arr)):
#     fig.add_trace(go.Scatter3d(
#         x=x_grid[:, i], 
#         y=y_grid[:, i], 
#         z=z_grid[:, i],
#         mode='lines',
#         line=dict(#color=rho_arr, 
#                   width=2)
#     ))
    
# for i in range(len(rho_arr)):
#     fig.add_trace(go.Scatter3d(
#         x=np.r_[x_grid[i,:], x_grid[i,:]], 
#         y=np.r_[y_grid[i,:], y_grid[i,:]], 
#         z=np.r_[z_grid[i,:], z_grid[i,:]],
#         mode='lines',
#         line=dict(#color=rho_arr[i], 
#             width=2)
#     ))
# fig.show()
```

## Not yet done


```python
from sympy import series, Eq
t0 = symbols('t_0', real=True)
```


```python
t0 = 0
exprs[0].subs(t, t0) + (t-t0) * exprs[0].diff(t, 1).subs(t, t0)
exprs[1].subs(t, t0) + (t-t0) * exprs[1].diff(t, 1).subs(t, t0)
```


```python
exprs[0].evalf(subs={t:0}) + exprs[0].diff(t, 1).evalf(subs={t:0})
```


```python
from sympy import Eq
import sympy.solvers.ode as ode
ode.systems.dsolve_system([
    Eq(linearized_exprs[0], 0), 
    Eq(linearized_exprs[1], 0)], funcs=[u1, u2])
```


```python

```


```python
def curvature_curve(surface):
    from sympy import Matrix, Array, Eq
    from sympy import Function, symbols
    import sympy.solvers.ode as ode
    t = symbols('t', real=True)
    # u1, u2 = symbols('u1, u2', real=True, cls=Function)
    u1 = Function(surface.sym(0), real=True)(t)
    u2 = Function(surface.sym(1), real=True)(t)
    curvature_curve_mat = Matrix([
        [u1.diff(t)**2, -u1.diff(t) * u2.diff(t), u2.diff(t)**2], 
        Array(surface.E_F_G).subs(surface.sym(0), u1), 
        Array(surface.L_M_N).subs(surface.sym(1), u2)])
    # typically there would be two solutions
    sol_with_u1_equal_t = ode.systems.dsolve_system(
        [Eq(curvature_curve_mat.det(), 0 ), Eq(u1.diff(t), 1)])[0]
    sol_with_u2_equal_t = ode.systems.dsolve_system(
        [Eq(curvature_curve_mat.det(), 0 ), Eq(u2.diff(t), 1)])[0]
    return [sol_with_u1_equal_t, sol_with_u2_equal_t]
```


```python
curvature_curve(s)
```


```python

```

c76e18bae4d5b9c52142f939e4c5bb37c578fc2d
ipynb
Jupyter Notebook
nb/construct_surface.ipynb
WenyinWei/silkpy
2e3773e15fbebc4914f5627d8361a24d35fb5935
MIT
2021-05-13T09:23:03.000Z
2021-05-13T09:23:03.000Z
nb/construct_surface.ipynb
WenyinWei/silkpy
2e3773e15fbebc4914f5627d8361a24d35fb5935
MIT
nb/construct_surface.ipynb
WenyinWei/silkpy
2e3773e15fbebc4914f5627d8361a24d35fb5935
MIT
2021-05-27T14:13:24.000Z
2021-05-27T14:13:24.000Z
Qwen/Qwen-72B
1. YES
2. YES