一次元静電粒子コードについて・・(つづき1)
' 電荷密度の計算
For I = 1 To IM - 1
RO(I) = -1
Next I
For IP = 1 To NPT
XNM = X(IP) / DX:IXP = Int(XNM)
B = XNM - IXP :If IXP = 0 Then IXP = IM - 1
IXP1 = IXP + 1 :If IXP1 = IM Then IXP1 = 1
RO(IXP) = RO(IXP) + PW * (1 - B):RO(IXP1) = RO(IXP1) + PW * B
Next IP
RO(IM) = RO(1)
For I = 1 To IM - 1
RO(I) = RO(I) * DX * DX
Next I
' 静電ポテンシャルの計算
For I = 2 To IM - 1
A(I, 1) = 1:A(I, 2) = -2:A(I, 3) = 1:A(I, 4) = RO(I)
Next I
A(1, 4) = 0 :A(1, 2) = 1:A(IM, 4) = 0:A(IM, 2) = 1:BO = -2
' -------ガウスの消去法-------
NN = 1:N = IM:NNP1 = NN + 1:NN2 = 2 * NNP1
For I = 1 To N - 1
For J = NN2 To NNP1 Step -1
A(I, J) = A(I, J) / A(I, NNP1)
Next J
IPNN = I + NN:If IPNN > N Then IPNN = N
For L = I + 1 To IPNN
For K = NNP1 + I - L + 1 To 2 * NN + 1 + I - L
A(L, K) = A(L, K) - A(I, K + L - I) * A(L, NNP1 + I - L)
Next K
A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L)
Next L
Next
A(N, NN2) = A(N, NN2) / A(N, NNP1)
For I = N To 2 Step -1
IMNN = I - NN:If IMNN < 1 Then IMNN = 1
For L = I - 1 To IMNN Step -1
A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L)
Next L
Next I
A(1, NN2) = A(1, NN2) / A(1, NNP1)
For I = 1 To IM
PHI(I) = A(I, 4)
Next I
つづく
関連URL:http://www.okweb.ne.jp/kotaeru.php3?q=448303