電気力線描画で余分な線の原因

 投稿者：島村1243  投稿日：2016年 3月15日(火)20時22分52秒
  一直線状に電荷が均等分布する導体が、空間中に複数本、平行に置かれた場合の静電場の電気力線方程式（厳密解）は
  （１） x軸-y軸平面状の任意点(x,y)
　（２） Q(i)：第i番目導体上の電荷密度
  （３）(Xq(i),Yq(i))は第i番目導体の位置座標
　（４）θ(i)＝ANGLE(x-Xq(i),y-Yq(i))
とすると
　（５）∑{Q(i)*θ(i)}＝Const
で表されます。

そこで、すべての(x,y)点の値Aq(x,y)を
　　Aq((x,y)=∑{Q(i)*θ(i)}
       i=１～導体総数
で求めておき、このAq((x,y)が等しい値C(パラメータ）を持つx-y座標点を連ねると、電気力線を描く事ができます。

点を連ねる判断は
IF Aq(x,y)<>Aq(x,y+1) OR Aq(x,y)>Aq(x+1,y) OR Aq(x,y)<>Aq(x+1,y+1) THEN
      PLOT POINTS:x,y
で行っています。
このIF文の考え方は
　Aq(x,y)=Aq(x,y+1) =Aq(x+1,y)=Aq(x+1,y+1)
であるなら、点を連ねる線は一回りしてしまうので不合理、つまり、一回りしないことを
　Not｛Aq(x,y)=Aq(x,y+1) =Aq(x+1,y)=Aq(x+1,y+1)｝＝TRUE
でチェックしたのが上記のIF文で、OKなら描画面の(x,y)点を打点する、というものです。

この考え方を元に、３本の直線導体が有る場合の電気力線作図プログラムが後期するコードですが、RUNすると、最上・下の導体に、「明らかに余分と分かる水平線が２本描かれ」てしまい、その原因がどうにも掴めません。

誤りの原因がお分かりになりましたらご教示頂けると有り難いです。

*****トラぶっているプログラム***********
!****指定値(整数)セット（その1/2）***
LET w=500       !描画画面の横幅
LET h=500       !描画画面の高さ
LET Nq=3        !扱う総電線数
LET Esen=20 !電気力線の描画数[本/C]
LET Vstep=6 !等電位線の描画本数係数
!**************************

DIM Aq(w,h),Vq(w,h)
DIM Q(6),Xq(6),Yq(6)

!****指定値セット（その2/2）****
LET Q(1)=1     !第1番目電線に分布する電荷密度[C/m]
LET Xq(1)=250  !第1番目電線のx座標（1以上でw以下）
LET Yq(1)=350  !第1番目電線のy座標（1以上でh以下）

LET Q(2)=-2
LET Xq(2)=250
LET Yq(2)=250

LET Q(3)=1
LET Xq(3)=250
LET Yq(3)=150
!***************************

!**** Main ProgramStart ******
SET WINDOW 0,w,0,h
Draw Grid (50,50)
SET POINT STYLE 1
LET hankei=10

LET Nriki=Esen/2
LET hankei2=hankei*hankei
CALL Marks
PRINT "作図データの計算中です。"
CALL PointData
PRINT "電気力線を作図します。"
Call Rikisen
PRINT "作図完了です。"

!****各関数の内容**********
Sub Marks !全電荷の位置をマークする
   FOR id=1 TO Nq
      LET x1=Xq(id)
      LET y1=Yq(id)
      LET q1=Q(id)
      IF Q(id)>0 THEN
         SET AREA COLOR "red"
      ELSEIF Q(id)<0 THEN
         SET AREA COLOR "blue"
      ELSE
         Exit sub
      END IF
      DRAW disk WITH SCALE(hankei)*SHIFT(Xq(id),Yq(id))
   NEXT id
End Sub

Sub PointData !作図データの作成
   FOR x=1 TO w
      FOR y=1 TO h
         LET Aq(x,y)=INT(Nriki*Qangle(x,y)/PI)
         !PRINT "x=";x;"y=";y;"Aq(x,y)=";Aq(x,y)!<===Test
      Next Y
   Next     X
End Sub

FUNCTION Volt(x,y) !<---(x,y)点の合成電位計算
   LET Vd=0.0
   FOR  iv=1 TO Nq    !Nqは総電線数
      LET x1=Xq(iv)   !Q(i)のX座標
      LET y1=Yq(iv)   !Q(i)のY座標
      LET dx=x-x1
      LET dy=y-y1
      IF dx<>0 OR dy<>0 THEN
         LET Vd =Vd-Q(iv)*LOG(dx*dx+dy*dy)
      END IF
   NEXT iv
   LET Volt=Vd
End Function

FUNCTION Qangle(x,y)  !(x,y)点の∑(Qi×θi)計算
   LET Sigma=0.0
   FOR iq=1 TO Nq
      LET dx=x-Xq(iq)
      LET dy=y-Yq(iq)
      IF ichiCheck(x,y)=0 THEN !<--ANGLE(0,0)エラー対策
         IF dx<0 AND dy=0 THEN
            LET angle1=-PI
         ELSE
            LET angle1=ANGLE(dx,dy)
         END IF
         LET Sigma=Sigma+Q(iq)*angle1
      END IF
   NEXT iq
   LET Qangle=Sigma
End Function

FUNCTION ichiCheck(x,y) !点(x,y)がどれかの電線内に入ったかの判定
   FOR ii=1 TO Nq
      LET xi=Xq(ii)
      LET yi=Yq(ii)
      IF (xi-x)*(xi-x)+(yi-y)*(yi-y) <= hankei2 THEN
         LET ichiCheck=1 !位置(x,y)が電線内と認識
         Exit Function
      END IF
   NEXT ii
   LET ichiCheck=0 !位置(x,y)が電線外と認識
End Function

Sub  Rikisen !電気力線を描く
   FOR x=1 TO w-1
      FOR y=1 TO h-1
         IF ichiCheck(x,y)=0 THEN
            IF Aq(x,y)<>Aq(x,y+1) OR Aq(x,y)>Aq(x+1,y) OR Aq(x,y)<>Aq(x+1,y+1) THEN
               PLOT POINTS:x,y
            END IF
         End IF
      NEXT Y
   NEXT X
END SUB

END 		 

Re: 電気力線描画で余分な線の原因

 投稿者：山中和義  投稿日：2016年 3月16日(水)12時55分59秒
  島村1243さんへのお返事です。

> 誤りの原因がお分かりになりましたらご教示頂けると有り難いです。


表示部分で不必要なデータを削除しました。
根本的なところは、「なぜ、このような構造のデータになるか」ですが、今のところ不明です。

Sub Rikisen !電気力線を描く
   FOR x=1 TO w-1
      FOR y=1 TO h-1
         IF ichiCheck(x,y)=0 THEN
            IF Aq(x,y)<>Aq(x,y+1) OR Aq(x,y)<>Aq(x+1,y) OR Aq(x,y)<>Aq(x+1,y+1) THEN
               IF NOT( Aq(x,y)=Aq(x+1,y) AND Aq(x,y+1)=Aq(x+1,y+1) AND MOD(ABS(Aq(x,y)-Aq(x,y+1)),Esen)=0 ) THEN
                  PLOT POINTS:x,y

                  !!PRINT x;y,Aq(x,y);Aq(x,y+1);Aq(x+1,y);Aq(x+1,y+1)
               END IF
            END IF
         End IF
      NEXT Y
   NEXT X
END SUB



 		 

Re: 電気力線描画で余分な線の原因

 投稿者：島村1243  投稿日：2016年 3月16日(水)18時25分35秒
  山中和義 様へ

> 表示部分で不必要なデータを削除しました。
> 根本的なところは、「なぜ、このような構造のデータになるか」ですが、今のところ不明です。

修正コードのご教示有難う御座いました。
こんなに速く解決策をご教示頂けるとは、と驚いています。
早速修正コードを適用し、計算条件を色々変えて試用しましたが、完璧です！！！
厳密解の線図を知りたいと思っていたので感激しています。
重ねて御礼申し上げます。
根本的な理由については、無い頭ですが推測してみます。


 		 

Re: 電気力線描画で余分な線の原因

 投稿者：山中和義  投稿日：2016年 3月17日(木)21時07分6秒
  島村1243さんへのお返事です。

> 根本的な理由については、無い頭ですが推測してみます。

angle関数が(-π,π]の範囲なので、Σの値が負になる可能性があります。
これにより、-πからπに移るとき、Esen分（その倍数）の段差が発生するので、ここに線が描かれます。
これが余計な水平な線になります。
したがって、負と0の場合、げたを履かせて補正するといいでしょう。

Sub PointData !作図データの作成
   FOR x=1 TO w
      FOR y=1 TO h
         LET Aq(x,y)=INT(Nriki*Qangle(x,y)/PI)
         !PRINT "x=";x;"y=";y;"Aq(x,y)=";Aq(x,y)!<===Test

         DO WHILE Aq(x,y)<Esen
            LET Aq(x,y)=Aq(x,y)+Esen
         LOOP
      NEXT Y
   NEXT X
End Sub


FUNCTION Qangle(x,y)  !(x,y)点の∑(Qi×θi)計算
   LET Sigma=0.0
   FOR iq=1 TO Nq
      IF ichiCheck(x,y)=0 THEN !<--ANGLE(0,0)エラー対策
         LET Sigma=Sigma+Q(iq)*ANGLE(x-Xq(iq),y-Yq(iq))
      END IF
   NEXT iq
   LET Qangle=Sigma
End Function


こちらは元に戻します。

Sub Rikisen !電気力線を描く
   FOR x=1 TO w-1
      FOR y=1 TO h-1
         IF ichiCheck(x,y)=0 THEN
            IF Aq(x,y)<>Aq(x,y+1) OR Aq(x,y)<>Aq(x+1,y) OR Aq(x,y)<>Aq(x+1,y+1) THEN
               PLOT POINTS: x,y
               !!!PRINT x;y,Aq(x,y);Aq(x,y+1);Aq(x+1,y);Aq(x+1,y+1) !debug
            END IF
         END IF
      NEXT y
   NEXT x
End Sub

 		 

Re: 電気力線描画で余分な線の原因

 投稿者：島村1243  投稿日：2016年 3月18日(金)15時50分42秒
  > No.4003[元記事へ]

山中和義様へのお返事です。

> angle関数が(-π,π]の範囲なので、Σの値が負になる可能性があります。
> これにより、-πからπに移るとき、Esen分（その倍数）の段差が発生するので、ここに線が描かれます。
> これが余計な水平な線になります。
> したがって、負と0の場合、げたを履かせて補正するといいでしょう。

数学の知識・プログラミングの知識が浅い私には、上記が原因とは全く思いも浮かびませんでした。
度重なる丁寧なご教示、有難う御座いました。

 		 

戻る