packages = ["numpy","scipy"]

Magnetic Field and Vector Potential from Rectangular Arc Segment

2025-4-28 (R.2)
Kiyoshi Yoshida
kiyoshi.yoshida@eagle.ocn.ne.jp

Radius a (m)
Radial Width Wr (m)
Axial Length Wz (m)
Current (MA)
Start Angle th1 (deg)
End Angle th2 (deg)
Position x (m)
Position y (m)
Position z (m)
from pyscript import document, display import js from numpy import abs, pi, sqrt, log, arctan, sin, cos, arctan2 from scipy import integrate small_angl = 1.00e-6 small_dist = 1.00e-6 small_value = 1.00e-12 def sgn(x): if x >= 0.0: return 1.0 else: return -1.0 def trans2(uo, vo, t): s = sin(t) c = cos(t) up = c * uo - s * vo vp = s * uo + c * vo return up, vp def ba0x(t1, t2, r1, r2, lg, xp, zp): r0 = 0.5*(r1 + r2) co1 = cos(t1) si1 = sin(t1) si2 = sin(t2) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co1 + zj**2) sg = sgn(z*(r0 - r)) s += sg*rj b = (si2 - si1)*s return b def ba0y(t1, t2, r1, r2, lg, xp, zp): r0 = 0.5 * (r1 + r2) co1 = cos(t1) co2 = cos(t2) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co1 + zj**2) sg = sgn(z*(r0 - r)) s += sg*rj b = -(co2 - co1)*s return b def ba0z(t1, t2, r1, r2, lg, xp, zp): r0 = 0.5 * (r1 + r2) co1 = cos(t1) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co1 + zj**2) sg = sgn(z*(r0 - r)) s += sg*zj*log(rj + r) b = (t1 - t2)*s return b def bax(t, r1, r2, lg, xp, zp): r0 = 0.5 * (r1 + r2) if abs(t) < small_angl: t = small_angl co = cos(t) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg if abs(z) < small_dist: z = small_dist # to avoid divergence zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co + zj**2) sg = sgn(z*(r0 - r)) s += sg*(rj + xp*co*log(rj + r - xp*co))*co return s def bay(t, r1, r2, lg, xp, zp): r0 = 0.5 * (r1 + r2) if abs(t) < small_angl: t = small_angl co = cos(t) si = sin(t) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg if abs(z) < small_dist: z = small_dist # to avoid divergence zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co + zj*zj) sg = sgn(z*(r0 - r)) s += sg*(rj + xp*co*log(rj + r - xp*co))*si return s def baz(t, r1, r2, lg, xp, zp): r0 = 0.50 * (r1 + r2) if abs(t) < small_angl: t = small_angl co = cos(t) si = sin(t) s = 0.0 if abs(xp) < small_dist: xp = small_dist for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg if abs(z) < small_dist: z = small_dist # to avoid divergence zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co + zj*zj) a = zj*(r - xp*co)/(rj*xp*si) sg = sgn(z*(r0 - r)) s += sg*(xp*si*arctan(a) - zj*log(rj + r - xp*co)+xp*co*log(rj + zj)) return s def b_arc(xp, yp, zp, r1, r2, lg, cd, zh, ts, te): pi2 = pi * 2.0 cc = -1.0e-7*cd # x MU/4PI, adjust right-hand screw low zq = zp - zh # convert grovel to local coordination rq = sqrt(xp**2 + yp**2) if r1 > r2: print('**ERROR at b_arc(), r1 < r2: ', r1, r2) return 0.0, 0.0, 0.0 if abs(xp) < small_dist and abs(yp) < small_dist: tha = 0.0 else: tha = arctan2(yp, xp) t1 = ts - tha t2 = te - tha if t1 > pi2 or t2 > pi2: t1 = t1 - pi2 t2 = t2 - pi2 if t1 < -pi2 or t2 < -pi2: t1 = t1 + pi2 t2 = t2 + pi2 if rq < small_dist: bxx = cc*ba0x(t1, t2, r1, r2, lg, xp, zp) byy = cc*ba0y(t1, t2, r1, r2, lg, xp, zp) bzz = cc*ba0z(t1, t2, r1, r2, lg, xp, zp) else: b, e = integrate.quad(bax, t1, t2, args=(r1, r2, lg, rq, zq)) bx = cc * b if abs(t1 - t2) >= pi2: by = 0.0 else: b, e = integrate.quad(bay, t1, t2, args=(r1, r2, lg, rq, zq)) by = cc * b b, e = integrate.quad(baz, t1, t2, args=(r1, r2, lg, rq, zq)) bzz = cc*b bxx, byy = trans2(bx, by, tha) return bxx, byy, bzz def aa0x(t1, t2, r1, r2, lg, xp, zp): r0 = 0.50*(r1 + r2) if abs(t1) < small_angl: t1 = small_angl if abs(t2) < small_angl: t2 = small_angl co1 = cos(t1) co2 = cos(t2) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co1 + zj**2) sg = sgn(z*(r0 - r)) s += sg*(zj*rj + r**2*log(rj + zj)) ax = (co1 - co2)*0.5*s return ax def aa0y(t1, t2, r1, r2, lg, xp, zp): r0 = 0.50*(r1 + r2) if abs(t1) < small_angl: t1 = small_angl if abs(t2) < small_angl: t2 = small_angl co1 = cos(t1) si1 = sin(t1) si2 = sin(t2) s = 0.0 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co1 + zj**2) sg = sgn(z*(r0 - r)) s += sg*(zj*rj + r**2*log(rj + zj)) ay = (si1 - si2)*0.5*s return ay def aax(t, r1, r2, lg, xp, zp): r0 = 0.50 * (r1 + r2) if abs(t) < small_angl: t = small_angl if abs(xp) < small_dist: xp = small_dist co = cos(t) si = sin(t) s = 0.00 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co + zj**2) b = ((xp*co - r)*rj - xp**2 - r**2 + 2.0*xp*r*co)/(zj*xp*si) sg = sgn(z*(r0 - r)) s += sg*(zj*rj*0.50 + zj*xp*co*log(rj + r - xp * co) + (r**2*0.50 + xp**2*0.50 - xp**2*co**2)*log(rj + zj) + xp**2*si*co*arctan(b))*si return s def aay(t, r1, r2, lg, xp, zp): r0 = 0.50 * (r1 + r2) if abs(t) < small_angl: t = small_angl if abs(xp) < small_dist: xp = small_dist co = cos(t) si = sin(t) s = 0.00 for i in range(4): if i == 0: r = r1; z = -0.5*lg elif i == 1: r = r1; z = 0.5*lg elif i == 2: r = r2; z = 0.5*lg else: r = r2; z = -0.5*lg zj = zp - z rj = sqrt(xp**2 + r**2 - 2.0*xp*r*co + zj**2) b = ((xp*co - r)*rj - xp**2 - r**2 + 2.0*xp*r*co)/(zj*xp*si) sg = sgn(z*(r0 - r)) s += sg*(zj*rj*0.50 + zj*xp*co*log(rj + r - xp*co) + (r**2*0.50 + xp**2*0.50 - xp**2*co**2)*log(rj + zj) + xp**2*si*co*arctan(b))*co return -s def a_arc(xp, yp, zp, r1, r2, lg, cd, zh, ts, te): pi2 = pi*2.0 zq = zp - zh # convert grovel to local coordination cc = -1.0e-7*cd # x MU/4PI, adjust right-hand screw low if r1 > r2: print('**ERROR at a_arc(), r1 < r2: ', r1, r2) return 0.0, 0.0 if abs(xp) < small_dist and abs(yp) < small_dist: tha = 0.0 else: tha = arctan2(yp, xp) t1 = ts - tha t2 = te - tha if t1 > pi2 or t2 > pi2: t1 = t1 - pi2 t2 = t2 - pi2 if t1 < -pi2 or t2 < -pi2: t1 = t1 + pi2 t2 = t2 + pi2 rq = sqrt(xp ** 2 + yp ** 2) if rq < small_dist: axx = cc*aa0x(t1, t2, r1, r2, lg, xp, zp) ayy = cc*aa0y(t1, t2, r1, r2, lg, xp, zp) else: a, e = integrate.quad(aax, t1, t2, args=(r1, r2, lg, rq, zq)) ax = cc*a a, e = integrate.quad(aay, t1, t2, args=(r1, r2, lg, rq, zq)) ay = cc*a axx, ayy = trans2(ax, ay, tha) return axx, ayy def bacalc(event): ra = float(js.document.getElementById("radius").value) wr = float(js.document.getElementById("wr").value) wz = float(js.document.getElementById("wz").value) ia = float(js.document.getElementById("current").value)*1.0e6 t1 = float(js.document.getElementById("th1").value)*pi/180.0 t2 = float(js.document.getElementById("th2").value)*pi/180.0 xp = float(js.document.getElementById("xp").value) yp = float(js.document.getElementById("yp").value) zp = float(js.document.getElementById("zp").value) r1 = ra - wr*0.5 r2 = ra + wr*0.5 lg = wz ja = ia/(wr*wz) zh = 0.0 bx, by, bz = b_arc(xp, yp, zp, r1, r2, lg, ja, zh, t1, t2) ax, ay = a_arc(xp, yp, zp, r1, r2, lg, ja, zh, t1, t2) js.document.getElementById("msg1").innerHTML = "Bx (T): "+ str(bx) js.document.getElementById("msg2").innerHTML = "By (T): "+ str(by) js.document.getElementById("msg3").innerHTML = "Bz (T): "+ str(bz) js.document.getElementById("msg4").innerHTML = "Ax (Tm): "+ str(ax) js.document.getElementById("msg5").innerHTML = "Ay (Tm): "+ str(ay)