Loading pysen/physics/biotsavart.py +96 −3 Original line number Diff line number Diff line Loading @@ -14,14 +14,15 @@ u0 = 4*pi * 10^-7 """ import numpy as np from numpy import (pi, radians, asarray, ones, zeros_like, linspace, sqrt, sin, cos, arctan2, inner, diag, cross) from numpy import (pi, radians, asarray, ones, zeros_like, linspace, sign, log, sqrt, sin, cos, arctan2, inner, diag, cross, arctan, dot) #pylint: disable=no-name-in-module from scipy.special import ellipk, ellipe #pylint: enable=no-name-in-module _four_pi = 4*pi _EPSILON = 1e-6 # ============================================================================ # SPECIAL CASES Loading Loading @@ -200,6 +201,98 @@ def solenoid(r, radius=1, start=-0.5, end=+0.5, nturns=10): bfield = bfield + db return bfield # ============================================================================ # CURRENT SHEETS, FLIPPERS # ============================================================================ def bsheet_kernel(r, width, height): """ Core of the field calculation: * current sheet is centered at the origin of the coordinate system. * current sheet lies in the x-z plane (width in x-direction, height in z-direction) * current flows in the direction of the z-axis. * observation point is at (x,y,z) with y!=0 ^ (z) Coordinate system for the current sheet in the x-z plane | |------|------+ | ^ ^ ^ ^ ^ | | | | | | | | | | | | | | | ---------+----------> (x) | | | | | | | | | | | | | | | | | | | | | +------|------+ The result must be mutliplied by mu0 and the current density (A/m) """ bx, by = 0, 0 r = r.T x = r[0] y = r[1] z = r[2] y = np.where(abs(y)<_EPSILON, _EPSILON, y) # avoid singularity in the plane of the sheet for dx in (-width/2, +width/2): for dz in (-height/2, +height/2): k = sign(dx)*sign(dz) xx = x + dx zz = z + dz rr = sqrt(xx*xx + zz*zz + y*y) bx = bx - k*arctan(xx*zz/(y*rr)) by = by + k*log((rr-zz)/(rr+zz))/2 return bx/_four_pi, by/_four_pi def current_sheet_finite(r, s1, s2, s3): """ Compute the magnetic field of a rectangular current sheet with a current density of 1 A/m. The reference point is given at r(1..3) s1|------------------|s2 | -------> | s3|------------------| A rectangular current sheet in an arbitrary orientation is defined by the three corners: s1, s2, s3. The current flows from s1 to s2. Translated from bfelder.c:bsheet subroutine by M. Monkenbusch """ ex = s3 - s1 ez = s2 - s1 bb = sqrt(sum(ex**2)) hh = sqrt(sum(ez**2)) ex = ex/bb ez = ez/hh if abs(dot(ex,ez))>_EPSILON: raise RuntimeError("current_sheet: sheet is not rectangular!") ey = cross( ez, ex ) M = np.vstack((ex,ey,ez)) r = r - (s2+s3)/2 r = np.dot(M,r) bx, by = bsheet_kernel(r,bb,hh) # Transformation of the field back into the original coordinate system return bx*ex + by*ey def current_sheet_infinite(r): """ Compute the magnetic field of an infinite current sheet with a current density of 1 A/m. The reference point is given at r(1..3) """ return None # FIXME: not yet implemented def current_sheet(r, s1, s2, s3, finite=True): """Calculate the magnetic field of a current sheet""" if finite: return current_sheet_finite(r, s1, s2, s3) return current_sheet_infinite(r) # ============================================================================ # OBSOLETE # ============================================================================ Loading test/test_physics_biotsavart.py +58 −4 Original line number Diff line number Diff line #!/usr/bin/env python """ Python unit test cases for mathutil Python unit test cases for Biot-Savart module """ import unittest import numpy as np from numpy import pi, sqrt, cos, allclose from numpy import pi, sqrt, cos, allclose, array from scipy.constants import mu_0 import pysen.physics.biotsavart as bs Loading @@ -18,7 +18,7 @@ L = 1.0 # coil separation in m class BiotSavartTestCase(unittest.TestCase): "Biot-Savart Test Cases" x1,x2 = np.array((-L, 0, 0)), np.array(( L, 0, 0)) # two point segment x1,x2 = array((-L, 0, 0)), array(( L, 0, 0)) # two point segment x = np.linspace(-2*L,2*L, 11) _r = np.ones_like (x)*R Loading Loading @@ -175,7 +175,61 @@ class BiotSavartTestCase(unittest.TestCase): act = bs.solenoid_axis(0, R=R, L=LS, N=N) self.assertAlmostEqual(exp, act, 5) def test_current_sheet_simple(self): "test simple cases for current sheet" xs1 = array((-1, 0,-1)) xs2 = array((-1, 0,+1)) xs3 = array((+1, 0,-1)) r = array((0.1,0.2,0.3)) bexp = array((-0.40581227, 0.02079833, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) xs1 = array((-100, 0,-100)) xs2 = array((-100, 0,+100)) xs3 = array((+100, 0,-100)) r = array((0.0,2e-6,0.0)) bexp = array((-0.5, 0, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) r = array((0.0,-2e-6,0.0)) bexp = array((0.5, 0, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) xs1 = array((-2, 0,-1)) xs2 = array((-2, 0,+1)) xs3 = array((+1, 0,-1)) r = array((0.4,0.6,0.5)) bexp = array((-0.23694913, 0.08906743, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) def test_current_sheet_rotated(self): "test rotated sheet" xs1 = array((-1, 0, 0)) xs2 = array(( 0, 0,+1)) xs3 = array(( 0, 0,-1)) r = array((0.1,0.2,0.3)) bexp = array((-0.255651, -0.0390186, 0.255651)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) def xtest_current_sheet_vectorized(self): "test vectorized current sheet" xs1 = array((-1, 0,-1)) xs2 = array((-1, 0,+1)) xs3 = array((+1, 0,-1)) r = np.array([[0.1,0.2,0.3], [0.0,2e-6,0.0], [0.0,-2e-6,0.0], [0.4,0.6,0.5]]) bexp = np.array([[-0.40581227, 0.02079833, 0], [-0.5 , 0 , 0], [ 0.5 , 0 , 0], [-0.23694913, 0.08906743, 0]]) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) if __name__ == "__main__": unittest.main(exit=False, verbosity=2) Loading
pysen/physics/biotsavart.py +96 −3 Original line number Diff line number Diff line Loading @@ -14,14 +14,15 @@ u0 = 4*pi * 10^-7 """ import numpy as np from numpy import (pi, radians, asarray, ones, zeros_like, linspace, sqrt, sin, cos, arctan2, inner, diag, cross) from numpy import (pi, radians, asarray, ones, zeros_like, linspace, sign, log, sqrt, sin, cos, arctan2, inner, diag, cross, arctan, dot) #pylint: disable=no-name-in-module from scipy.special import ellipk, ellipe #pylint: enable=no-name-in-module _four_pi = 4*pi _EPSILON = 1e-6 # ============================================================================ # SPECIAL CASES Loading Loading @@ -200,6 +201,98 @@ def solenoid(r, radius=1, start=-0.5, end=+0.5, nturns=10): bfield = bfield + db return bfield # ============================================================================ # CURRENT SHEETS, FLIPPERS # ============================================================================ def bsheet_kernel(r, width, height): """ Core of the field calculation: * current sheet is centered at the origin of the coordinate system. * current sheet lies in the x-z plane (width in x-direction, height in z-direction) * current flows in the direction of the z-axis. * observation point is at (x,y,z) with y!=0 ^ (z) Coordinate system for the current sheet in the x-z plane | |------|------+ | ^ ^ ^ ^ ^ | | | | | | | | | | | | | | | ---------+----------> (x) | | | | | | | | | | | | | | | | | | | | | +------|------+ The result must be mutliplied by mu0 and the current density (A/m) """ bx, by = 0, 0 r = r.T x = r[0] y = r[1] z = r[2] y = np.where(abs(y)<_EPSILON, _EPSILON, y) # avoid singularity in the plane of the sheet for dx in (-width/2, +width/2): for dz in (-height/2, +height/2): k = sign(dx)*sign(dz) xx = x + dx zz = z + dz rr = sqrt(xx*xx + zz*zz + y*y) bx = bx - k*arctan(xx*zz/(y*rr)) by = by + k*log((rr-zz)/(rr+zz))/2 return bx/_four_pi, by/_four_pi def current_sheet_finite(r, s1, s2, s3): """ Compute the magnetic field of a rectangular current sheet with a current density of 1 A/m. The reference point is given at r(1..3) s1|------------------|s2 | -------> | s3|------------------| A rectangular current sheet in an arbitrary orientation is defined by the three corners: s1, s2, s3. The current flows from s1 to s2. Translated from bfelder.c:bsheet subroutine by M. Monkenbusch """ ex = s3 - s1 ez = s2 - s1 bb = sqrt(sum(ex**2)) hh = sqrt(sum(ez**2)) ex = ex/bb ez = ez/hh if abs(dot(ex,ez))>_EPSILON: raise RuntimeError("current_sheet: sheet is not rectangular!") ey = cross( ez, ex ) M = np.vstack((ex,ey,ez)) r = r - (s2+s3)/2 r = np.dot(M,r) bx, by = bsheet_kernel(r,bb,hh) # Transformation of the field back into the original coordinate system return bx*ex + by*ey def current_sheet_infinite(r): """ Compute the magnetic field of an infinite current sheet with a current density of 1 A/m. The reference point is given at r(1..3) """ return None # FIXME: not yet implemented def current_sheet(r, s1, s2, s3, finite=True): """Calculate the magnetic field of a current sheet""" if finite: return current_sheet_finite(r, s1, s2, s3) return current_sheet_infinite(r) # ============================================================================ # OBSOLETE # ============================================================================ Loading
test/test_physics_biotsavart.py +58 −4 Original line number Diff line number Diff line #!/usr/bin/env python """ Python unit test cases for mathutil Python unit test cases for Biot-Savart module """ import unittest import numpy as np from numpy import pi, sqrt, cos, allclose from numpy import pi, sqrt, cos, allclose, array from scipy.constants import mu_0 import pysen.physics.biotsavart as bs Loading @@ -18,7 +18,7 @@ L = 1.0 # coil separation in m class BiotSavartTestCase(unittest.TestCase): "Biot-Savart Test Cases" x1,x2 = np.array((-L, 0, 0)), np.array(( L, 0, 0)) # two point segment x1,x2 = array((-L, 0, 0)), array(( L, 0, 0)) # two point segment x = np.linspace(-2*L,2*L, 11) _r = np.ones_like (x)*R Loading Loading @@ -175,7 +175,61 @@ class BiotSavartTestCase(unittest.TestCase): act = bs.solenoid_axis(0, R=R, L=LS, N=N) self.assertAlmostEqual(exp, act, 5) def test_current_sheet_simple(self): "test simple cases for current sheet" xs1 = array((-1, 0,-1)) xs2 = array((-1, 0,+1)) xs3 = array((+1, 0,-1)) r = array((0.1,0.2,0.3)) bexp = array((-0.40581227, 0.02079833, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) xs1 = array((-100, 0,-100)) xs2 = array((-100, 0,+100)) xs3 = array((+100, 0,-100)) r = array((0.0,2e-6,0.0)) bexp = array((-0.5, 0, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) r = array((0.0,-2e-6,0.0)) bexp = array((0.5, 0, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) xs1 = array((-2, 0,-1)) xs2 = array((-2, 0,+1)) xs3 = array((+1, 0,-1)) r = array((0.4,0.6,0.5)) bexp = array((-0.23694913, 0.08906743, 0)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) def test_current_sheet_rotated(self): "test rotated sheet" xs1 = array((-1, 0, 0)) xs2 = array(( 0, 0,+1)) xs3 = array(( 0, 0,-1)) r = array((0.1,0.2,0.3)) bexp = array((-0.255651, -0.0390186, 0.255651)) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) def xtest_current_sheet_vectorized(self): "test vectorized current sheet" xs1 = array((-1, 0,-1)) xs2 = array((-1, 0,+1)) xs3 = array((+1, 0,-1)) r = np.array([[0.1,0.2,0.3], [0.0,2e-6,0.0], [0.0,-2e-6,0.0], [0.4,0.6,0.5]]) bexp = np.array([[-0.40581227, 0.02079833, 0], [-0.5 , 0 , 0], [ 0.5 , 0 , 0], [-0.23694913, 0.08906743, 0]]) bact = bs.current_sheet(r, xs1, xs2, xs3) self.assertTrue(np.allclose(bexp, bact)) if __name__ == "__main__": unittest.main(exit=False, verbosity=2)
