Skip to content

Commit 07c24f0

Browse files
author
Elizabeth
committed
Initial commit.
1 parent cb377be commit 07c24f0

File tree

2 files changed

+96
-17
lines changed

2 files changed

+96
-17
lines changed

‎src/simsopt/field/magneticfieldclasses.py‎

Lines changed: 39 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -245,15 +245,47 @@ def _B_impl(self, B):
245245
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
246246
z = points[:, 2]
247247
phi = np.arctan2(points[:, 1], points[:, 0])
248-
B[:] = np.array(self.Blambdify(r, z, phi)).T
248+
B_cyl = np.array(self.Blambdify(r, z, phi)).T
249+
# Bx = Br cos(phi) - Bphi sin(phi)
250+
B[:, 0] = B_cyl[:, 0] * np.cos(phi) - B_cyl[:, 1] * np.sin(phi)
251+
# By = Br sin(phi) + Bphi cos(phi)
252+
B[:, 1] = B_cyl[:, 0] * np.sin(phi) + B_cyl[:, 1] * np.cos(phi)
253+
B[:, 2] = B_cyl[:, 2]
249254

250255
def _dB_by_dX_impl(self, dB):
251256
points = self.get_points_cart_ref()
252257
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
253-
r = np.sqrt(np.square(points[:, 0]) + np.square(points[:, 1]))
254258
z = points[:, 2]
255259
phi = np.arctan2(points[:, 1], points[:, 0])
256-
dB[:] = np.array(self.dBlambdify_by_dX(r, z, phi)).transpose((2, 0, 1))
260+
dB_cyl = np.array(self.dBlambdify_by_dX(r, z, phi)).transpose((2, 0, 1))
261+
dBrdx = dB_cyl[:, 0, 0]
262+
dBrdy = dB_cyl[:, 1, 0]
263+
dBrdz = dB_cyl[:, 2, 0]
264+
dBphidx = dB_cyl[:, 0, 1]
265+
dBphidy = dB_cyl[:, 1, 1]
266+
dBphidz = dB_cyl[:, 2, 1]
267+
dB[:, 0, 2] = dB_cyl[:, 0, 2]
268+
dB[:, 1, 2] = dB_cyl[:, 1, 2]
269+
dB[:, 2, 2] = dB_cyl[:, 2, 2]
270+
dcosphidx = -points[:, 0]**2/r**3 + 1/r
271+
dsinphidx = -points[:, 0]*points[:, 1]/r**3
272+
dcosphidy = -points[:, 0]*points[:, 1]/r**3
273+
dsinphidy = -points[:, 1]**2/r**3 + 1/r
274+
B_cyl = np.array(self.Blambdify(r, z, phi)).T
275+
Br = B_cyl[:, 0]
276+
Bphi = B_cyl[:, 1]
277+
# Bx = Br cos(phi) - Bphi sin(phi)
278+
dB[:, 0, 0] = dBrdx * np.cos(phi) + Br * dcosphidx - dBphidx * np.sin(phi) \
279+
- Bphi * dsinphidx
280+
dB[:, 1, 0] = dBrdy * np.cos(phi) + Br * dcosphidy - dBphidy * np.sin(phi) \
281+
- Bphi * dsinphidy
282+
dB[:, 2, 0] = dBrdz * np.cos(phi) - dBphidz * np.sin(phi)
283+
# By = Br sin(phi) + Bphi cos(phi)
284+
dB[:, 0, 1] = dBrdx * np.sin(phi) + Br * dsinphidx + dBphidx * np.cos(phi) \
285+
+ Bphi * dcosphidx
286+
dB[:, 1, 1] = dBrdy * np.sin(phi) + Br * dsinphidy + dBphidy * np.cos(phi) \
287+
+ Bphi * dcosphidy
288+
dB[:, 2, 1] = dBrdz * np.sin(phi) + dBphidz * np.cos(phi)
257289

258290

259291
class CircularCoil(MagneticField):
@@ -292,13 +324,13 @@ def __init__(self, r0=0.1, center=[0, 0, 0], I=5e5/np.pi, normal=[0, 0]):
292324

293325
self.rotMatrix = np.array([
294326
[np.cos(phi) * np.cos(theta)**2 + np.sin(theta)**2,
295-
-np.sin(phi / 2)**2 * np.sin(2 * theta),
327+
-np.sin(phi / 2)**2 * np.sin(2 * theta),
296328
np.cos(theta) * np.sin(phi)],
297-
[-np.sin(phi / 2)**2 * np.sin(2 * theta),
298-
np.cos(theta)**2 + np.cos(phi) * np.sin(theta)**2,
329+
[-np.sin(phi / 2)**2 * np.sin(2 * theta),
330+
np.cos(theta)**2 + np.cos(phi) * np.sin(theta)**2,
299331
np.sin(phi) * np.sin(theta)],
300332
[-np.cos(theta) * np.sin(phi),
301-
-np.sin(phi) * np.sin(theta),
333+
-np.sin(phi) * np.sin(theta),
302334
np.cos(phi)]
303335
])
304336

‎tests/field/test_magneticfields.py‎

Lines changed: 57 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -116,16 +116,63 @@ def test_scalarpotential_Bfield(self):
116116
Bscalar.set_points(points)
117117
B1 = np.array(Bscalar.B())
118118
dB1_by_dX = np.array(Bscalar.dB_by_dX())
119+
119120
# Analytical Formula for B
120121
rphiz = [[np.sqrt(np.power(point[0], 2) + np.power(point[1], 2)), np.arctan2(point[1], point[0]), point[2]] for point in points]
121122
B2 = np.array([[0.2*point[2]+0.8*point[0], (0.1+0.3*point[2])/point[0], 0.2*point[0]+0.3*point[1]+point[2]] for point in rphiz])
123+
# Convert to Cartesian coordinates
124+
r = np.sqrt(np.power(points[:, 0], 2) + np.power(points[:, 1], 2))
125+
phi = np.arctan2(points[:, 1], points[:, 0])
126+
z = points[:, 2]
127+
B2_cart = np.zeros_like(B2)
128+
# Bx = Br cos(phi) - Bphi sin(phi)
129+
B2_cart[:, 0] = B2[:, 0] * np.cos(phi) - B2[:, 1] * np.sin(phi)
130+
# By = Br sin(phi) + Bphi cos(phi)
131+
B2_cart[:, 1] = B2[:, 0] * np.sin(phi) + B2[:, 1] * np.cos(phi)
132+
B2_cart[:, 2] = B2[:, 2]
122133
dB2_by_dX = np.array([
123134
[[0.8*np.cos(point[1]), -(np.cos(point[1])/point[0]**2)*(0.1+0.3*point[2]), 0.2*np.cos(point[1])-0.3*np.sin(point[1])/point[0]],
124135
[0.8*np.sin(point[1]), -(np.sin(point[1])/point[0]**2)*(0.1+0.3*point[2]), 0.2*np.sin(point[1])+0.3*np.cos(point[1])/point[0]],
125136
[0.2, 0.3/point[0], 1]] for point in rphiz])
137+
dBxdx = dB1_by_dX[:, 0, 0]
138+
dBxdy = dB1_by_dX[:, 1, 0]
139+
dBxdz = dB1_by_dX[:, 2, 0]
140+
dBydx = dB1_by_dX[:, 0, 1]
141+
dBydy = dB1_by_dX[:, 1, 1]
142+
dBydz = dB1_by_dX[:, 2, 1]
143+
dB1_by_dX_cyl = np.zeros_like(dB2_by_dX)
144+
dcosphidx = -points[:, 0]**2/r**3 + 1/r
145+
dsinphidx = -points[:, 0]*points[:, 1]/r**3
146+
dcosphidy = -points[:, 0]*points[:, 1]/r**3
147+
dsinphidy = -points[:, 1]**2/r**3 + 1/r
148+
Bx = B1[:, 0]
149+
By = B1[:, 1]
150+
# Br = Bx cos(phi) + By sin(phi)
151+
dB1_by_dX_cyl[:, 0, 0] = dBxdx * np.cos(phi) + Bx * dcosphidx + dBydx * np.sin(phi) \
152+
+ By * dsinphidx
153+
dB1_by_dX_cyl[:, 1, 0] = dBxdy * np.cos(phi) + Bx * dcosphidy + dBydy * np.sin(phi) \
154+
+ By * dsinphidy
155+
dB1_by_dX_cyl[:, 2, 0] = dBxdz * np.cos(phi) + dBydz * np.sin(phi)
156+
# Bphi = - sin(phi) Bx + cos(phi) By
157+
dB1_by_dX_cyl[:, 0, 1] = - dBxdx * np.sin(phi) - Bx * dsinphidx + dBydx * np.cos(phi) \
158+
+ By * dcosphidx
159+
dB1_by_dX_cyl[:, 1, 1] = - dBxdy * np.sin(phi) - Bx * dsinphidy + dBydy * np.cos(phi) \
160+
+ By * dcosphidy
161+
dB1_by_dX_cyl[:, 2, 1] = - dBxdz * np.sin(phi) + dBydz * np.cos(phi)
162+
dB1_by_dX_cyl[:, :, 2] = dB1_by_dX[:, :, 2]
126163
# Verify
127-
assert np.allclose(B1, B2)
128-
assert np.allclose(dB1_by_dX, dB2_by_dX)
164+
assert np.allclose(B1, B2_cart)
165+
assert np.allclose(dB1_by_dX_cyl, dB2_by_dX)
166+
167+
# Check for divergence-free condition for dipole field
168+
# Set up magnetic field scalar potential
169+
PhiStr = "Z/(R*R + Z*Z)**(3/2)"
170+
# Set up scalar potential B
171+
Bscalar = ScalarPotentialRZMagneticField(PhiStr)
172+
Bscalar.set_points(points)
173+
dB1_by_dX = np.array(Bscalar.dB_by_dX())
174+
divB = dB1_by_dX[:, 0, 0] + dB1_by_dX[:, 1, 1] + dB1_by_dX[:, 2, 2]
175+
assert np.allclose(np.abs(divB), 0)
129176

130177
def test_circularcoil_Bfield(self):
131178
current = 1.2e7
@@ -284,8 +331,8 @@ def test_circularcoil_Bfield_toroidal_arrangement(self):
284331
N_coils = 30
285332

286333
N_turns = 3
287-
a1 = 10 / 2 * 0.0254
288-
a2 = 19.983 / 2 * 0.0254
334+
a1 = 10 / 2 * 0.0254
335+
a2 = 19.983 / 2 * 0.0254
289336
r_array = np.linspace(a1, a2, N_turns)
290337
I_amp = 433 * (33/N_turns)
291338

@@ -420,16 +467,16 @@ def test_Reiman(self):
420467
x = points[:, 0]
421468
y = points[:, 1]
422469
z = points[:, 2]
423-
Bx = (y*np.sqrt(x**2 + y**2) + x*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
424-
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
470+
Bx = (y*np.sqrt(x**2 + y**2) + x*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
471+
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
425472
0.06*x*(1 - np.sqrt(x**2 + y**2))*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2 *
426473
np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/(x**2 + y**2)
427-
By = (-1.*x*np.sqrt(x**2 + y**2) + y*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
428-
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
474+
By = (-1.*x*np.sqrt(x**2 + y**2) + y*z*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
475+
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2))))) +
429476
0.06*y*(1 - np.sqrt(x**2 + y**2))*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2 *
430477
np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/(x**2 + y**2)
431-
Bz = (-((-1 + np.sqrt(x**2 + y**2))*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
432-
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))) -
478+
Bz = (-((-1 + np.sqrt(x**2 + y**2))*(0.15 + 0.38*((-1 + np.sqrt(x**2 + y**2))**2 + z**2) -
479+
0.06*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.cos(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))) -
433480
0.06*z*((-1 + np.sqrt(x**2 + y**2))**2 + z**2)**2*np.sin(np.arctan2(y, x) - 6*np.arctan(z/(-1 + np.sqrt(x**2 + y**2)))))/np.sqrt(x**2 + y**2)
434481
B2 = np.array(np.vstack((Bx, By, Bz)).T)
435482
assert np.allclose(B1, B2)

0 commit comments

Comments
 (0)