Bugs and improvement: yangdikun@gmail.com
In [1]:
import numpy as np
from ipywidgets import interact, interactive, fixed, Layout
import ipywidgets as widgets
In [2]:
import plotly.plotly as py
import plotly.graph_objs as go
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
print(__version__) # requires version >= 1.9.0
init_notebook_mode(connected=True)
In [3]:
# x, y, z: location of the observing point
# Bx, By, Bz: vector B field at the observing point in vaccum
def VerticalMagneticDipoleField(x, y, z):
r = (x**2 + y**2 + z**2)**(0.5)
Bx = 1e-7/r**3 * 3*x*z/r**2
By = 1e-7/r**3 * 3*y*z/r**2
Bz = 1e-7/r**3 * (3*z**2/r**2 - 1)
return Bx, By, Bz
In [4]:
# radius: a float defining how far the field lines expand
# stepSize: length of the segments that make the field lines
# yloc, zloc: a single field line path in the y-z plane
def VerticalMagneticDipoleLine(radius,stepSize=1.):
x, y, z = 0., [radius], [-stepSize/4.]
while (y[-1]>=0.) & (z[-1] <= 0.):
_, By, Bz = VerticalMagneticDipoleField(x, y[-1], z[-1])
B = (By**2 + Bz**2)**(0.5)
By, Bz = By/B*stepSize, Bz/B*stepSize
y.append(y[-1]+By)
z.append(z[-1]+Bz)
y.pop()
z.pop()
yloc = np.r_[y[-1:0:-1], y]
zloc = np.r_[-np.array(z[-1:0:-1]),z]
return yloc, zloc
In [5]:
# dipoleLoc: (x,y,z) of a magnetic dipole
# dipoleDec, dipoleInc: declination and inclination angles of the dipole direction in degree
# radii: a float vector defining the radii of the field line shells
# naz: number of azimuths for the plotting of field lines
# stepSize: length of the segments that make the field lines
# xloc, yloc, zloc: a collection of all field line paths
def MagneticDipoleLine(dipoleLoc,dipoleDec,dipoleInc,radii,naz=10,stepSize=1.):
x0, y0, z0 = dipoleLoc[0], dipoleLoc[1], dipoleLoc[2]
# rotation matrix
theta, alpha = -np.pi*(dipoleInc+90.)/180., -np.pi*dipoleDec/180.
Rx = np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])
Rz = np.array([[np.cos(alpha),-np.sin(alpha),0.],[np.sin(alpha),np.cos(alpha),0.],[0.,0.,1.]])
R = np.dot(Rz,Rx) # Rz @ Rx
azimuth = np.linspace(0.,2*np.pi,num=naz,endpoint=False)
xloc, yloc, zloc = [], [], []
for r in radii:
hloc, vloc = VerticalMagneticDipoleLine(r,stepSize)
for a in azimuth:
x, y, z = np.sin(a)*hloc, np.cos(a)*hloc, vloc
xyz = np.dot(R, np.vstack((x,y,z)))
xloc.append(xyz[0]+x0)
yloc.append(xyz[1]+y0)
zloc.append(xyz[2]+z0)
return xloc, yloc, zloc
In [6]:
# dipoleLoc: (x,y,z) of a magnetic dipole
# dipoleDec, dipoleInc: declination and inclination angles of the dipole direction in degree
# radius: radius of loop surface
# nseg: number of line segments that make a loop
# xloc, yloc, zloc: a collection of line paths for a loop and its normal direction
def MagneticDipoleLoop(dipoleLoc,dipoleDec,dipoleInc,radius,nseg=20):
x0, y0, z0 = dipoleLoc[0], dipoleLoc[1], dipoleLoc[2]
# rotation matrix
theta, alpha = -np.pi*(dipoleInc+90.)/180., -np.pi*dipoleDec/180.
Rx = np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])
Rz = np.array([[np.cos(alpha),-np.sin(alpha),0.],[np.sin(alpha),np.cos(alpha),0.],[0.,0.,1.]])
R = np.dot(Rz,Rx) # Rz @ Rx
# loop path
azimuth = np.linspace(0.,2*np.pi,num=nseg+1,endpoint=True)
x, y, z = np.sin(azimuth)*radius, np.cos(azimuth)*radius, azimuth*0.
xyz = np.dot(R, np.vstack((x,y,z)))
xloc1, yloc1, zloc1 = xyz[0]+x0, xyz[1]+y0, xyz[2]+z0
# normal direction pointer
x, y, z = [0.,0.],[0.,0.],[0.,radius]
xyz = np.dot(R, np.vstack((x,y,z)))
xloc2, yloc2, zloc2 = xyz[0]+x0, xyz[1]+y0, xyz[2]+z0
return [xloc1,xloc2], [yloc1,yloc2], [zloc1,zloc2]
In [7]:
def plot_magDipoleEM(x1,y1,z1,dec1,inc1,nazimuth,nshell,rmax,x2,y2,z2,dec2,inc2,loop1flag,line1flag,loop2flag):
# internal parameters
loopRadius1, loopRadius2 = 0.5, 0.5
stepSize = 0.2 # length of field line segment
# get line path for the 1st loop
loopx1, loopy1, loopz1 = MagneticDipoleLoop((x1,y1,z1),dec1,inc1,loopRadius1)
# get line path for the field lines of the 1st loop
radii = np.linspace(rmax,0.,nshell,endpoint=False)
linex1, liney1, linez1 = MagneticDipoleLine((x1,y1,z1),dec1,inc1,radii,nazimuth,stepSize)
# get line path for the 2nd loop
loopx2, loopy2, loopz2 = MagneticDipoleLoop((x2,y2,z2),dec2,inc2,loopRadius2)
# create plot
data = []
if loop1flag:
for xx,yy,zz in zip(loopx1,loopy1,loopz1):
data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='red',width=3)))
if line1flag:
for xx,yy,zz in zip(linex1,liney1,linez1):
data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='blue',width=3)))
if loop2flag:
for xx,yy,zz in zip(loopx2,loopy2,loopz2):
data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='green',width=3)))
layout = dict(width=1000, height=500, autosize=False, showlegend=False,
margin=go.Margin(l=5,r=5,b=5,t=5,pad=4),paper_bgcolor='#efefef',
scene = dict(aspectmode = 'data', aspectratio = dict(x=1, y=1, z=1)))
fig = dict(data=data, layout=layout)
iplot(fig)
return
interact(plot_magDipoleEM,
x1=fixed(0.), y1=fixed(0.), z1=fixed(0.),
dec1=widgets.IntSlider(min=0,max=360,step=1,value=0,description='Loop 1 dec'),
inc1=widgets.IntSlider(min=-90,max=90,step=1,value=-90,description='Loop 1 inc'),
nazimuth=fixed(8),
nshell=widgets.IntSlider(min=1,max=5,step=1,value=2,description='Number of line layers',style = {'description_width': 'initial'}),
rmax=widgets.FloatSlider(min=1.,max=20.,step=1.,value=5.,description='Max. radius of lines',style = {'description_width': 'initial'}),
x2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=10.,description='Loop 2 X'),
y2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=0.,description='Loop 2 Y'),
z2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=0.,description='Loop 2 Z'),
dec2=widgets.IntSlider(min=0,max=360,step=1,value=0,description='Loop 2 dec'),
inc2=widgets.IntSlider(min=-90,max=90,step=1,value=-90,description='Loop 2 inc'),
loop1flag=widgets.Checkbox(value=True,disabled=False,description='Show Loop 1'),
line1flag=widgets.Checkbox(value=False,disabled=False,description='Show field lines of Loop 1'),
loop2flag=widgets.Checkbox(value=False,disabled=False,description='Show Loop 2')
)
Out[7]:
In [ ]: