Field Integration Techniques¶

Many analysis techniques for vector fields require solving an initial value problem for an arbitrary set of seed points and evaluating such solutions at a chosen resolution. Kamodo makes it easy to generate fieldline solutions by providing a function decorator that wraps scipy's powerful solve_ivp function. Each family of solutions is represented by a single function of a complex parameter. We illustrate the flexibility of this approach in the example below.

In [1]:

Copied!





# initialize
from plotly.offline import iplot, plot, init_notebook_mode
init_notebook_mode(connected = True)

from kamodo import Kamodo, event, pointlike, kamodofy, solve
import numpy as np
import pandas as pd
# initialize
from plotly.offline import iplot, plot, init_notebook_mode
init_notebook_mode(connected = True)

from kamodo import Kamodo, event, pointlike, kamodofy, solve
import numpy as np
import pandas as pd

Dipole field model¶

We use the following dipole field model that can accept (m,) and (1,m), and (n,m) arrays.

In [2]:

Copied!





def Bdip(rvec):
    """Need math to work in a variety of arg shapes"""
    muvec = Bdip.muvec    
    r = np.linalg.norm(rvec, axis = 1)
    r[r==0] = np.nan

    try:
        rhat = rvec/r
    except:
        rhat = (rvec.T/r).T

    try:
        result = 3*np.dot(rhat, muvec.T)
    except:
        result = 3*np.dot(rhat.T, muvec.T).T


    result = (rhat.T*result).T

    try:
        result = result - muvec
    except:
        result = (result - muvec.T).T

    try:
        result = result/r**3
    except:
        result = (result.T/r**3).T

    return result

# set dipole moment
Bdip.muvec = np.array([0, 0, -1]) 

# pointlike enforces dimensionality
Bdip = pointlike(Bdip, '(n,m)->(n,m)', [float], squeeze = 0)
def Bdip(rvec):
    """Need math to work in a variety of arg shapes"""
    muvec = Bdip.muvec    
    r = np.linalg.norm(rvec, axis = 1)
    r[r==0] = np.nan

    try:
        rhat = rvec/r
    except:
        rhat = (rvec.T/r).T

    try:
        result = 3*np.dot(rhat, muvec.T)
    except:
        result = 3*np.dot(rhat.T, muvec.T).T


    result = (rhat.T*result).T

    try:
        result = result - muvec
    except:
        result = (result - muvec.T).T

    try:
        result = result/r**3
    except:
        result = (result.T/r**3).T

    return result

# set dipole moment
Bdip.muvec = np.array([0, 0, -1]) 

# pointlike enforces dimensionality
Bdip = pointlike(Bdip, '(n,m)->(n,m)', [float], squeeze = 0)

In [3]:

Copied!

kamodo = Kamodo()
kamodo = Kamodo()

In [4]:

Copied!

kamodo['Bvec'] = Bdip # register the dipole field
kamodo
kamodo['Bvec'] = Bdip # register the dipole field
kamodo

Out[4]:

\begin{equation}\vec{B}{\left(\vec{r} \right)} = \lambda{\left(\vec{r} \right)}\end{equation}

Bvec works on a list of points and on individual points

In [5]:

Copied!





kamodo.Bvec([[1,0,0], # x,y,z
             [2,0,0],
             [0,0,1],
             [0,0,2],])
kamodo.Bvec([[1,0,0], # x,y,z
             [2,0,0],
             [0,0,1],
             [0,0,2],])

Out[5]:

array([[ 0.   ,  0.   ,  1.   ],
       [ 0.   ,  0.   ,  0.125],
       [-0.   , -0.   , -2.   ],
       [-0.   , -0.   , -0.25 ]])

In [6]:

Copied!

kamodo.Bvec([1,0,0])
kamodo.Bvec([1,0,0])

Out[6]:

array([0., 0., 1.])

Normalization¶

Instead of solving the initial value problem on the original field, we will be solving on the normalized field. This will mean that the integral path is the same as the arclength, allowing us to control the visual fidelity of the resulting field.

Create a normalization function to be applied to our field

In [7]:

Copied!





@kamodofy(equation = "\\vec{y}/\\sqrt{\\vec{y} \\cdot \\vec{y}}")
@pointlike(signature = '(m,n)->(m,n)', squeeze = 0)
def normalized(yvec):   
    r = np.linalg.norm(yvec, axis = 1)
    r[r==0] = np.nan

    try:
        return yvec/r
    except:
        return (yvec.T/r).T


kamodo['nhat'] = normalized
@kamodofy(equation = "\\vec{y}/\\sqrt{\\vec{y} \\cdot \\vec{y}}")
@pointlike(signature = '(m,n)->(m,n)', squeeze = 0)
def normalized(yvec):   
    r = np.linalg.norm(yvec, axis = 1)
    r[r==0] = np.nan

    try:
        return yvec/r
    except:
        return (yvec.T/r).T


