松下先生電磁気教科書
図1.20
電気双極子の等電位面と電気力線の様子

ここでは教科書図1.20にある電気双極子の等電位面と電気力線の様子を計算した。
Er, Ethをr, thで表すようにした。

等電位線

a = 
        0.01;   
dr = 
        a*2.5;    
ep0 = 1;   
Q1 = 
        1;   

Clear[v];

v[r_, th_] := 
      Q1 a Cos[
            th]/(4 Pi ep0 r^2);     
r[x_, z_] := Sqrt[x^2 + z^2]; 
    th[x_, z_] := 
      Which[x > 0, ArcTan[x/z], z < 0, Pi + ArcTan[x/z], x == 0, 
        0];

次の等電位線の方向を探る。1度ずつ角度を振ってdrの距離進み、そのときに電圧が等しいところを探す。そしてその角度を返す。

caltheta[x_, z_] := 
      Module[{theta, ret = 0}, 
vst = 
          v[r[x, z], th[x, z]]; 
For[theta = 0, 
          theta > -Pi, 
          theta -= Pi/360, 
xnew = 
            x + dr*Cos[theta]; 
znew = 
            z + dr*Sin[theta]; 
If[
            Abs[(v[r[xnew, znew], th[xnew, znew]] - vst)/vst] < 
              0.005, Break[];];
]; 

 Return[theta];
]

範囲をはみ出したところで計算をストップする。

limit[x_, z_] := 
      Module[{}, 
Return[(x > 1) || (x < 
                  0) || (z > 1) || (z < 
                  0)];
]

計算は1/4面だけやって後は対象性を利用する。

cal[] := 
      Module[{theta}, 
lines = {Dashing[{0.01, 
                0.01}]}; 
For[zst = 1, zst > 1/10, 
          zst -= 2/10, 
Print[zst]; 
x = 
            0; 
z = 
            zst; 
ts1 = {{x, 
                z}}; 
ts2 = {{-x, 
                z}}; 
ts3 = {{x, -z}}; 
ts4 = {{-x, -z}}; 
i = 
            0; 
While[
            Not[limit[x, z]] && i++ < 200, 
theta = 
              caltheta[x, z]; 

x = 
              x + dr*Cos[theta]; 
z = 
              z + dr*Sin[theta]; 
ts1 = 
              Append[ts1, {x, z}]; 
ts2 = 
              Append[ts2, {-x, z}]; 
ts3 = 
              Append[ts3, {x, -z}]; 
ts4 = 
              Append[ts4, {-x, -z}];
]; 
lines = 
            Append[lines, {Line[ts1], Line[ts2], Line[ts3], 
                Line[ts4]}];
]; 
Show[
          Graphics[lines], AspectRatio -> 1, 
          PlotRange -> {{-1.1, 1.1}, {-1.1, 
                1.1}}];
]

gra1 = cal[]

電界の計算

Er[r_, th_] := 
      Q1 a Cos[th]/(2 Pi ep0 r^3);    
Eth[r_, th_] := 
      Q1 a Sin[th]/(4 Pi ep0 r^3);  
Ex[x_, z_] := 
        Er[r[x, z], th[x, z]]*Sin[th[x, z]] + 
          Eth[r[x, z], th[x, z]]*Cos[th[x, z]];
Ez[x_, z_] := 
        Er[r[x, z], th[x, z]]*Cos[th[x, z]] - 
          Eth[r[x, z], th[x, z]]*Sin[th[x, z]];

電界の次の方向を求める。電界ベクトルから方向を決めてdrだけ進んだところのx,zの位置を計算する。

\(\(calx[x_, z_] := 
        Module[{length, ret = 0}, 
length = 
            Sqrt[Ex[x, z]^2 + Ez[x, z]^2]; 
	
          ret = x + Ex[x, z]/length*dr; 
Return[
            ret];
];\)\)
\(calz[x_, z_] := 
      Module[{length, ret = 0}, 
length = 
          Sqrt[Ex[x, z]^2 + Ez[x, z]^2]; 
	
        ret = z + Ez[x, z]/length*dr; 
Return[
          ret];
]\)

同様に1/4面のみ計算。スタートは点電荷の周りから行う。

cal2[] := 
      Module[{}, 
lines2 = {}; 
For[
          theta = 0/180*Pi, theta < Pi/2, 
          theta += 15/180*Pi, 
x = 
            dr*Cos[theta]; 
z = 
            a + dr*Sin[theta]; 
ts1 = {{x, 
                z}}; 
ts2 = {{-x, 
                z}}; 
ts3 = {{x, -z}}; 
ts4 = {{-x, -z}}; 
While[
            Not[limit[x, z]], 
x = 
              calx[x, z]; 
z = 
              calz[x, z]; 
ts1 = 
              Append[ts1, {x, z}]; 
ts2 = 
              Append[ts2, {-x, z}]; 
ts3 = 
              Append[ts3, {x, -z}]; 
ts4 = 
              Append[ts4, {-x, -z}];
]; 
lines2 = 
            Append[lines2, {Line[ts1], Line[ts2], Line[ts3], 
                Line[ts4]}];
]; 
Show[
          Graphics[lines2], AspectRatio -> 1, 
          PlotRange -> {{-1, 1}, {-1, 
                1}}];

]

dr = a/3; gra2 = cal2[]

Show[Graphics[{lines, lines2}], AspectRatio -> 1, 
        PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}];

この図と次のただの正負等量の点電荷で作る「電気双極子もどき」とをよく比較して欲しい。