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

# # How to compute a gradient, a divergence or a curl?
# 
# 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/VectorCalculus/vector_calc_cartesian.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')


# ## First stage: introduce the Euclidean 3-space
# 
# Before evaluating some vector-field operators, one needs to define the arena
# in which vector fields live, namely the *3-dimensional Euclidean space*
# $\mathbb{E}^3$. In SageMath, we declare it, along with the standard
# *Cartesian coordinates* $(x,y,z)$, as follows:

# In[3]:


E.<x,y,z> = EuclideanSpace()
print(E)
E


# Thanks to the notation `<x,y,z>` in the above declaration, the coordinates
# $(x,y,z)$ are immediately available as three symbolic variables ``x``,
# ``y`` and ``z`` (there is no need to type ``x, y, z = var('x y z')``):

# In[4]:


x is E.cartesian_coordinates()[1]


# In[5]:


y is E.cartesian_coordinates()[2]


# In[6]:


z is E.cartesian_coordinates()[3]


# In[7]:


type(z)


# Besides, $\mathbb{E}^3$ is endowed with the *orthonormal vector frame*
# $(e_x, e_y, e_z)$ associated with Cartesian coordinates:

# In[8]:


E.frames()


# At this stage, this is the *default* vector frame on $\mathbb{E}^3$
# (being the only vector frame introduced so far):

# In[9]:


E.default_frame()


# ## Defining a vector field

# We define a vector field on $\mathbb{E}^3$ from its components in the vector frame $(e_x,e_y,e_z)$:

# In[10]:


v = E.vector_field(-y, x, sin(x*y*z), name='v')
v.display()


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

# In[11]:


v[1]


# In[12]:


v[:]


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

# In[13]:


p = E((3,-2,1), name='p')
print(p)


# In[14]:


p.coordinates()


# In[15]:


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


# In[16]:


vp.display()


# Vector fields can be plotted (see the [list of options](http://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/differentiable/vectorfield.html#sage.manifolds.differentiable.vectorfield.VectorField.plot) for customizing the plot)

# In[17]:


v.plot(max_range=1.5, scale=0.5)


# We may define a vector field with generic components $(u_x, u_y, u_z)$:

# In[18]:


u = E.vector_field(function('u_x')(x,y,z),
                   function('u_y')(x,y,z),
                   function('u_z')(x,y,z),
                   name='u')
u.display()


# In[19]:


u[:]


# Its value at the point $p$ is then

# In[20]:


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


# ## How to compute various vector products?
# 
# ### Dot product
# 
# The dot (or scalar) product $u\cdot v$ of the vector fields $u$ and $v$ is obtained by the method `dot_product`, which admits `dot` as a shortcut alias:

# In[21]:


u.dot(v) == u[1]*v[1] + u[2]*v[2] + u[3]*v[3]


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

# In[22]:


s = u.dot(v)
s


# In[23]:


print(s)


# In[24]:


s.display()


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

# In[25]:


s(p)


# Its coordinate expression is

# In[26]:


s.expr()


# ### Norm
# 
# The norm $\|u\|$ of the vector field $u$ is defined in terms
# of the dot product by $\|u\| = \sqrt{u\cdot u}$:

# In[27]:


norm(u) == sqrt(u.dot(u))


# It is a scalar field on $\mathbb{E}^3$:

# In[28]:


s = norm(u)
print(s)


# In[29]:


s.display()


# In[30]:


s.expr()


# For $v$, we have

# In[31]:


norm(v).expr()


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

# In[32]:


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


# In[33]:


s.display()


# We can check the standard formulas expressing the cross product in terms of
# the components:

# In[34]:


all([s[1] == u[2]*v[3] - u[3]*v[2],
     s[2] == u[3]*v[1] - u[1]*v[3],
     s[3] == u[1]*v[2] - u[2]*v[1]])


# ### Scalar triple product
# 
# Let us introduce a third vector field, $w$ say. 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[35]:


w = E.vector_field(name='w')
w[1] = x*z
w[2] = y*z
w.display()


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

# In[36]:


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


# In[37]:


s.expr()


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

# In[38]:


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


# ## How to evaluate the standard 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[39]:


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 some unspecified function of $(x,y,z)$:

# In[40]:


F = E.scalar_field(function('f')(x,y,z), name='F')
F.display()


# The value of $F$ at a point:

# In[41]:


F(p)


# The gradient of $F$:

# In[42]:


print(grad(F))


# In[43]:


grad(F).display()


# In[44]:


norm(grad(F)).display()


# ### Divergence 
# 
# The divergence of the vector field $u$:

# In[45]:


s = div(u)
s.display()


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

# In[46]:


div(v).expr()


# In[47]:


div(w).expr()


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

# In[48]:


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


# ### Curl
# 
# The curl of the vector field $u$:

# In[49]:


s = curl(u)
print(s)


# In[50]:


s.display()


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

# In[51]:


rot = curl


# An alternative is

# In[52]:


from sage.manifolds.operators import curl as rot


# We have then

# In[53]:


rot(u).display()


# In[54]:


rot(u) == curl(u)


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

# In[55]:


curl(v).display()


# In[56]:


curl(w).display()


# The curl of a gradient is always zero:

# In[57]:


curl(grad(F)).display()


# The divergence of a curl is always zero:

# In[58]:


div(curl(u)).display()


# An identity valid for any scalar field $F$ and any vector field $u$ is
# $$ \mathrm{curl}(Fu) = \mathrm{grad}\, F\times u + F\,  \mathrm{curl}\, u, $$
# as we can check:

# In[59]:


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


# ### Laplacian

# The Laplacian $\Delta F$ of a scalar field $F$ is another scalar
# field:

# In[60]:


s = laplacian(F)
s.display()


# The following identity holds:
# $$\Delta F = \mathrm{div}\left(\mathrm{grad}\, F\right)$$
# as we can check:

# In[61]:


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


# The Laplacian $\Delta u$ of a vector field $u$ is another vector field:

# In[62]:


Du = laplacian(u)
print(Du)


# whose components are

# In[63]:


Du.display()


# In the Cartesian frame, the components of $\Delta u$ are nothing but the
# (scalar) Laplacians of the components of $u$, as we can check:

# In[64]:


e = E.cartesian_frame()
Du == sum(laplacian(u[[i]])*e[i] for i in E.irange())


# In the above formula, `u[[i]]` return the $i$-th component of $u$ as a scalar field, while `u[i]` would have returned the coordinate expression of this scalar field; besides, `e` is the Cartesian frame:

# In[65]:


e[:]


# For the vector fields $v$ and $w$, we have 

# In[66]:


laplacian(v).display()


# In[67]:


laplacian(w).display()


# We have

# In[68]:


curl(curl(u)).display()


# In[69]:


grad(div(u)).display()


# A famous identity is
# 
# $$ \mathrm{curl}\left(\mathrm{curl}\, u\right) =
# \mathrm{grad}\left(\mathrm{div}\, u\right) - \Delta u.$$
# Let us check it:

# In[70]:


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


# ## How to customize various symbols?
# 
# ### Customizing the symbols of the orthonormal frame vectors
# 
# By default, the vectors of the orthonormal frame associated with Cartesian coordinates are denoted $(e_x,e_y,e_z)$:

# In[71]:


frame = E.cartesian_frame()
frame


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

# In[72]:


frame.set_name('a', indices=('x', 'y', 'z'))
frame


# In[73]:


v.display()


# In[74]:


frame.set_name(('hx', 'hy', 'hz'), 
               latex_symbol=(r'\mathbf{\hat{x}}', r'\mathbf{\hat{y}}', 
                             r'\mathbf{\hat{z}}'))
frame


# In[75]:


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[76]:


E.<x,y,z> = EuclideanSpace()


# which resulted in the coordinate symbols $(x,y,z)$ and in the corresponding Python variables `x`, `y` and `z` (SageMath symbolic expressions). To use other symbols, for instance $(X,Y,Z)$, it suffices to create `E` as

# In[77]:


E.<X,Y,Z> = EuclideanSpace()


# We have then:

# In[78]:


E.atlas()


# In[79]:


E.cartesian_frame()


# In[80]:


v = E.vector_field(-Y, X, sin(X*Y*Z), name='v')
v.display()


# By default the LaTeX symbols of the coordinate coincide with the letters given within the angle brackets. But this can be adjusted through the optional argument `symbols` of the function `EuclideanSpace`, which 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[81]:


E.<xi,et,ze> = EuclideanSpace(symbols=r"xi:\xi et:\eta ze:\zeta")
E.atlas()


# In[82]:


v = E.vector_field(-et, xi, sin(xi*et*ze), name='v')
v.display()