Field Integration Techniques¶

# 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

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)

kamodo = Kamodo()

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

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

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

@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

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

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

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)]))

@event
def boundary(s, rvec):
    r = np.linalg.norm(rvec)

    if np.isnan(r):
        result = 0
    else:
        result = r - 1
    return result

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().head()

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

seeds[[0,1,2]]

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

import plotly.io as pio

fig = kamodo.plot('svec')

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

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

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]

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