http://nptel.ac.in/courses/115101005/downloads/lectures-doc/Lecture-25.pdf and Bhooshan https://www.amazon.com/Fundamentals-Engineering-Electromagnetics-Sunil-Bhooshan/dp/0198077947/ref=sr_1_4?ie=UTF8&qid=1519255078&sr=8-4&keywords=bhooshan page 268 ex 8.17
Also thanks to Stephan Hoyer @shoyer for help with https://github.com/pydata/xarray/issues/1773
from sympy import *
init_printing()
import numpy as np
import xarray as xr
import itertools as Iter
import yt
yt.toggle_interactivity()
from scipy import constants as consts
import matplotlib.pyplot as plt
%matplotlib notebook
Using matplotlib backend: Qt5Agg Warning: Cannot change to a different GUI toolkit: notebook. Using qt instead.
I, x, y, z, r, rho0, theta=symbols('I, x, y, z, r, rho_0, theta')
rdef=sqrt(x**2+y**2+z**2)
thetaDef=atan2(z, rho0)
rdef, thetaDef
Ay=(consts.mu_0*rho0**2*I/4)*sin(theta)/r**2
Ay
Ay=simplify(Ay.subs({theta:thetaDef, r: rdef})); Ay
sups={rho0:.5, I:1} #current loop with radidis .5m with 1A
sups
AyN=lambdify((x, y, z), Ay.subs(sups), dummify=False)
AyN(1,1,1)
Bx=simplify(-Derivative(Ay, z).doit()); Bx
BxN=lambdify((x, y, z), Bx.subs(sups), dummify=False)
BxN(1,1,1)
Bz=simplify(Derivative(Ay, x).doit()); Bz
BzN=lambdify((x, y, z), Bz.subs(sups), dummify=False)
BzN(1,1,1)
DomainSpaceSize = 6 # using cartesian 3D
SpaceDensity = [10] # 100 divisions in space 5 in time
x_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])
y_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])
z_coord = np.linspace(-DomainSpaceSize, +DomainSpaceSize, SpaceDensity[0])
CoorEnterFunc = lambda x, y, z: 1+0*x+0*y+0*z
dx=xr.DataArray(x_coord, dims='x')
dy=xr.DataArray(y_coord, dims='y')
dz=xr.DataArray(z_coord, dims='z')
SimDataSet = xr.Dataset({'Ay':(['x', 'y', 'z'], AyN(*np.meshgrid(x_coord, y_coord, z_coord)))},
coords={'x':x_coord, 'y':y_coord, 'z':z_coord})
SimDataSet.attrs['units']={'Ay':'V/s/m'}
SimDataSet
<xarray.Dataset> Dimensions: (x: 10, y: 10, z: 10) Coordinates: * x (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * y (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * z (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... Data variables: Ay (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ... Attributes: units: {'Ay': 'V/s/m'}
SimDataSet['Bx']=(['x', 'y', 'z',],BxN(*np.meshgrid(x_coord, y_coord, z_coord)))
SimDataSet.attrs['units']['Bx']='T'
SimDataSet
<xarray.Dataset> Dimensions: (x: 10, y: 10, z: 10) Coordinates: * x (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * y (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * z (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... Data variables: Ay (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ... Bx (x, y, z) float64 7.969e-11 8.085e-11 6.879e-11 2.328e-11 ... Attributes: units: {'Ay': 'V/s/m', 'Bx': 'T'}
SimDataSet['Bz']=(['x', 'y', 'z',],BzN(*np.meshgrid(x_coord, y_coord, z_coord)))
SimDataSet.attrs['units']['Bz']='T'
SimDataSet
<xarray.Dataset> Dimensions: (x: 10, y: 10, z: 10) Coordinates: * x (x) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * y (y) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... * z (z) float64 -6.0 -4.667 -3.333 -2.0 -0.6667 0.6667 2.0 3.333 ... Data variables: Ay (x, y, z) float64 -7.247e-10 -8.327e-10 -9.345e-10 -1.003e-09 ... Bx (x, y, z) float64 7.969e-11 8.085e-11 6.879e-11 2.328e-11 ... Bz (x, y, z) float64 -8.052e-11 -1.066e-10 -1.349e-10 -1.583e-10 ... Attributes: units: {'Ay': 'V/s/m', 'Bx': 'T', 'Bz': 'T'}
bbox=[]
for i in ['x', 'y', 'z']:
bbox.append([float(SimDataSet.coords[i].min()), float(SimDataSet.coords[i].max())])
bbox=np.array(bbox)
bbox
array([[-6., 6.], [-6., 6.], [-6., 6.]])
YTDataObj= dict(MagPoty = (np.array(SimDataSet['Ay']), 'g/m**3'),
MagFieldx = (np.array(SimDataSet['Bx']), 'g/m**3'),
MagFieldz = (np.array(SimDataSet['Bz']), 'g/m**3')
)
#would not except unit in SimDataSet.attrs['units']['ChargeDen'] had to fake with
#'g/m**3'
YTDataObj = yt.load_uniform_grid(YTDataObj, tuple(dict(SimDataSet.dims).values()), length_unit="m", bbox=bbox, nprocs=64)
yt : [INFO ] 2018-02-21 16:29:44,209 Parameters: current_time = 0.0 yt : [INFO ] 2018-02-21 16:29:44,211 Parameters: domain_dimensions = [10 10 10] yt : [INFO ] 2018-02-21 16:29:44,212 Parameters: domain_left_edge = [-6. -6. -6.] yt : [INFO ] 2018-02-21 16:29:44,212 Parameters: domain_right_edge = [ 6. 6. 6.] yt : [INFO ] 2018-02-21 16:29:44,214 Parameters: cosmological_simulation = 0.0
Trying to do transfer function but giving up
tf = yt.ColorTransferFunction((SimDataSet['Ay'].min(), SimDataSet['Ay'].max()), grey_opacity=False)
tf.map_to_colormap(SimDataSet['Ay'].min(), SimDataSet['Ay'].max(), scale=15.0, colormap='algae')
MagPotyScene = yt.create_scene(YTDataObj, 'MagPoty')
source2 = MagPotyScene[0]
MagPotyScene.camera.set_width(MagPotyScene.quan(7, 'm'))
MagPotyScene.show(sigma_clip=3)
# the transfer function should have a defualt to the max min of the data cube
#being renderd
yt : [INFO ] 2018-02-21 16:29:45,197 Rendering scene (Can take a while). yt : [INFO ] 2018-02-21 16:29:45,200 Creating volume /home/iridium/anaconda3/lib/python3.6/site-packages/yt/units/yt_array.py:1293: RuntimeWarning: invalid value encountered in log10 out_arr = func(np.asarray(inp), out=out, **kwargs) yt : [INFO ] 2018-02-21 16:29:45,315 Creating transfer function yt : [INFO ] 2018-02-21 16:29:45,316 Calculating data bounds. This may take a while. Set the TranferFunctionHelper.bounds to avoid this.
from yt.visualization.volume_rendering.api import Scene, VolumeSource
sc = Scene()
cam = sc.add_camera(YTDataObj, lens_type='perspective')
cam.resolution = [400, 400]
cam.position = YTDataObj.arr([9, 9, 9], 'm')
cam.switch_orientation()
BFx = VolumeSource(YTDataObj, field='MagFieldx')
BFx.use_ghost_zones = True
sc.add_source(BFx)
sc.show( sigma_clip=1)
yt : [INFO ] 2018-02-21 16:29:49,025 Rendering scene (Can take a while). yt : [INFO ] 2018-02-21 16:29:49,027 Creating volume /home/iridium/anaconda3/lib/python3.6/site-packages/yt/units/yt_array.py:1293: RuntimeWarning: invalid value encountered in log10 out_arr = func(np.asarray(inp), out=out, **kwargs) yt : [INFO ] 2018-02-21 16:29:49,692 Creating transfer function yt : [INFO ] 2018-02-21 16:29:49,693 Calculating data bounds. This may take a while. Set the TranferFunctionHelper.bounds to avoid this.
# add rendering of x-velocity field
BFz = VolumeSource(YTDataObj, field='MagFieldz')
BFz.use_ghost_zones = True
sc.add_source(BFz)
sc.show( sigma_clip=1)
Not sure if this consulted the Bx and Bz fields into one image