vectorized version of Circular current

Author: OrtnerMichael

CHANGELOG.md

@@ -10,10 +10,10 @@ All notable changes to magpylib are documented here.
 
 ### Changed
 - Improved performance of getB for diametral magnetized Cylinders by 20%
+- IMPORTANT: position arguments of `getBv` functions have been flipped! First comes the source position POSm THEN the observer position POSo!
 
 ### Added
-- getBv_Cylinder
-- added private vectorized implementation ellipticV to mathLib_vector
+- completed the library vector functionality adding magnet Cylinder, moment Dipole, current Circular and Line. This includes adding several private vectorized functions (e.g. ellipticV) to mathLib_vector, adding respective tests and docu examples.
 
 ## [2.1.0b] - 2019-12-06
 

magpylib/_lib/fields/Current_CircularLoop_vector.py

@@ -0,0 +1,80 @@
+# -------------------------------------------------------------------------------
+# magpylib -- A Python 3 toolbox for working with magnetic fields.
+# Copyright (C) Silicon Austria Labs, https://silicon-austria-labs.com/,
+#               Michael Ortner <magpylib@gmail.com>
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+# PARTICULAR PURPOSE.  See the GNU Affero General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Affero General Public License along
+# with this program.  If not, see <https://www.gnu.org/licenses/>.
+# The acceptance of the conditions of the GNU Affero General Public License are
+# compulsory for the usage of the software.
+#
+# For contact information, reach out over at <magpylib@gmail.com> or our issues
+# page at https://www.github.com/magpylib/magpylib/issues.
+# -------------------------------------------------------------------------------
+
+from numpy import sqrt, array, cos, sin, NaN, arctan2, empty, pi
+from magpylib._lib.mathLib_vector import ellipticKV, ellipticEV, angleAxisRotationV_priv
+from warnings import warn
+import numpy as np
+
+
+# %% CIRCULAR CURRENT LOOP
+# Describes the magnetic field of a circular current loop that lies parallel to
+#    the x-y plane. The loop has the radius r0 and carries the current i0, the
+#    center of the current loop lies at posCL
+
+# i0   : float  [A]     current in the loop
+# d0   : float  [mm]    diameter of the current loop
+# posCL: arr3  [mm]     Position of the center of the current loop
+
+# VECTORIZATION
+
+def Bfield_CircularCurrentLoopV(I0, D, POS):
+
+    R = D/2  #radius
+
+    N = len(D)  # vector size
+
+    X,Y,Z = POS[:,0],POS[:,1],POS[:,2]
+
+    RR, PHI = sqrt(X**2+Y**2), arctan2(Y, X)      # cylindrical coordinates
+
+
+    deltaP = sqrt((RR+R)**2+Z**2)
+    deltaM = sqrt((RR-R)**2+Z**2)
+    kappa = deltaP**2/deltaM**2
+    kappaBar = 1-kappa
+
+    # allocate solution vector
+    field_R = empty([N])
+
+    # avoid discontinuity of Br on z-axis
+    maskRR0 = RR == np.zeros([N])
+    field_R[maskRR0] = 0
+    # R-component computation
+    notM = np.invert(maskRR0)
+    field_R[notM] = -2*1e-4*(Z[notM]/RR[notM]/deltaM[notM])*(ellipticKV(kappaBar[notM]) -
+                        (2-kappaBar[notM])/(2-2*kappaBar[notM])*ellipticEV(kappaBar[notM]))
+    
+    # Z-component computation
+    field_Z = -2*1e-4*(1/deltaM)*(-ellipticKV(kappaBar)+(2-kappaBar -
+                                                      4*(R/deltaM)**2)/(2-2*kappaBar)*ellipticEV(kappaBar))
+
+    # transformation to cartesian coordinates
+    Bcy = np.array([field_R,np.zeros(N),field_Z]).T
+    AX = np.zeros([N,3])
+    AX[:,2] = 1
+    Bkart = angleAxisRotationV_priv(PHI/pi*180,AX,Bcy)
+
+    return (Bkart.T * I0).T * 1000 # to mT
+

magpylib/_lib/fields/PM_Cylinder_vector.py

@@ -169,23 +169,8 @@ def Bfield_CylinderV(MAG, POS, DIM, Nphi0):  # returns arr3
         mask0Inside = mask0 * (RR<R) * (abs(Z)<H/2)
         Bfield[mask0Inside,:2] += MAG[mask0Inside,:2]
 
