PIC法によるプラズマシミュレーションコードのための外挿ができたがここから内挿をどうやってするかわからない


わからない点  
cahrgeのデータ保存場所の名前。外挿した後のcharge(`HUT`,`HST`,`MST`,`MUT`)名前はそれぞれ変えたほうがいいのか。各所でオブジェクトを表示する処理を行ったがないほうがいいのか。各格子点に割り振られた電荷を用いて電位計算を行う。計算した電位により発生した電界を内挿により荷電粒子に加わる静電気力として与える。  
下記に現時点のプログラムを記載します。

[code]

    from matplotlib import pyplot as plt
    from matplotlib import cm
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    import matplotlib.pyplot as plt
    import itertools
    import csv
    import threading
    import random
    import time
    from numba import jit  # Better to write at the beginning of the code
    import decimal
    import math
    
    
    
    # 条件設定
    delta_L = 1.0
    delta = 1.0
    n_inter = 0
    nx = 1000         #x方向のメッシュの数
    ny = 1000       #y方向のメッシュの数
    xmin = 0         #xの最小値
    xmax = 1000*delta #xの最大値
    ymin = 0         #yの最小値
    ymax = 1000*delta #yの最大値
    n0 = 1e12;       #density in [/m^3]
    phi0 = 0;        #reference potential
    Te = 1;          #electron temperature in eV
    v_drift = 7000;  #ion injection velocity, 7[km/s]
    phi_p = -7;      #plate potential
    
    
    
    convergence_criterion = 1.0e-5   #前回との差
    
    
    
    dx = (xmax - xmin) / (nx - 1)
    dy = (ymax - ymin) / (ny - 1)
    
    # 初期状態
    phi  = np.zeros((ny, nx))
    electrode  = np.zeros((ny, nx))
    x  = np.linspace(xmin, xmax, nx)
    y  = np.linspace(xmin, xmax, ny)
    
    
    
    # for SOR method
    aa_recta=0.5*(np.cos(np.pi/nx)+np.cos(np.pi/ny)) #
    omega_SOR_recta = 2/(1+np.sqrt(1-aa_recta**2)) # 長方形領域に対する最適加速パラメータ
    print("omega_SOR_rect=",omega_SOR_recta)
    
    
    
    T=1 #粒子計算における時間ステップ
    
    # 粒子データを格納するための辞書を初期化
    particle_data = {}
    
    # 粒子数
    n_particles = 100
    
    # 粒子の座標と名前をランダムに生成して粒子データに追加
    for i in range(n_particles):
        a = (random.randint(0, 1000000)) / 1000
        b = (random.randint(0, 1000000)) / 1000
        particle_name = f"Particle_{i}"  # 粒子の名前を設定
        particle_data[particle_name] = {"a": a, "b": b, "charge": []}  # chargeキーを追加し初期化
    
    # 粒子ごとに計算と結果の保存を行う
    for name, data in particle_data.items():
        a = data["a"]
        b = data["b"]
    
        # 粒子ごとの計算
        aa = math.floor(a)
        bb = math.floor(b)
        
        a = math.floor(a * 1000)
        b = math.floor(b * 1000)
        
        aa = aa * 1000
        bb = bb * 1000
        
        a1 = (a - aa)
        a2 = ((aa + 1000) - a)
        b1 = (b - bb) * 1000
        b2 = ((bb + 1000) - b)
        
        HU = b2 * a1
        MU = a2 * b2
        HS = a1 * b1
        MS = a2 * b2
        MEN = HU + HS + MU + MS
        
        HUT = (HU / MEN * T)
        HST = (HS / MEN * T)
        MUT = (MU / MEN * T)
        MST = (HS / MEN * T)
        
        aa = math.floor(aa / 1000)
        bb = math.floor(bb / 1000)
        
        charge = np.zeros((ny, nx))
        
        charge[aa, bb+1] = HUT
        charge[aa, bb] = HST
        charge[aa+1, bb] = MST
        charge[aa+1, bb+1] = MUT
        
        # 粒子の計算結果をデータに保存
        data["charge"] = charge
        
        # 計算結果をファイルに保存
        with open("{}.csv".format(name), "w") as f:
            writer = csv.writer(f)
            writer.writerows(charge)
        
    # 計算結果の出力（任意の方法で）
    for name, data in particle_data.items():
        charge = data["charge"]
        print("Particle Name:", name)
        print("Charge Matrix:")
        print(charge)
        
        
    # 電荷置
    eps0 = 1
    
    
    # 電位計算　外挿？
    @jit(nopython=True)
    def calc_checkered(offset, ph):
        for i in range(nx):
            for j in range((i + offset) % 2, ny, 2):
                if i == 0 or i == nx - 1 or j == 0 or j == ny - 1:
                    ph[i, j] = 0
                else:
                    ph[i, j] = ph[i, j] + omega_SOR_recta * (
                        (ph[i+1, j] + ph[i-1, j] + ph[i, j+1] 
                            + ph[i, j-1])/4- ph[i, j] 
                            - (dx*dy/(4*eps0))*charge[i, j])
    
    
    
    conv_check = []
    while delta > convergence_criterion:
        phi_in = phi.copy()
        if n_inter % 50 == 0:
            print("interaction No =", n_inter, "delta=", delta)
            conv_check.append([n_inter, delta])
        calc_checkered(0, phi)  # red group
        calc_checkered(1, phi)  # black group
        delta = np.max(abs(phi - phi_in))
        n_inter += 1
            #red-black-SOR法
    
    
    # 電場計算　内挿？
    
    
    
    # ファイルを書き込みモードでオープン
    body = (phi)
    header = ('Electrical potential of 2D array')
    with open('Deni1.csv', 'w') as f:
    
        writer = csv.writer(f)  # writerオブジェクトを作成
        writer.writerows(header) # ヘッダーを書き込む
        writer.writerows(body)  # 内容を書き込む
    
    
    
    body2 = (x)
    header2 = ('Value1 of axis')
    with open('jikunoatai1.csv', 'w') as f1:
        writer = csv.writer(f1)
        writer.writerow(header2)
        writer.writerow(body2)
    
    
    
    def plot2D(x, y, phi):
        fig = plt.figure(figsize=(11, 7), dpi=100)
        ax = fig.add_subplot(111, projection='3d')
        X, Y = np.meshgrid (x , y)
        surf = ax.plot_surface(X, Y, phi[:], rstride=1, cstride=1, 
                cmap=cm.viridis,linewidth=0, antialiased=False)
        ax.view_init(30, 225)
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_zlabel('phi')
    
    plot2D(x, y, phi)
    plt.show()
    
    
    
    #オブジェクトを生成
    fig = plt.figure(figsize=(11,7), dpi=100)
    fig
    
    xmin = 0
    xmax = 1000
    ymin = 0
    ymax = 1000
    
    x  = np.linspace(xmin, xmax, nx)
    y  = np.linspace(xmin, xmax, ny)
    
    X, Y = np.meshgrid(x, y)
    plt.contourf(X, Y, phi, alpha=0.5, cmap=cm.viridis)
    plt.colorbar()
    plt.xlabel('X')
    plt.ylabel('Y')
    
    
    
    #電位計算
    print("Keisan V")
    for i in range(1,6):
        print("x=",i*delta_L,"V=",phi[500+i,500] 
                - phi[500+int(1/delta_L),500])
    
    
    
    #電場計算
    L = 999
    Ey = np.zeros([L,L])
    for i in range(1, L-1):
        for j in range(1, L-1):
            Ey[i,j] = -(phi[i+1,j]-phi[i-1,j])/2/delta_L
    print("Electric Field Calculation")
    print("Keisan E")
    for i in range(1,6):
        print("x=",i*delta_L,"E=",Ey[500+i,500])
    
    # 電場の可視化
    fig = plt.figure(figsize=(11, 7), dpi=100)
    ax = fig.add_subplot(111, projection='3d')
    
    X, Y = np.meshgrid(x[1:L-1], y[1:L-1])
    surf = ax.plot_surface(X, Y, Ey[1:L-1, 1:L-1], rstride=1, cstride=1, 
    cmap=cm.viridis, linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('Electric Field (Ey)')
    
    plt.show()
    
[/code]

