BAwirearc

packages = ["numpy","scipy"]

Magnetic Field and Vector Potential from Circular Wire Arc

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

Current		(MA)
Radius		(m)
Start th1		(deg)
End th1		(deg)
Position X		(m)
Position Y		(m)
Position z		(m)

from pyscript import document, display import js from scipy.special import ellipkinc, ellipeinc from numpy import abs, pi, sqrt, sin, cos, arctan2 small_dist = 1.00e-6 def ba_wirarc(x, y, z, ia, ra, p1, p2): # Babic r = sqrt(x**2 + y**2) if r == ra: return 0.0, 0.0, 0.0, 0.0, 0.0 if abs(r) < small_dist: r = small_dist pp = arctan2(y, x) alp1 = (pi + pp - p1)*0.50 alp2 = (pi + pp - p2)*0.50 ca = 1.0e-7*ia a2 = (ra + r)**2 + z**2 k2 = 4.0*ra*r/a2 k1 = sqrt(k2) kd2 = 1.0 - k2 ke2 = k2 - 2.0 fe1 = ellipkinc(alp1, k2) # Incomplete elliptic integral first fe2 = ellipkinc(alp2, k2) fed = fe2 - fe1 ee1 = ellipeinc(alp1, k2) # Incomplete elliptic integral second ee2 = ellipeinc(alp2, k2) eed = ee2 - ee1 ih1 = sqrt(1.0 - k2*sin(alp1)**2) ih2 = sqrt(1.0 - k2*sin(alp2)**2) ihd = ih2 - ih1 ihe = 1.0/ih2 - 1.0/ih1 ig1 = sin(alp1)*cos(alp1)/ih1 ig2 = sin(alp2)*cos(alp2)/ih2 igd = ig2 - ig1 kfe = ke2*fed + 2.0*eed ia = sqrt(ra/r)/(k1*r) ax = -ca*ia*(y*kfe + 2.0*x*ihd) ay = ca*ia*(x*kfe - 2.0*y*ihd) kef = ke2*eed + 2.0*kd2*fed - k2*ke2*igd ib = z*k1/(4.0*r**2*sqrt(ra*r)*kd2) bx = ca*ib*(x*kef + 2.0*y*kd2*ihe) by = ca*ib*(y*kef - 2.0*x*kd2*ihe) ibz = (k2*(ra + r) - 2.0*r)*eed + 2.0*r*kd2*fed + k2*(2.0*r-(ra+r)*k2)*igd bz = -ca*k1/(4.0*r*sqrt(ra*r)*kd2)*ibz return bx, by, bz, ax, ay # main() def bacalc(x): #display("bacalc") ia = float(js.document.getElementById("current").value)*1.0e6 ra = float(js.document.getElementById("radius").value) p1 = float(js.document.getElementById("th1").value)*pi/180.0 p2 = 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) bx, by, bz, ax, ay = ba_wirarc(xp, yp, zp, ia, ra, p1, p2) 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)

Questions to Kiyoshi Yoshida
以下にPython Codeを示す。

# BA Wire Arc
# 2025-4-25 Kiyoshi Yoshida
from scipy.special import ellipkinc, ellipeinc
from numpy import abs, pi, sqrt, sin, cos, arctan2
small_dist = 1.00e-6

def main():
  print("Magnetic Field and Vector Potential of Wire Arc Current")
  ia0=1; ra=0.5; p01=0; p02=180.0; xp=0.01; yp=0.01; zp=0.01
  while True:
    ia0 = inpfuncf("Current (MA), stop<0", ia0)
    if ia0 <0: break
    ra = inpfuncf("Radius (m)", ra)
    p01 = inpfuncf("Th1 (deg)", p01)
    p02 = inpfuncf("Th2 (deg)", p02)
    xp = inpfuncf("xp (m)", xp)
    yp = inpfuncf("yp (m)", yp)
    zp = inpfuncf("zp (m)", zp)
    ia= ia0*1.0e6
    p1 = p01*pi/180.0
    p2 = p02*pi/180.0
    bx, by, bz, ax, ay = ba_wirarc(xp, yp, zp, ia, ra, p1, p2)
    print("Bx (T): {:.6f}".format(bx))
    print("By (T): {:.6f}".format(by))
    print("Bz (T): {:.6f}".format(bz))
    print("Ax (Tm): {:.6f}".format(ax))
    print("Ay (Tm): {:.6f}".format(ay))

def ba_wirarc(x, y, z, ia, ra, p1, p2):    # Babic
  r = sqrt(x**2 + y**2)
  if r == ra:
    return 0.0, 0.0, 0.0, 0.0, 0.0
  if abs(r) < small_dist:
    r = small_dist
  pp = arctan2(y, x)
  alp1 = (pi + pp - p1)*0.50
  alp2 = (pi + pp - p2)*0.50
  ca = 1.0e-7*ia
  a2 = (ra + r)**2 + z**2
  k2 = 4.0*ra*r/a2
  k1 = sqrt(k2)
  kd2 = 1.0 - k2
  ke2 = k2 - 2.0
  fe1 = ellipkinc(alp1, k2)  # Incomplete elliptic integral first
  fe2 = ellipkinc(alp2, k2)
  fed = fe2 - fe1
  ee1 = ellipeinc(alp1, k2)  # Incomplete elliptic integral second
  ee2 = ellipeinc(alp2, k2)
  eed = ee2 - ee1
  ih1 = sqrt(1.0 - k2*sin(alp1)**2)
  ih2 = sqrt(1.0 - k2*sin(alp2)**2)
  ihd = ih2 - ih1
  ihe = 1.0/ih2 - 1.0/ih1
  ig1 = sin(alp1)*cos(alp1)/ih1
  ig2 = sin(alp2)*cos(alp2)/ih2
  igd = ig2 - ig1
  kfe = ke2*fed + 2.0*eed
  ia = sqrt(ra/r)/(k1*r)
  ax = -ca*ia*(y*kfe + 2.0*x*ihd)
  ay = ca*ia*(x*kfe - 2.0*y*ihd)
  kef = ke2*eed + 2.0*kd2*fed - k2*ke2*igd
  ib = z*k1/(4.0*r**2*sqrt(ra*r)*kd2)
  bx = ca*ib*(x*kef + 2.0*y*kd2*ihe)
  by = ca*ib*(y*kef - 2.0*x*kd2*ihe)
  ibz = (k2*(ra + r) - 2.0*r)*eed + 2.0*r*kd2*fed + k2*(2.0*r-(ra+r)*k2)*igd
  bz = -ca*k1/(4.0*r*sqrt(ra*r)*kd2)*ibz
  return bx, by, bz, ax, ay

# keyboard input float
def inpfuncf(msg,x):
  s=input(msg+" : "+str(x)+" : ")
  if s=='':
    y=x
  else:  
    y=float(s)
  return y

if __name__ == "__main__":
  main()