-        print(mask0Inside)
     ### END TRANS MAG CONTRIBUTION ###########################
 
     return(Bfield)
 
 
-MAG = np.array([[0,0,-44],[0,0,55],[11,22,33],[-14,25,36],[17,-28,39],[-10,-21,32],[0,12,23],[0,-14,25],[16,0,27],[-18,0,29]])
-POS=MAG*0.1*np.array([.8,-1,-1.3])
-DIM = np.array([[1,1]]*10)
-
-Bv = Bfield_CylinderV(MAG,POS,DIM,50)
-print(Bv)
-
-MAG = np.array([[0,0,1],[0,1,0],[1,0,0],[0,1,1],[1,0,1],[1,1,0],[1,1,1]])
-POS = np.zeros([7,3])-.1
-DIM = np.ones([7,2])
-
-Bv = Bfield_CylinderV(MAG,POS,DIM,50)
-
-

magpylib/_lib/getBvector.py

@@ -25,13 +25,13 @@
 from magpylib._lib.fields.PM_Cylinder_vector import Bfield_CylinderV
 from magpylib._lib.fields.PM_Sphere_vector import Bfield_SphereV
 from magpylib._lib.fields.Moment_Dipole_vector import Bfield_DipoleV
+from magpylib._lib.fields.Current_CircularLoop_vector import Bfield_CircularCurrentLoopV
 from magpylib._lib.mathLib_vector import QconjV, QrotationV, QmultV, getRotQuatV, angleAxisRotationV_priv
 import numpy as np
 
-def getBv_magnet(type,MAG,DIM,POSo,POSm,ANG=[],AX=[],ANCH=[],Nphi0=50):
+def getBv_magnet(type,MAG,DIM,POSm,POSo,ANG=[],AX=[],ANCH=[],Nphi0=50):
     """
     Calculate the field of magnets using vectorized performance code.
-    This function is faster than getB when more than ~5 fields are evaluated.
 
     Parameters
     ----------
@@ -106,16 +106,77 @@ def getBv_magnet(type,MAG,DIM,POSo,POSm,ANG=[],AX=[],ANCH=[],Nphi0=50):
 # -------------------------------------------------------------------------------
 
 
-def getBv_current():
-    return 0
+def getBv_current(type,CURR,DIM,POSm,POSo,ANG=[],AX=[],ANCH=[]):
+    """
+    Calculate the field of currents using vectorized performance code.
+
+    Parameters
+    ----------
+
+    type : string
+        source type either 'circular' or 'line'
+
+    MAG : Nx3 numpy array float [mT]
+        vector of N magnetizations.
+
+    DIM : NxY numpy array float [mm]
+        vector of N dimensions for each evaluation. The form of this vector depends
+        on the source type. Y=3/2/1 for box/cylinder/sphere
+
+    POSo : Nx3 numpy array float [mm]
+        vector of N positions of the observer.
+    
+    POSm : Nx3 numpy array float [mm]
+        vector of N initial source positions. These positions will be adjusted by
+        the given rotation parameters.
+
+    ANG=[] : length M list of size N numpy arrays float [deg]
+       Angles of M subsequent rotation operations applied to the N-sized POSm and
+       the implicit source orientation.
+    
+    AX=[] : length M list of Nx3 numpy arrays float []
+        Axis vectors of M subsequent rotation operations applied to the N-sized
+        POSm and the implicit source orientation.
+    
+    ANCH=[] : length M list of Nx3 numpy arrays float [mm]
+        Anchor positions of M subsequent rotations applied ot the N-sized POSm and
+        the implicit source orientation.
+    """
+    N = len(POSo)
+
+    Q = np.zeros([N,4])
+    Q[:,0] = 1              # init orientation
+    
+    Pm = POSm               #initial position
+
+    #apply rotation operations
+    for ANGLE,AXIS,ANCHOR in zip(ANG,AX,ANCH):
+        Q = QmultV(getRotQuatV(ANGLE,AXIS),Q)
+        Pm = angleAxisRotationV_priv(ANGLE,AXIS,Pm-ANCHOR)+ANCHOR
+
+    # transform into CS of source
+    POSrel = POSo-Pm        #relative position
+    Qc = QconjV(Q)          #orientierung
+    POSrot = QrotationV(Qc,POSrel)  #rotation der pos in das CS der Quelle
+
+    # calculate field
+    if type == 'circular':
+        Brot = Bfield_CircularCurrentLoopV(CURR, DIM, POSrot)
+    else:
+        print('Bad type')
+        return 0
+    
+    # transform back
+    B = QrotationV(Q,Brot)  #rückrotation des feldes
+    
+    return B
 
 # -------------------------------------------------------------------------------
 # -------------------------------------------------------------------------------
 
-def getBv_moment(type,MOM,POSo,POSm,ANG=[],AX=[],ANCH=[]):
+def getBv_moment(type,MOM,POSm,POSo,ANG=[],AX=[],ANCH=[]):
     """
     Calculate the field of magnetic moments using vectorized performance code.
-    This function is faster than getB when more than ~5 fields are evaluated.
 
     Parameters
     ----------

magpylib/_lib/mathLib_vector.py

@@ -210,6 +210,7 @@ def anglesFromAxisV(AXIS):
     return np.array([PHI, TH]).transpose()
 
 
+
 def angleAxisRotationV(POS,ANG,AXIS,ANCHOR):
     """
     This is the vectorized version of angleAxisRotation(). Each entry 
@@ -258,8 +259,12 @@ def angleAxisRotationV(POS,ANG,AXIS,ANCHOR):
     return POSnew
 
 
-# vectorized version of elliptic integral
+
 def ellipticV(kc,p,c,s):
+    '''
+    vectorized version of the elliptic integral
+    original implementation from paper Kirby
+    '''
 
     #if kc == 0:
     #    return NaN
@@ -319,4 +324,39 @@ def ellipticV(kc,p,c,s):
         #   entries are reiterated
         mask = (np.abs(g-k) > g*errtol)
 
-    return(np.pi/2)*(ss+cc*em)/(em*(em+pp))
+    return(np.pi/2)*(ss+cc*em)/(em*(em+pp))
+
+
+
+def ellipticKV(x):
+    '''
+    special case complete elliptic integral of first kind ellipticK
+    0 <= x <1
+    '''
+    N = len(x)
+    onez = np.ones([N])
+    return ellipticV((1-x)**(1/2.), onez, onez, onez)
+
+
+
+def ellipticEV(x):
+    '''
+    special case complete elliptic integral of second kind ellipticE
+    E(x) = int_0^pi/2 (1-x sin(phi)^2)^(1/2) dphi
+    requires x < 1 ! 
+    '''
+    N = len(x)
+    onez = np.ones([N])
+    return ellipticV((1-x)**(1/2.), onez, onez, 1-x)
+
+
+
+def ellipticPiV(x, y):
+    '''
+    special case complete elliptic integral of third kind ellipticPi
+    E(x) = int_0^pi/2 (1-x sin(phi)^2)^(1/2) dphi
+    requires x < 1 ! 
+    '''
+    N = len(x)
+    onez = np.ones([N])
+    return ellipticV((1-y)**(1/2.), 1-x, onez, onez)

tests/test_MathVector.py

@@ -4,6 +4,8 @@
 from magpylib._lib.mathLib import getRotQuat, axisFromAngles, anglesFromAxis, angleAxisRotation, elliptic
 from magpylib.math import randomAxisV
 import numpy as np
+from magpylib._lib.mathLib import ellipticK, ellipticE, ellipticPi
+from magpylib._lib.mathLib_vector import ellipticKV, ellipticEV, ellipticPiV
 
 
 def test_QV():
@@ -105,4 +107,22 @@ def test_ellipticV():
     #vector solution
     solV = ellipticV(INP[:,0],INP[:,1],INP[:,2],INP[:,3])
 
-    assert np.amax(abs(solC-solV)) < 1e-10
+    assert np.amax(abs(solC-solV)) < 1e-10
+
+
+
+def test_ellipticSpecialCases():
+    X = np.linspace(0,.9999,33)
+    Y = np.linspace(0,.9999,33)[::-1]
+
+    solV = ellipticKV(X)
+    solC = np.array([ellipticK(x) for x in X])
+    assert np.amax(abs(solV-solC)) < 1e-15
+
+    solV = ellipticEV(X)
+    solC = np.array([ellipticE(x) for x in X])
+    assert np.amax(abs(solV-solC)) < 1e-15
+
+    solV = ellipticPiV(X,Y)
+    solC = np.array([ellipticPi(x,y) for x,y in zip(X,Y)])
+    assert np.amax(abs(solV-solC)) < 1e-15

tests/test_magpyVector.py

@@ -3,6 +3,7 @@
 import numpy as np
 from magpylib.source.magnet import Box, Cylinder, Sphere
 from magpylib.source.moment import Dipole
+from magpylib.source.current import Circular
 from magpylib.vector import getBv_magnet, getBv_current, getBv_moment
 from magpylib.math import axisFromAngles
 from magpylib.math import angleAxisRotationV
@@ -59,7 +60,7 @@ def getB(phi,th,psi):
 
     # vector
     DIM = np.array([dim]*NN)
-    Bv = getBv_magnet('box',MAG,DIM,POSo,POSm,[ANG1,ANG2],[AX1,AX2],[ANCH1,ANCH2])
+    Bv = getBv_magnet('box',MAG,DIM,POSm,POSo,[ANG1,ANG2],[AX1,AX2],[ANCH1,ANCH2])
     Bv = Bv.reshape([Nth,Npsi,Nphi,3])
 
     # assert
@@ -78,7 +79,7 @@ def getB2(phi,th,psi):
 
     # vector
     DIM2 = np.array([dim2]*NN)
-    Bv = getBv_magnet('sphere',MAG,DIM2,POSo,POSm,[ANG1,ANG2],[AX1,AX2],[ANCH1,ANCH2])
+    Bv = getBv_magnet('sphere',MAG,DIM2,POSm,POSo,[ANG1,ANG2],[AX1,AX2],[ANCH1,ANCH2])
     Bv = Bv.reshape([Nth,Npsi,Nphi,3])
 
     #assert
@@ -92,7 +93,7 @@ def test_vectorMagnetCylinder():
     POSO = MAG*0.1*np.array([.8,-1,-1.3])+POSM
     DIM = np.ones([10,2])
 
-    Bv = getBv_magnet('cylinder',MAG,DIM,POSO,POSM)
+    Bv = getBv_magnet('cylinder',MAG,DIM,POSM,POSO)
 
     Bc = []
     for mag,posM,posO,dim in zip(MAG,POSM,POSO,DIM):
@@ -109,7 +110,7 @@ def test_vectorMagnetCylinder():
     DIM = np.ones([7,2])
     POSM = np.zeros([7,3])
 
-    Bv = getBv_magnet('cylinder',MAG,DIM,POSO,POSM,Nphi0=11)
+    Bv = getBv_magnet('cylinder',MAG,DIM,POSM,POSO,Nphi0=11)
 
     Bc = []
     for mag,posM,posO,dim in zip(MAG,POSM,POSO,DIM):
@@ -127,7 +128,7 @@ def test_vectorMomentDipole():
     POSM = np.ones([10,3])
     POSO = MOM*0.1*np.array([.8,-1,-1.3])+POSM
 
-    Bv = getBv_moment('dipole',MOM,POSO,POSM)
+    Bv = getBv_moment('dipole',MOM,POSM,POSO)
 
     Bc = []
     for mom,posM,posO in zip(MOM,POSM,POSO):
@@ -136,3 +137,22 @@ def test_vectorMomentDipole():
     Bc = np.array(Bc)
 
     assert np.amax(abs(Bv-Bc)) < 1e-15
+
+
+
+def test_vectorCurrentCircular():
+    
+    I = np.ones([10])
+    D = np.ones([10])*4
+    Pm = np.zeros([10,3])
+    Po = np.array([[0,0,1],[0,0,-1],[1,1,0],[1,-1,0],[-1,-1,0],[-1,1,0],[5,5,0],[5,-5,0],[-5,-5,0],[-5,5,0]])
+
+    Bc = []
+    for i,d,pm,po in zip(I,D,Pm,Po):
+        s = Circular(curr=i,dim=d,pos=pm)
+        Bc += [s.getB(po)]
+    Bc = np.array(Bc)
+
+    Bv = getBv_current('circular',I,D,Pm,Po)
+
+    assert np.amax(abs(Bc-Bv))<1e-10