# 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.

# 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.

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)', [np.float], squeeze = 0)

kamodo = Kamodo()

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


## 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

@kamodofy(equation = "$$\\hat{n}(\\vec{y}) = \\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

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


## Solving the initial value problem¶

Generate a set of seed points for integration

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

@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

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


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.

kamodo.svec().head()

0 1 2
seed integral
0.0 -6.610169 -0.347554 -0.347554 -0.924068
-6.101695 -0.615865 -0.615865 -1.261579
-5.593220 -0.922720 -0.922720 -1.525495
-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:

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

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


compare with original seeds:

seeds[[0,1,2]]

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.

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

0 1 2
seed integral
0.0 -6.000 -0.674467 -0.674467 -1.320541
-5.000 -1.314574 -1.314574 -1.737228
0.000 -3.141593 -3.141593 1.000000
5.000 -0.120683 -0.120683 0.491788
6.000 0.115024 0.115024 -0.404619
4.0 4.000 0.094188 -0.156979 0.481816
-5.777 0.234822 -0.391371 -0.826931

## Plotting Fieldlines¶

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

import plotly.io as pio

fig = kamodo.plot('svec')
pio.write_image(fig, 'images/fieldlines.svg')


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

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


## Integration totals¶

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

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

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

totals[:5]

[10.677966101694915,
8.64406779661017,
7.627118644067796,
7.627118644067796,
8.64406779661017]


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

pd.DataFrame.groupby.a

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

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