Commit 04465e15 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

documentation fixed in biotsavart

parent d1ff3f78

pysen/physics/biotsavart.py

+74 −71
Original line number Diff line number Diff line
@@ -9,9 +9,10 @@ Note: the original paper by M. Monkenbusch has a typo 
Current loop calculated analytically,

see Batygin, Toptygin, problem 249.





u0 = 4*pi * 10^-7



Note: the result of each function must be multiplied by mu0 and I,

or NI or N/I as noted below in order to get the magnetic field.

"""

import numpy as np

from   numpy import (pi, radians, asarray, ones, zeros_like, linspace, sign, log,
@@ -25,7 +26,7 @@ _four_pi = 4*pi 
_EPSILON  = 1e-6



# ============================================================================

# SPECIAL CASES

# Special (simple) cases

# ============================================================================

def straight_wire(x, R, L):

    """Special case of Biot-Savart for a finite wire section
@@ -41,6 +42,8 @@ def straight_wire(x, R, L): 
        4*pi                 R

    theta_1 = arctan2(R, x+L)

    theta_2 = arctan2(R, L-x)



    Note: the result must be multiplied by mu0 and I to get the magnetic field.

    """

    return 1/(4*pi) * (cos(arctan2(R, L-x)) + cos(arctan2(R, L+x))) / R
@@ -50,6 +53,8 @@ def current_loop_axis(x, R=1): 
        u0*I      R^2

    B = ----- *  ---------------

          2      (R^2+x^2)^(3/2)



    Note: the result must be multiplied by mu0 and I to get the magnetic field.

    """

    return 1/2 * R**2 / sqrt((R**2+x**2)**3)
@@ -59,6 +64,8 @@ def helmholtz_pair_axis(x, R=1, L=None, anti=False): 
        u0*I

    B = ----- * R^2 * { [(x-L/2)^2+R^2]^3/2 +  [(x+L/2)^2+R^2]^3/2}

          2



    Note: the result must be multiplied by mu0 and I to get magnetic field.

    """

    L = L or R

    x1 = sqrt((x-L/2)**2+R**2)**3
@@ -69,20 +76,24 @@ def helmholtz_pair_axis(x, R=1, L=None, anti=False): 
        x = 1/x1 + 1/x2

    return R**2*x/2



def solenoid_axis(x,  R=1, L=None, N=None):

def solenoid_axis(x,  R=1, L=None):

    """

    Field of a finite solenoid on the symmetry axis



        u0*N*I

    B = ----- * {  (L/2-x)/[L SQRT(R^2 + (L/2-x)^2)] +

          2        (L/2+x)/[L SQRT(R^2 + (L/2+x)^2)] }



    Note: the result must be multiplied by mu0 and N*I to get the magnetic field.

    """

    L = L or R

    N = N or 100

    #

    x1 = (L/2-x)/sqrt(R**2 + (L/2-x)**2)

    x2 = (L/2+x)/sqrt(R**2 + (L/2+x)**2)

    return N/2/L*(x1+x2)

    return 1/2/L*(x1+x2)



# ============================================================================

# General (more difficult) cases

# ============================================================================

def biot_savart_section(r, x1, x2):

    """General case of Biot-Savart formula for a straight wire section
@@ -95,7 +106,7 @@ def biot_savart_section(r, x1, x2): 
    Note: one needs to multiply the result by current

    and the appropriate magnetic permeability. E.g. for SI::



    db = (mu_0/4pi * I) * ds x (r-s) /|r-s|^3

    dB = (mu_0/4pi * I) * db

    """

    x1 = asarray(x1)

    x2 = asarray(x2)
@@ -119,12 +130,18 @@ def biot_savart_section(r, x1, x2): 
        db = db[:, np.newaxis]

    return db*xp/_four_pi



def _current_loop(x, r, theta, a):

def current_loop_kernel(x, r, theta, a):

    """Calulates magnetic field due to a current loop.

    The loop is assumed to be symmetric arount x-axis

    x,r,theta - observation point

    a         - loop radius



    See Batygin, Toptygin, problem 249.



    Notes:

    1) Batygin, Toptygin use Gauss units, so their result must be converted

            2/c -> mu0/2pi

    2) The result must be multiplied by mu0 and I to get the magnetic field.

    """

    assert a>0

    assert r>=0
@@ -146,63 +163,28 @@ def _current_loop(x, r, theta, a): 
    f  = 2*pi*a*sqrt(q) # re-scale factor

    return bx/f,by/f,bz/f



_current_loop = np.vectorize(_current_loop)

current_loop_kernel = np.vectorize(current_loop_kernel)



def current_loop(r, radius=1.0, xloop=0.0):

    """Calulates magnetic field due to a current loop.

    The loop is assumed to be symmetric arount x-axis

    r         - observation point

    radius    - loop radius

    xloop     - loop position along the z axis

    wrapper around _current_loop

    xloop     - loop offset along the z axis



    This is a "vectorized" wrapper around current_loop_kernel



    The result must be multiplied by mu0 and I to get the magnetic field.

    """

    #bfield = zeros_like(r)

    r = np.atleast_2d(r)

    x, y, z = r[...,0], r[...,1], r[...,2]

    rho = sqrt(y**2+z**2)

    theta = arctan2(z,y)

    bx, by, bz = _current_loop(x - xloop, rho, theta, radius)

    bx, by, bz = current_loop_kernel(x - xloop, rho, theta, radius)

    return np.squeeze(np.vstack((bx,by,bz)).T)





def helmholtz_pair(r, radius=1, separation=-1, xloop=0, anti=False):

    """Calulates magnetic field due to a Helmholtz pair

    The loop is approximated as a series of straight sections

    The loop is assumed to be perpendicular to z-axis.

    r          - observation point

    radius     - loop radius

    separation - loop separation (default is radius)

    xloop      - loop position along the z axis

    nsegments  - number of segments

    """

    #bfield = zeros_like(r)    # resulting field

    if separation < 0:

        separation = radius



    b1 = current_loop(r, radius=radius, xloop=xloop-separation/2)

    b2 = current_loop(r, radius=radius, xloop=xloop+separation/2)



    if anti:

        return b2-b1

    return b2+b1



def solenoid(r, radius=1, start=-0.5, end=+0.5, nturns=10):

    """Calulates magnetic field due to a solenoid

    The solenoid axis is assumed to be parallel to x-axis.

        r       - observation point

        radius  - solenoid radius

        start   - solenoid start

        xcenter - solenoid end

        nturns  - number of turns

    """

    bfield = zeros_like(r)    # resulting field

    for xloop in linspace(start, end, nturns):

        db = current_loop(r, radius, xloop)

        bfield = bfield + db

    return bfield



# ============================================================================

# CURRENT SHEETS, FLIPPERS

# Current sheets and flippers

# ============================================================================

def bsheet_kernel(r, width, height):

    """
@@ -258,6 +240,8 @@ def current_sheet_finite(r, s1, s2, s3): 
    by the three corners: s1, s2, s3. The current flows from s1 to s2.



    Translated from bfelder.c:bsheet subroutine by M. Monkenbusch



    Note: the result must be multiplied by mu0 and I/N to get the magnetic field.

    """

    ex  = s3 - s1

    ez  = s2 - s1
@@ -286,7 +270,6 @@ def current_sheet_infinite(_r): 
    """

    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:
@@ -294,9 +277,46 @@ def current_sheet(r, s1, s2, s3, finite=True): 
    return current_sheet_infinite(r)



# ============================================================================

# OBSOLETE

# Helmholtz pairs and solenoids

# ============================================================================

def helmholtz_pair(r, radius=1, separation=-1, xloop=0, anti=False):

    """Calulates magnetic field due to a Helmholtz pair



    The loop is assumed to be perpendicular to z-axis.

    r          - observation point

    radius     - loop radius

    separation - loop separation (default is radius)

    xloop      - loop position along the z axis

    """

    #bfield = zeros_like(r)    # resulting field

    if separation < 0:

        separation = radius



    b1 = current_loop(r, radius=radius, xloop=xloop-separation/2)

    b2 = current_loop(r, radius=radius, xloop=xloop+separation/2)



    if anti:

        return b2-b1

    return b2+b1



def solenoid(r, radius=1, start=-0.5, end=+0.5, nturns=10):

    """Calulates magnetic field due to a solenoid

    The solenoid axis is assumed to be parallel to x-axis.

        r       - observation point

        radius  - solenoid radius

        start   - solenoid start

        xcenter - solenoid end

        nturns  - number of turns

    """

    bfield = zeros_like(r)    # resulting field

    for xloop in linspace(start, end, nturns):

        db = current_loop(r, radius, xloop)

        bfield = bfield + db

    return bfield



# ============================================================================

# Obsolete/Testing

# ============================================================================

def current_loop_section(r, radius=1.0, xloop=0.0, nsegments=8):

    """Calulates magnetic field due to a current loop.

    The loop is approximated as a series of straight sections
@@ -318,25 +338,6 @@ def current_loop_section(r, radius=1.0, xloop=0.0, nsegments=8): 
        bfield = bfield + db

    return bfield



def helmholtz_pair_section(r, radius=1, separation=-1, xloop=0, nsegments=8, anti=False):

    """Calulates magnetic field due to a Helmholtz pair

    The loop is approximated as a series of straight sections

    The loop is assumed to be perpendicular to z-axis.

    r         - observation point

    radius    - loop radius

    zloop     - loop position along the z axis

    nsegments - number of segments

    """

    if separation < 0:

        separation = radius



    db1 = current_loop_section(r, radius, xloop-separation/2, nsegments)

    db2 = current_loop_section(r, radius, xloop+separation/2, nsegments)

    if anti:

        return db2 - db1

    return db2 + db1





def solenoid_section(r, radius=1, length=1, xcenter=0, nturns=10, nsegments=8):

    """Calulates magnetic field due to a solenoid

    The solenoid axis is assumed to be to x-axis.
@@ -346,6 +347,8 @@ def solenoid_section(r, radius=1, length=1, xcenter=0, nturns=10, nsegments=8): 
        xcenter - position of the solenoid center along the z-axis

        nturns  - number of

        nsegments - "



    This is used for testing only

    """

    bfield = zeros_like(r)    # resulting field

    for xloop in linspace(xcenter-length/2.0, xcenter+length/2.0, nturns):

test/test_physics_biotsavart.py

+32 −4
Original line number Diff line number Diff line
@@ -58,10 +58,40 @@ class BiotSavartTestCase(unittest.TestCase): 
        act = bs.biot_savart_section(self.r1, self.x1, self.x2)[:,-1]

        self.assertTrue(allclose(exp, act))



    def test_straight_wire_4(self):

        "compare directly to Michael's code"

        #TODO: check units

        s1 = array((-1, 0,-1))

        s2 = array((+1, 0,-1))

        r  = array(( 0.1, 0.2, 0.3))

        bexp = array((0, -0.0721311, 0.0110971))

        bact = bs.biot_savart_section(r, s1, s2)

        self.assertTrue(np.allclose(bexp, bact), f"{bexp} != {bact}")

        #

        s1 = array((-1, 0, 0))

        s2 = array((+1, 0, 0))

        r  = array(( 0.0, 1.0, 0.0))

        bexp = array((0, 0, sqrt(2)/4/pi))

        bact = bs.biot_savart_section(r, s1, s2)

        self.assertTrue(np.allclose(bexp, bact), f"{bexp} != {bact}")

        r  = array(( 0.0,-1.0, 0.0))

        bexp = array((0, 0, -sqrt(2)/4/pi))

        bact = bs.biot_savart_section(r, s1, s2)

        self.assertTrue(np.allclose(bexp, bact), f"{bexp} != {bact}")

        #

        s1 = array((-1, 0, 0))

        s2 = array(( 0, 0,+1))

        r  = array(( 0.1, 0.2, 0.3))

        bexp = array((-0.0447825, 0.17913, 0.0447825))

        bact = bs.biot_savart_section(r, s1, s2)

        self.assertTrue(np.allclose(bexp, bact), f"{bexp} != {bact}")





    def test_current_loop_1(self):

        "current loop in the center"

        exp=bs.current_loop_axis(0, R)

        act=bs.current_loop((0,0,0), radius=R)

        print(exp, act)

        self.assertTrue(allclose(exp, act[0]))



    def test_current_loop_2(self):
@@ -125,7 +155,7 @@ class BiotSavartTestCase(unittest.TestCase): 
        "solenoid field"

        LS=100*R

        N=1000

        exp = bs.solenoid_axis(0,  R=R, L=LS, N=N)

        exp = bs.solenoid_axis(0,  R=R, L=LS)*N

        act = bs.solenoid((0,0,0), radius=R, start=-LS/2, end=+LS/2, nturns=N)

        accuracy = abs(exp-act[0])/exp

        self.assertAlmostEqual(exp, act[0], 1)
@@ -141,8 +171,6 @@ class BiotSavartTestCase(unittest.TestCase): 
            act = bs.solenoid        (r, radius=R, start=-LS/2, end=+LS/2, nturns=N)

            self.assertTrue(allclose(exp, act, atol=1e-3))







    def test_absolute_values_segment(self):

        "wire segment value"

        exp = 1/(2*pi)*cos(pi/4)/R
@@ -173,7 +201,7 @@ class BiotSavartTestCase(unittest.TestCase): 
        N   = 1000

        LS  = 1000*R

        exp = N/LS

        act = bs.solenoid_axis(0, R=R, L=LS, N=N)

        act = bs.solenoid_axis(0, R=R, L=LS)*N

        self.assertAlmostEqual(exp, act, 5)



    def test_current_sheet_simple(self):