kamodo['nhat'] = normalized

Create a normalized field

In [8]:

Copied!

kamodo['bhat'] = "nhat(Bvec)"
kamodo
kamodo['bhat'] = "nhat(Bvec)"
kamodo

Out[8]:

\begin{equation}\vec{B}{\left(\vec{r} \right)} = \lambda{\left(\vec{r} \right)}\end{equation} \begin{equation}\hat{n}{\left(\vec{y} \right)} = \vec{y}/\sqrt{\vec{y} \cdot \vec{y}}\end{equation} \begin{equation}\hat{b}{\left(\vec{r} \right)} = \hat{n}{\left(\vec{B}{\left(\vec{r} \right)} \right)}\end{equation}

In [9]:

Copied!





kamodo.bhat([[1,0,0], # x,y,z
             [2,0,0],
             [0,0,1],
             [0,0,2]])
kamodo.bhat([[1,0,0], # x,y,z
             [2,0,0],
             [0,0,1],
             [0,0,2]])

Out[9]:

array([[ 0.,  0.,  1.],
       [ 0.,  0.,  1.],
       [-0., -0., -1.],
       [-0., -0., -1.]])

Solving the initial value problem¶

Generate a set of seed points for integration

In [10]:

Copied!

x0 = np.linspace(-np.pi,np.pi,6)
y0 = np.linspace(-np.pi,np.pi,6)
z0 = 1

seeds = np.array(np.column_stack([c.ravel() for c in np.meshgrid(x0,y0,z0)]))
x0 = np.linspace(-np.pi,np.pi,6)
y0 = np.linspace(-np.pi,np.pi,6)
z0 = 1

seeds = np.array(np.column_stack([c.ravel() for c in np.meshgrid(x0,y0,z0)]))

Create a stopping boundary for field line integrator

In [11]:

Copied!





@event
def boundary(s, rvec):
    r = np.linalg.norm(rvec)
    
    if np.isnan(r):
        result = 0
    else:
        result = r - 1
    return result
@event
def boundary(s, rvec):
    r = np.linalg.norm(rvec)
    
    if np.isnan(r):
        result = 0
    else:
        result = r - 1
    return result

Solve the initial value problem for the normalized field

In [12]:

Copied!





kamodo['svec'] = solve(kamodo.bhat, # the field to be solved
                       seeds, # the initial positions
                       's', # the name of the integration parameter
                       (0,30), # the span to integrate over
                       npoints = 60, # the number of points to evaluate the solution
                       events = boundary, # stop at the boundary
                      )
kamodo
kamodo['svec'] = solve(kamodo.bhat, # the field to be solved
                       seeds, # the initial positions
                       's', # the name of the integration parameter
                       (0,30), # the span to integrate over
                       npoints = 60, # the number of points to evaluate the solution
                       events = boundary, # stop at the boundary
                      )
kamodo

Out[12]:

\begin{equation}\vec{B}{\left(\vec{r} \right)} = \lambda{\left(\vec{r} \right)}\end{equation} \begin{equation}\hat{n}{\left(\vec{y} \right)} = \vec{y}/\sqrt{\vec{y} \cdot \vec{y}}\end{equation} \begin{equation}\hat{b}{\left(\vec{r} \right)} = \hat{n}{\left(\vec{B}{\left(\vec{r} \right)} \right)}\end{equation} \begin{equation}\vec{s}{\left(s \right)} = \lambda{\left(s \right)}\end{equation}

The solver returns a family of solutions, represented as a single function of a complex array, $\vec{s}(s)$ where $s$ is a complex array.

Evaluating the Solutions¶

On evaluation, $\vec{s}(s)$ returns a pandas dataframe.

In [13]:

Copied!

kamodo.svec().head()
kamodo.svec().head()

Out[13]:

		0	1	2
seed	integral
0.0	-6.610169	-0.347547	-0.347547	-0.924155
	-6.101695	-0.615886	-0.615886	-1.261575
	-5.593220	-0.922735	-0.922735	-1.525472
	-5.084746	-1.256963	-1.256963	-1.713145
	-4.576271	-1.608411	-1.608411	-1.822192

When using the default argument above, the solution evaluates at a resolution of npoints/span, stopping at the boundary.

Complex parameterization¶

Kamodo represents the family of solutions to the initial value problem as a single function of a complex array.

The floor of the real part of the input parameter corresponds to the original seed array:

In [14]:

Copied!

kamodo.svec([0,1,2]).values
kamodo.svec([0,1,2]).values

Out[14]:

array([[-3.14159265, -3.14159265,  1.        ],
       [-1.88495559, -3.14159265,  1.        ],
       [-0.62831853, -3.14159265,  1.        ]])

compare with original seeds:

In [15]:

Copied!

seeds[[0,1,2]]
seeds[[0,1,2]]

Out[15]:

