円電流による磁界

ベクトルポテンシャルの計算

r0 = 2;
r[x_, z_] := Sqrt[x^2 + z^2];
th[x_, z_] := ArcTan[z, x]

パソコン電磁学 足立武彦ら　朝倉書店1986年の(6.60), (6.61)よりベクトルポテンシャルの式がある。ただし、K(k), E(k)->K(k^2), E(k^2)の間違いがあるので注意する。

k[r_, th_] := 
        Sqrt[4*r*r0*
            Sin[th]/(r^2 + r0^2 + 
                  2*r*r0*Sin[th])];
a[r_, th_] := 
        r0*k[r, th]/
            Sqrt[r*r0*Sin[th]]*((2/k[r, th]^2 - 1)*
                EllipticK[k[r, th]^2] - (2/k[r, th]^2)*
                EllipticE[k[r, th]^2]);

一応、(6.57)を直接積分してベクトルポテンシャルを決定したものと一致していることを数値的に確認した。

Integrate[
      Cos[thd]/Sqrt[r^2 + r0^2 - 2*r*r0*Sin[th]*Cos[thd]], {thd, 0, 
        2  Pi}]

a2[r_, 
        th_] := -(Csc[
                th] (-(4 + 
                          (r^2)) EllipticK[((8 r Sin[th])/(4 + (r^2) + 4 r Sin[th]))] (Sqrt[4 + (r^2) - 4 r Sin[th]]) + 
                  EllipticE[((8 r Sin[th])/(4 + (r^2) + 4 r Sin[th]))] (Sqrt[4 + (r^2) - 4 r Sin[th]]) (4 + (r^2) + 
                        4 r Sin[
                            th]) + (Sqrt[4 + (r^2) + 4 r Sin[th]]) (-(4 + (r^2)) EllipticK[-((8 r Sin[th])/(4 + (r^2) - 
                                    4 r Sin[th]))] + 
                        EllipticE[-((8 r Sin[th])/(4 + (r^2) - 
                                    4 r Sin[th]))] (4 + (r^2) - 
                              4 r Sin[
                                  th]))))/(2 r (Sqrt[4 + (r^2) - 4 r Sin[th]]) (Sqrt[4 + (r^2) + 4 r Sin[th]]))

こちらは(6.61)式

gra1 = 
      ContourPlot[a[r[x, z], th[x, z]], {x, 0.001, 20}, {z, -10, 10}, 
        PlotPoints -> 100, ContourShading -> False]

こちらは直接積分した式

ContourPlot[a2[r[x, z], th[x, z]], {x, 0.001, 20}, {z, -10, 10}, 
      PlotPoints -> 100, ContourShading -> False]

磁場の計算

br[r_, th_] = 1/r/Sin[th]*D[Sin[th]*a[r, th], th];

bth[r_, th_] = -1/r*D[r*a[r, th], r];

bx[x_, z_] := 
        br[r[x, z], th[x, z]] Sin[th[x, z]] + 
          bth[r[x, z], th[x, z]] Cos[th[x, z]];
bz[x_, z_] := 
        br[r[x, z], th[x, z]] Cos[th[x, z]] - 
          bth[r[x, z], th[x, z]] Sin[th[x, z]];

Needs["Graphics`PlotField`"];

gra2 = 
      PlotVectorField[{bx[x, z], bz[x, z]}, {x, 0.001, 20}, {z, -10, 10}, 
        ScaleFunction -> (0.5 &)]

Show[gra1, gra2]

若干狭い領域で表示させる。

gra1 = 
      ContourPlot[a[r[x, z], th[x, z]], {x, 0.0001, 4}, {z, -2, 2}, 
        PlotPoints -> 100, ContourShading -> False]

gra2 = 
      PlotVectorField[{bx[x, z], bz[x, z]}, {x, 0.001, 4}, {z, -2, 2}, 
        ScaleFunction -> (0.5 &)]

Show[gra1, gra2]

磁束線の表示

\(\(\( (*\ 
      x軸方向にどのくらい移動するか\ *) \)\(
\)\(\(calx[x_, y_] := 
        Module[{length, ret = 0}, 
	
          length = 
            Sqrt[bx[x, y]^2 + 
                bz[x, y]^2]; 
\ \ \ \ \ \ \ \ \ ret = 
            x + bx[x, y]/length*dr; 
	
          Return[ret];
];\)
 (*\ 
      y軸方向にどのくらい移動するか\ *) 

    \(caly[x_, y_] := 
        Module[{r, rdx, rdy, ret = 0}, 
	
          r = Sqrt[
              bx[x, y]^2 + 
                bz[x, y]^2]; 
\ \ \ \ \ \ \ \ \ ret = 
            y + bz[x, y]/r*dr; 
	
          Return[ret];
];\)\)\)\)

limit[x_, y_] := 
        Module[{r}, 
Return[
              y < 0  ||  y >  4];];

dr = 0.01;
\(cal[] := 
      Module[{x, 
          y}, 
li = {}; 	
For[xs = 1.8, 
            xs ≥ 0.2, xs -= 0.2, 
	
            x = xs; 
	y = 0.0; 
	
            tsl = {{x, 
                  y}}; \  (*\ 上からの計算*) 
	
            tsr = {{x, \(-y\)}};  (*\ 下からの計算*) 
	
            While[Not[limit[x, y]], 
	
              x = calx[x, y]; 
	
              y = caly[x, y]; 
	
	
              tsl = Append[tsl, {x, y}]; 
	
              tsr = Append[
                  tsr, {x, \(-y\)}];
\ \ \ \ \ \ ]; 
	
            li = Append[li, {Line[tsl], Line[tsr]}]; 
	
            Print[xs];
\ ]
 (*\ 計算が終わったので作画する\ *) \ 
\ \ Show[
            Graphics[{li}, AspectRatio -> 1, 
              PlotRange -> {{0, 8}, {\(-4\), 
                    4}}]];
]\)

cal[]

gra3 = 
      Graphics[{li}, AspectRatio -> 1, 
        PlotRange -> {{0, 8}, {-4, 4}}]

Show[gra1, gra3]