<jupyter_start><jupyter_text>Test Case 4 Calculating interaction of randomly generated sphecial gaussian charge distribution with similar r_ext in 3D space<jupyter_code>import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import numpy as np
from fast_multipole_method import operation as op 
from scipy.special import erf
from scipy.special import erfc
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
plt.style.use('ggplot')
def plot_3d(x):
    """plot particles in 3 dimentional"""
    y = np.transpose(x)
    fig = plt.figure(figsize=(8,8))
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(y[0], y[1], y[2])
    ax = ax.view_init(30)
    plt.show()
    return 
# btm_level == 3, l = 1/8, for WS_index <=2
a_min = 2 * ((erfc(1-1e-16) * 8) ** 2)
a_min
#case 4.1 construction: random sphecial distributions, similar extend to make WS<=2
num_distribution = 100
x_i = np.ndarray(shape=(3, num_distribution))
x_i[0] = 20 * np.random.rand(num_distribution) - 10
x_i[1] = 20 * np.random.rand(num_distribution) - 10
x_i[2] = 20 * np.random.rand(num_distribution) - 10
x_i = np.transpose(x_i)
K_i = np.ones(num_distribution)
a_i = 10 * np.random.rand(num_distribution) + a_min

#case 4.2 constuction: uniform distributed sphecial distributions, same WS index
num_distribution_in_a_box = 1
num_distribution_1D = 3
num_distribution = num_distribution_in_a_box * num_distribution_1D **3
x_i = np.zeros(shape=(num_distribution,3))

for i in range(0,num_distribution_1D):
    for j in range(0,num_distribution_1D):
        for k in range(0,num_distribution_1D):
            x_i[i*num_distribution_1D*num_distribution_1D+j*num_distribution_1D+k] = [i,j,k]

K_i = np.ones(num_distribution)
a_i = 10 * np.random.rand(num_distribution) + a_min
[x0_i, scale_factor] = op.cartesian_scaling_to_unit_range(x_i)
plot_3d(x0_i)
a_i
# analytical answer
pair_potential = np.zeros(shape=(num_distribution,num_distribution)) 
pre_factor = np.power(np.pi, 3) 
for i in range(0, num_distribution):
    for j in range(i+1, num_distribution):
        pre_factor2 = K_i[i] * K_i[j] / ( np.power(a_i[i]*a_i[j], 1.5) * op.distance_cal(x0_i[i], x0_i[j]))
        t_sqrt = np.sqrt(a_i[i]*a_i[j]/(a_i[i]+a_i[j])) * op.distance_cal(x0_i[i], x0_i[j]) * scale_factor[1]
        pair_potential[i][j] = pre_factor * pre_factor2 * erf(t_sqrt)
        
pair_potential /= scale_factor[1]
pair_potential
J_analytic = np.zeros(num_distribution)
for i in range(0, num_distribution):
    for j in range(0, num_distribution):
        if j<i:
            J_analytic[i] += pair_potential[j][i]
        if j>i:
            J_analytic[i] += pair_potential[i][j]

J_analytic
total_energy = 0.5 * sum(J_analytic)
total_energy
from fast_multipole_method import fmm
from fast_multipole_method import fmm_q_gaussain_distribution as fq
# build list of q_source
q_source = np.ndarray(shape=(len(x0_i)), dtype=fq) 
for i in range(0, len(x0_i)):
    q_source[i] = fq(x0_i[i], a_i[i], K_i[i])
btm_level = 3
p = 10
ws_index = 3
[J_far_field, J_near_field] = fmm(q_source, btm_level, p, scale_factor[1], ws_index)
J_far_field
J_near_field
J_total = J_far_field + J_near_field
J_total
total_energy = 0.5 * sum(J_total)
total_energy
J_error = np.abs(J_total-J_analytic) / J_analytic
J_error<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup