Bfield and Bgrad for rectangular prism magnets, unit test, and LaTeX …

…document

Codes/Bcube.py

@@ -0,0 +1,146 @@
+import numpy as np
+import sympy as sp
+
+mu0 = 4*np.pi*10**-7
+dim = np.array([1,1,1])
+
+#FOR CUBIC MAGNETS
+
+# def Pd(phi,theta): #goes from global to local
+#     return np.array([
+#         [np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)],
+#         [-np.sin(phi), np.cos(phi), 0],
+#         [np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)]
+#     ])
+
+def Pdsym(phi,theta): #goes from global to local, sympy allows for exact trig values
+    return np.array([
+        [float((sp.cos(theta)*sp.cos(phi)).evalf()), float(-sp.cos(theta)*sp.sin(phi).evalf()), float((sp.sin(theta)).evalf())],
+        [float(sp.sin(phi).evalf()), float(sp.cos(phi).evalf()), 0.0],
+        [float((-sp.sin(theta)*sp.cos(phi)).evalf()), float((sp.sin(theta)*sp.sin(phi)).evalf()), float(sp.cos(theta).evalf())]
+    ])
+
+def Hd_i_prime(r, dims):
+    H = np.zeros((3,3))
+
+    xp, yp, zp = r
+
+    x = np.array([xp + dims[0]/2, xp - dims[0]/2])
+    y = np.array([yp + dims[1]/2, yp - dims[1]/2])
+    z = np.array([zp + dims[2]/2, zp - dims[2]/2])
+
+    lst = np.array([0,1])
+    for i in lst:
+        for j in lst:
+            for k in lst:
+                summa = (-1)**(i+j+k)
+                rijk = np.sqrt(x[i]**2 + y[j]**2 + z[k]**2)
+
+                H[0] += summa * np.array([np.arctan2(y[j]*x[i],z[k]*rijk) + np.arctan2(z[k]*x[i],y[j]*rijk), np.log(z[k] + rijk), np.log(y[j] + rijk)])
+                H[1] += summa * np.array([np.log(z[k] + rijk), np.arctan2(x[i]*y[j],z[k]*rijk) + np.arctan2(z[k]*y[j],x[i]*rijk), np.log(x[i] + rijk)])
+                H[2] += summa * np.array([np.log(y[j] + rijk), np.log(x[i] + rijk), np.arctan2(x[i]*z[k],y[j]*rijk) + np.arctan2(y[j]*z[k],x[i]*rijk)])
+    return H/(4*np.pi)
+
+def B_direct(points, magPos, M, dims, phiThetas):
+    N = len(points)
+    D = len(M)
+
+    B = np.zeros((N,3))
+    for n in range(N):
+        for d in range(D):
+            P = Pdsym(phiThetas[d,0],phiThetas[d,1])
+            r_loc = P @ (points[n] - magPos[d])
+
+            tx = np.heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5)
+            ty = np.heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5)
+            tz = np.heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5)       
+            tm = 2*tx*ty*tz
+
+            B[n] += mu0 * P.T @ (Hd_i_prime(r_loc,dims) @ (P @ M[d]) + tm*P@M[d])
+
+    return B
+
+def Bn_direct(points, magPos, M, norms, dims, phiThetas): #solve Bnorm using analytic formula
+    N = len(points)
+    B = B_direct(points, magPos, M, dims, phiThetas)
+    Bn = np.zeros(N)
+    for n in range(N):
+        Bn[n] = B[n] @ norms[n]
+    return Bn
+
+def gd_i(r_loc, n_i_loc, dims): #for matrix formulation
+    tx = np.heaviside(dims[0]/2 - np.abs(r_loc[0]),0.5)
+    ty = np.heaviside(dims[1]/2 - np.abs(r_loc[1]),0.5)
+    tz = np.heaviside(dims[2]/2 - np.abs(r_loc[2]),0.5)       
+    tm = 2*tx*ty*tz
+
+    return (Hd_i_prime(r_loc,dims).T + tm*np.eye(3)) @ n_i_loc
+
+def Acube(points, magPos, norms, dims, phiThetas):
+    N = len(points)
+    D = len(magPos)
+
+    A = np.zeros((N,3*D))
+    for n in range(N):
+        for d in range(D):
+            P = Pdsym(phiThetas[d,0],phiThetas[d,1])
+            r_loc = P @ (points[n] - magPos[d])
+            n_loc = P @ norms[n]
+
+            g = P.T@gd_i(r_loc,n_loc,dims)#P.T to make g global
+            A[n,3*d] = g[0]
+            A[n,3*d+1] = g[1]
+            A[n,3*d+2] = g[2]
+    return mu0 * A
+
+def Bn_fromMat(points, magPos, M, norms, dims, phiThetas): #solving Bnorm using matrix formulation
+    A = Acube(points, magPos, norms, dims, phiThetas)
+    Ms = np.concatenate(M)
+    return A @ Ms
+
+
+#FOR DIPOLES
+
+def Bdip_direct(points, magPos, m):
+    N = len(points)
+    D = len(m)
+
+    B = np.zeros((N,3))
+    for n in range(N):
+        for d in range(D):
+            r = points[n] - magPos[d]
+            B[n] += mu0/(4*np.pi) * (3*m[d]@r * r/np.linalg.norm(r)**5 - m[d]/np.linalg.norm(r)**3)
+    return B
+
+def Bndip_direct(points, magPos, m, norms): #solve Bnorm using analytic formula
+    N = len(points)
+    B = Bdip_direct(points, magPos, m)
+    Bn = np.zeros(N)
+
+    for n in range(N):
+        Bn[n] = B[n] @ norms[n]
+    return Bn
+
+def g_dip(r,n): #for matrix formulation
+    return mu0/(4*np.pi) * (3*r@n * r/np.linalg.norm(r)**5 - n/np.linalg.norm(r)**3)
+
+def Adip(points, magPos, norms):
+    N = len(points)
+    D = len(magPos)
+
+    A = np.zeros((N,3*D))
+    for n in range(N):
+        for d in range(D):
+            r = (points[n] - magPos[d])            
+            g = g_dip(r,norms[n])
+
+            A[n,3*d] = g[0]
+            A[n,3*d+1] = g[1]
+            A[n,3*d+2] = g[2]
+    return A
+
+def Bndip_fromMat(points, magPos, m, norms): #solve Bnorm using matrix formulation
+    A = Adip(points, magPos, norms)
+    ms = np.concatenate(m)
+    return A @ ms
+

Codes/Bdipole.py

@@ -0,0 +1,56 @@
+import numpy as np
+
+mu0 = 4*np.pi*10**-7
+
+#what to do with square root term?
+def g_d(r_i, n_i):
+    return mu0/(4*np.pi) * (3*r_i@n_i/np.linalg.norm(r_i)**5 * r_i - n_i/np.linalg.norm(r_i)**3)
+
+def g_d_sqrt(r_i, n_i, dphi_i, dtheta_i):
+    return mu0/(4*np.pi) * (3*r_i@n_i/np.linalg.norm(r_i)**5 * r_i - n_i/np.linalg.norm(r_i)**3) * np.sqrt(dphi_i*dtheta_i*np.linalg.norm(n_i))
+
+def constructA(points, norms, magLocs, dims):
+    assert len(points) == len(norms)
+    D = len(magLocs)
+    N = len(points)
+
+    A = np.zeros((N,3*D))
+
+    for d in range(D):
+        r = points - magLocs[d]
+        for i in range(N):
+            gd_i = g_d(r[i], norms[i])
+            A[i][3*d] = gd_i[0]
+            A[i][3*d+1] = gd_i[1]
+            A[i][3*d+2] = gd_i[2]
+
+    return A
+
+def Bnorm_solve(points, norms, magLocs, dims, m):
+    A = constructA(points, norms, magLocs, dims)
+    ms = np.concatenate(m)
+    return A@ms
+
+#solve directly, without A matrix
+def Bnorm_dip_simple(points, dip_locs, m, norms): 
+    D = len(dip_locs)
+    N = len(points)
+
+    B = np.zeros((N,3))
+    Bnorm = np.zeros(N)
+
+    for n in range(N):
+        mag = np.zeros((D,3))
+        r = points[n] - dip_locs
+
+        for d in range(D):
+            mag[d] = (3*m[d]@r[d]/(np.linalg.norm(r[d])**5)) * r[d] - m[d]/(np.linalg.norm(r[d])**3)
+
+        magSum = np.sum(mag,0)
+        B[n] = (mu0 / (4*np.pi)) * magSum
+        Bnorm[n] = B[n]@norms[n]
+
+    return Bnorm, B
+
+print(Bnorm_dip_simple(np.array([[3000,-30,10000]]),np.array([[0,0,0]]),np.array([[0,0,1]]),np.array([[0,0,1]]))[0])
+

Codes/Bgrad.py

@@ -0,0 +1,77 @@
+import numpy as np
+import sympy as sp
+from Bcube__7_1_24 import Pdsym
+
+
+#FOR CUBE
+mu0 = 4*np.pi*10**-7
+dim = np.array([1,1,1])
+
+
+def grad_r_H(r_loc, dims, P):
+
+    dH = np.zeros((3,3,3))
+
+    x = np.array([r_loc[0] + dims[0]/2, r_loc[0] - dims[0]/2])
+    y = np.array([r_loc[1] + dims[1]/2, r_loc[1] - dims[1]/2])
+    z = np.array([r_loc[2] + dims[2]/2, r_loc[2] - dims[2]/2])
+
+    lst = np.array([0,1])
+
+    for s in range(3):
+        for row in range(3):
+            for col in range(3):
+                for i in lst:
+                    for j in lst:
+                        for k in lst:
+                            summa = (-1)**(i+j+k)
+
+                            rp = np.array([x[i],y[j],z[k]])
+                            drp = np.array([P[0][s],P[1][s],P[2][s]])
+
+                            rijk = np.sqrt(rp[0]**2 + rp[1]**2 + rp[2]**2)
+                            drijk = (rp[0]*drp[0]+rp[1]*drp[1]+rp[2]*drp[2])/rijk
+
+                            if row == col:
+                                dH[row,col,s] += summa * ( (rp[row-1]*rijk*(drp[row]*rp[row-2]+rp[row]*drp[row-2]) - rp[row-2]*rp[row]*(drp[row-1]*rijk+rp[row-1]*drijk))/((rp[row-1]*rijk)**2 + (rp[row-2]*rp[row])**2) + (rp[row-2]*rijk*(drp[row]*rp[row-1]+rp[row]*drp[row-1]) - rp[row-1]*rp[row]*(drp[row-2]*rijk+rp[row-2]*drijk))/((rp[row-2]*rijk)**2 + (rp[row-1]*rp[row])**2) )
+                            else:
+                                if 2*row-col > 2:
+                                    index = 2*row-col - 3
+                                else:
+                                    index = 2*row-col
+                                dH[row,col,s] += summa * 1/(rp[index]+rijk) * (drp[index]+drijk)
+    return dH/(4*np.pi)
+
+
+def gradr_Bcube(points, magPos, M, phiThetas, dims):
+    D = len(M)
+    N = len(points)
+
+    dB = np.zeros((N,3,3))
+    for n in range(N):
+        for d in range(D):
+            P = Pdsym(phiThetas[d][0],phiThetas[d][1])
+            r_loc = P @ (points[n] - magPos[d])
+            drH = grad_r_H(r_loc, dims, P)
+            dB[n] += mu0 * P.T @ (drH @ M[d]) @ P
+    return dB
+
+
+#FOR DIPOLE
+
+def gradr_Bdip(points, magPos, m):
+    D = len(m)
+    N = len(points)
+    hat = np.array([[1,0,0],[0,1,0],[0,0,1]])
+
+    dB = np.zeros((N,3,3))
+    for n in range(N):
+        for d in range(D):
+            r = points[n] - magPos[d]
+            rmag = np.linalg.norm(r)
+            for s in range(3):
+                grad = 3*mu0/(4*np.pi) * (((m[d,s]*r+(m[d]@r)*hat[s])*rmag**2 - 5*r[s]*r*(m[d]@r))/rmag**7 + m[d]*r[s]/rmag**5)
+                dB[n,:,s] += grad
+
+    return dB
+