array([[-3.14159265, -3.14159265,  1.        ],
       [-1.88495559, -3.14159265,  1.        ],
       [-0.62831853, -3.14159265,  1.        ]])

The imaginary part denotes the integral along the corresponding solution. Here, we can choose evaluation points that were not in the original solution. Parameters outside the original span will be extrapolated.

In [16]:

Copied!

kamodo.svec([-6j, -5j, 0, 5j, 6j, 4 + 4j, 4 -5.777j])
kamodo.svec([-6j, -5j, 0, 5j, 6j, 4 + 4j, 4 -5.777j])

Out[16]:

		0	1	2
seed	integral
0.0	-6.000	-0.674502	-0.674502	-1.320504
	-5.000	-1.314574	-1.314574	-1.737228
	0.000	-3.141593	-3.141593	1.000000
	5.000	-0.120606	-0.120606	0.491892
	6.000	0.125472	0.125472	-0.393292
4.0	4.000	0.094223	-0.157038	0.481740
4.0	-5.777	0.234804	-0.391340	-0.827054

Plotting Fieldlines¶

We can quickly generate plots for all fieldlines at the default resolution by calling plot with the name of the fieldlines solution.

In [17]:

Copied!

import plotly.io as pio
import plotly.io as pio

In [18]:

Copied!

#fig = kamodo.plot('svec')
#fig = kamodo.plot('svec')

In [19]:

Copied!

#pio.write_image(fig, './images/fieldlines.svg')
#pio.write_image(fig, './images/fieldlines.svg')

images/fieldlines.svg

To show the direction of the field at each point, we can evaluate $\hat{B}(\vec{s}(s))$

In [20]:

Copied!

#fig = kamodo.plot('svec', 
#                  bhat = dict(rvec = kamodo.svec()))
#pio.write_image(fig,'./images/fieldlines_vectors.svg')
#fig = kamodo.plot('svec', 
#                  bhat = dict(rvec = kamodo.svec()))
#pio.write_image(fig,'./images/fieldlines_vectors.svg')

fieldlines

Integration totals¶

To compute the total integral for each fieldline individually, we need a function to subtract the integration results at the endpoints.

In [21]:

Copied!

def integral(fieldline):
    endpoints = fieldline.reset_index().integral.iloc[[0,-1]]
    return endpoints.values[-1] - endpoints.values[0]
def integral(fieldline):
    endpoints = fieldline.reset_index().integral.iloc[[0,-1]]
    return endpoints.values[-1] - endpoints.values[0]

In [22]:

Copied!

totals = []
for seed, fieldline in kamodo.svec().groupby(level = 'seed'):
    totals.append(integral(fieldline))
    
totals[:5]
totals = []
for seed, fieldline in kamodo.svec().groupby(level = 'seed'):
    totals.append(integral(fieldline))
    
totals[:5]

Out[22]:

[10.677966101694915,
 8.64406779661017,
 7.627118644067796,
 7.627118644067796,
 8.64406779661017]

Alternatively, we can use pandas' aggregation methods to apply our function on each fieldline.

In [23]:

Copied!

kamodo.svec().groupby(level='seed').aggregate(integral)
kamodo.svec().groupby(level='seed').aggregate(integral)

Out[23]:

	0	1	2
seed
0.0	10.677966	10.677966	10.677966
1.0	8.644068	8.644068	8.644068
2.0	7.627119	7.627119	7.627119
3.0	7.627119	7.627119	7.627119
4.0	8.644068	8.644068	8.644068
5.0	10.677966	10.677966	10.677966
6.0	8.644068	8.644068	8.644068
7.0	6.610169	6.610169	6.610169
8.0	5.084746	5.084746	5.084746
9.0	5.084746	5.084746	5.084746
10.0	6.610169	6.610169	6.610169
11.0	8.644068	8.644068	8.644068
12.0	7.627119	7.627119	7.627119
13.0	5.084746	5.084746	5.084746
14.0	5.593220	5.593220	5.593220
15.0	5.593220	5.593220	5.593220
16.0	5.084746	5.084746	5.084746
17.0	7.627119	7.627119	7.627119
18.0	7.627119	7.627119	7.627119
19.0	5.084746	5.084746	5.084746
20.0	5.593220	5.593220	5.593220
21.0	5.593220	5.593220	5.593220
22.0	5.084746	5.084746	5.084746
23.0	7.627119	7.627119	7.627119
24.0	8.644068	8.644068	8.644068
25.0	6.610169	6.610169	6.610169
26.0	5.084746	5.084746	5.084746
27.0	5.084746	5.084746	5.084746
28.0	6.610169	6.610169	6.610169
29.0	8.644068	8.644068	8.644068
30.0	10.677966	10.677966	10.677966
31.0	8.644068	8.644068	8.644068
32.0	7.627119	7.627119	7.627119
33.0	7.627119	7.627119	7.627119
34.0	8.644068	8.644068	8.644068
35.0	10.677966	10.677966	10.677966