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

# # Vector calculus with SageMath
# 
# ## Part 2: Using spherical coordinates
# 
# This notebook illustrates some vector calculus capabilities of SageMath within the 3-dimensional Euclidean space. The corresponding tools have been developed within
# the  [SageManifolds](https://sagemanifolds.obspm.fr) project.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Notebooks/SM_vector_calc_spherical.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, via the command `sage -n jupyter`

# *NB:* a version of SageMath at least equal to 8.3 is required to run this notebook:

# In[1]:


version()


# First we set up the notebook to display math formulas using LaTeX formatting:

# In[2]:


get_ipython().run_line_magic('display', 'latex')


# ## The 3-dimensional Euclidean space
# 
# We start by declaring the 3-dimensional Euclidean space $\mathbb{E}^3$, with $(r,\theta,\phi)$ as spherical coordinates:

# In[3]:


E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
print(E)
E


# $\mathbb{E}^3$ is endowed with the chart of spherical coordinates:

# In[4]:


E.atlas()


# as well as with the associated orthonormal vector frame $(e_r, e_\theta, e_\phi)$:

# In[5]:


E.frames()


# In the above output, $\left(\frac{\partial}{\partial r}, \frac{\partial}{\partial\theta}, \frac{\partial}{\partial \phi}\right)$ is the coordinate frame associated with $(r,\theta,\phi)$; it is not an orthonormal frame and will not be used below.

# ## Vector fields

# We define a vector field on $\mathbb{E}^3$ from its components in the orthonormal vector frame $(e_r,e_\theta,e_\phi)$:

# In[6]:


v = E.vector_field(r*sin(2*ph)*sin(th)^2 + r, 
                    r*sin(2*ph)*sin(th)*cos(th), 
                    2*r*cos(ph)^2*sin(th), name='v')
v.display()


# We can access to the components of $v$ via the square bracket operator:

# In[7]:


v[1]


# In[8]:


v[:]


# A vector field can evaluated at any point of $\mathbb{E}^3$:

# In[9]:


p = E((1, pi/2, pi), name='p')
print(p)


# In[10]:


p.coordinates()


# In[11]:


vp = v.at(p)
print(vp)


# In[12]:


vp.display()


# We may define a vector field with generic components:

# In[13]:


u = E.vector_field(function('u_r')(r,th,ph),
                   function('u_theta')(r,th,ph),
                   function('u_phi')(r,th,ph),
                   name='u')
u.display()


# In[14]:


u[:]


# In[15]:


up = u.at(p)
up.display()


# ## Algebraic operations on vector fields
# 
# ### Dot product
# 
# The dot (or scalar) product of the vector fields $u$ and $v$ is obtained by the method `dot_product`, which admits `dot` as a shortcut alias:

# In[16]:


s = u.dot(v)
s


# In[17]:


print(s)


# $s= u\cdot v$ is a *scalar field*, i.e. a map $\mathbb{E}^3 \rightarrow \mathbb{R}$:

# In[18]:


s.display()


# It maps points of $\mathbb{E}^3$ to real numbers:

# In[19]:


s(p)


# Its coordinate expression is

# In[20]:


s.expr()


# ### Norm
# 
# The norm of a vector field is 

# In[21]:


s = norm(u)
s


# In[22]:


s.display()


# In[23]:


s.expr()


# The norm is related to the dot product by $\|u\|^2 = u\cdot u$, as we can check:

# In[24]:


norm(u)^2 == u.dot(u)


# For $v$, we have

# In[25]:


norm(v).expr()


# ### Cross product
# 
# The cross product of $u$ by $v$ is obtained by the method `cross_product`, which admits `cross` as a shortcut alias:

# In[26]:


s = u.cross(v)
print(s)


# In[27]:


s.display()


# ### Scalar triple product
# 
# Let us introduce a third vector field. As a example, we do not pass the components as arguments of `vector_field`, as we did for $u$ and $v$; instead, we set them in a second stage, via the square bracket operator, any unset component being assumed to be zero:

# In[28]:


w = E.vector_field(name='w')
w[1] = r
w.display()


# The scalar triple product of the vector fields $u$, $v$ and $w$ is obtained as follows: 

# In[29]:


triple_product = E.scalar_triple_product()
s = triple_product(u, v, w)
print(s)


# In[30]:


s.expr()


# Let us check that the scalar triple product of $u$, $v$ and $w$ is $u\cdot(v\times w)$:

# In[31]:


s == u.dot(v.cross(w))


# ## Differential operators

# While the standard operators $\mathrm{grad}$, $\mathrm{div}$, $\mathrm{curl}$, etc. involved in vector calculus are accessible via the dot notation (e.g. `v.div()`), let us import functions `grad`, `div`, `curl`, etc. that allow for using standard mathematical notations
# (e.g. `div(v)`):

# In[32]:


from sage.manifolds.operators import *


# ### Gradient of a scalar field
# 
# We first introduce a scalar field, via its expression in terms of Cartesian coordinates; in this example, we consider a unspecified function of $(r,\theta,\phi)$:

# In[33]:


F = E.scalar_field(function('f')(r,th,ph), name='F')
F.display()


# The value of $F$ at a point:

# In[34]:


F(p)


# The gradient of $F$:

# In[35]:


print(grad(F))


# In[36]:


grad(F).display()


# In[37]:


norm(grad(F)).display()


# ### Divergence 
# 
# The divergence of a vector field:

# In[38]:


s = div(u)
s.display()


# In[39]:


s.expr().expand()


# For $v$ and $w$, we have

# In[40]:


div(v).expr()


# In[41]:


div(w).expr()


# An identity valid for any scalar field $F$ and any vector field $u$:

# In[42]:


div(F*u) == F*div(u) + u.dot(grad(F))


# ### Curl
# 
# The curl of a vector field:

# In[43]:


s = curl(u)
print(s)


# In[44]:


s.display()


# To use the notation `rot` instead of  `curl`, simply do

# In[45]:


rot = curl


# An alternative is

# In[46]:


from sage.manifolds.operators import curl as rot


# We have then

# In[47]:


rot(u).display()


# In[48]:


rot(u) == curl(u)


# For $v$ and $w$, we have

# In[49]:


curl(v).display()


# In[50]:


curl(w).display()


# The curl of a gradient is always zero:

# In[51]:


curl(grad(F)).display()


# The divergence of a curl is always zero:

# In[52]:


div(curl(u)).display()


# An identity valid for any scalar field $F$ and any vector field $u$:

# In[53]:


curl(F*u) == grad(F).cross(u) + F*curl(u)


# ### Laplacian

# The Laplacian of a scalar field:

# In[54]:


s = laplacian(F)
s.display()


# In[55]:


s.expr().expand()


# For a scalar field, the Laplacian is nothing but the divergence of the gradient:

# In[56]:


laplacian(F) == div(grad(F))


# The Laplacian of a vector field:

# In[57]:


Du = laplacian(u)
Du.display()


# Since this expression is quite lengthy, we may ask for a display component by component:

# In[58]:


Du.display_comp()


# We may expand each component:

# In[59]:


for i in E.irange():
    Du[i].expand()
Du.display_comp()


# In[60]:


Du[1]


# In[61]:


Du[2]


# In[62]:


Du[3]


# As a test, we may check that these formulas coincide with those of
# Wikipedia's article [*Del in cylindrical and spherical coordinates*](https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates#Del_formula).

# For $v$ and $w$, we have 

# In[63]:


laplacian(v).display()


# In[64]:


laplacian(w).display()


# We have

# In[65]:


curl(curl(u)).display()


# In[66]:


grad(div(u)).display()


# and we may check a famous identity:

# In[67]:


curl(curl(u)) == grad(div(u)) - laplacian(u)


# ## Customizations
# 
# ### Customizing the symbols of the orthonormal frame vectors
# 
# By default, the vectors of the orthonormal frame associated with spherical coordinates are denoted $(e_r,e_\theta,e_\phi)$:

# In[68]:


frame = E.spherical_frame()
frame


# But this can be changed, thanks to the method `set_name`:

# In[69]:


frame.set_name('a', indices=('r', 'th', 'ph'),
               latex_indices=('r', r'\theta', r'\phi'))
frame


# In[70]:


v.display()


# In[71]:


frame.set_name(('hr', 'hth', 'hph'), 
               latex_symbol=(r'\hat{r}', r'\hat{\theta}', r'\hat{\phi}'))
frame


# In[72]:


v.display()


# ### Customizing the coordinate symbols
# 
# The coordinates symbols are defined within the angle brackets `<...>` at the construction of the Euclidean space. Above we did 

# In[73]:


E.<r,th,ph> = EuclideanSpace(coordinates='spherical')


# which resulted in the coordinate symbols $(r,\theta,\phi)$ and in the corresponding Python variables `r`, `th` and `ph` (SageMath symbolic expressions). Using other symbols, for instance $(R,\Theta,\Phi)$, is possible through the optional argument `symbols` of the function `EuclideanSpace`. It has to be a string, usually prefixed by r (for raw string, in order to allow for the backslash character of LaTeX expressions). This string contains the coordinate fields separated by a blank space; each field contains the coordinate’s text symbol and possibly the coordinate’s LaTeX symbol (when the latter is different from the text symbol), both symbols being separated by a colon (`:`):

# In[74]:


E.<R,Th,Ph> = EuclideanSpace(coordinates='spherical', symbols=r'R Th:\Theta Ph:\Phi')


# We have then

# In[75]:


E.atlas()


# In[76]:


E.frames()


# In[77]:


E.spherical_frame()


# In[78]:


v = E.vector_field(R*sin(2*Ph)*sin(Th)^2 + R, 
                    R*sin(2*Ph)*sin(Th)*cos(Th), 
                    2*R*cos(Ph)^2*sin(Th), name='v')
v.display()