Bipole Current Propagation Direction

[XueJingru]

As far as I know, the function "empymod.bipole" can work as a current line source whose result is as a intergration of "nums= scripts" dipoles' results. Does the line current have a propagation direction or can I adjust the propogation by modify the code or any parameters. My code is below and I wonder whether the part of sources above the surface make sense for final results.

import time

import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import mu_0
from scipy.special import erf

import empymod

stime = time.time()

def calculate_segment_length(end_src):
    # Extract coordinates for line 1
    line1_x1, line2_x1 = end_src[0]
    line1_x2, line2_x2 = end_src[1]
    line1_y1, line2_y1 = end_src[2]
    line1_y2, line2_y2 = end_src[3]
    line1_z1, line2_z1 = end_src[4]
    line1_z2, line2_z2 = end_src[5]

    # Calculate the length of line 1
    segment1_length = np.sqrt(
        (line1_x2 - line1_x1) ** 2
        + (line1_y2 - line1_y1) ** 2
        + (line1_z2 - line1_z1) ** 2
    )

    # Calculate the length of line 2
    segment2_length = np.sqrt(
        (line2_x2 - line2_x1) ** 2
        + (line2_y2 - line2_y1) ** 2
        + (line2_z2 - line2_z1) ** 2
    )

    return segment1_length, segment2_length

def generate_line_segments(base_src, end_src, num_segments):
    # Calculate the increment for each coordinate
    x1_increment = (end_src[1][0] - base_src[1][0]) / num_segments
    x2_increment = (end_src[1][1] - base_src[1][1]) / num_segments
    y1_increment = (end_src[3][0] - base_src[3][0]) / num_segments
    y2_increment = (end_src[3][1] - base_src[3][1]) / num_segments
    z1_increment = (end_src[5][0] - base_src[5][0]) / num_segments
    z2_increment = (end_src[5][1] - base_src[5][1]) / num_segments

    # Generate the line segments
    sources = []
    lengths = []
    for i in range(num_segments + 1):
        src = [
            [base_src[0][0], base_src[0][1]],
            [base_src[1][0] + i * x1_increment, base_src[1][1] + i * x2_increment],
            [base_src[2][0] + i * y1_increment, base_src[2][1] + i * y2_increment],
            [base_src[3][0] + i * y1_increment, base_src[3][1] + i * y2_increment],
            [base_src[4][0], base_src[4][1]],
            [base_src[5][0] + i * z1_increment, base_src[5][1] + i * z2_increment],
        ]
        sources.append(src)

        # Calculate the length of the segment
        length = calculate_segment_length(src)
        lengths.append(length)

    for i in range(num_segments + 1):
        src = [
            [base_src[0][0], base_src[0][1]],
            [base_src[1][0] + i * x1_increment, base_src[1][1] + i * x2_increment],
            [-base_src[2][0] - i * y1_increment, -base_src[2][1] - i * y2_increment],
            [-base_src[3][0] - i * y1_increment, -base_src[3][1] - i * y2_increment],
            [base_src[4][0], base_src[4][1]],
            [base_src[5][0] + i * z1_increment, base_src[5][1] + i * z2_increment],
        ]
        sources.append(src)

        # Calculate the length of the segment
        length = calculate_segment_length(src)
        lengths.append(length)

    y_increment = (end_src[3][1] + end_src[3][1]) / (2 * num_segments)
    for i in range(2 * num_segments):
        src = [
            [end_src[0][0], end_src[0][1]],
            [end_src[1][0], end_src[1][1]],
            [end_src[2][0] - i * y_increment, end_src[2][1] - i * y_increment],
            [end_src[3][0] - i * y_increment, end_src[3][1] - i * y_increment],
            [end_src[4][0], end_src[4][1]],
            [end_src[5][0], end_src[5][1]],
        ]
        sources.append(src)

        # Calculate the length of the segment
        length = calculate_segment_length(src)
        lengths.append(length)

    return sources, lengths
base_src = [
    [0.5, -0.5],
    [0.45, -0.45],
    [0.5, 0.5],
    [0.5, 0.5],
    [0.1, 0.1],
    [-0.1, -0.1],
]

end_src = [
    [0.5, -0.5],
    [0.25, -0.25],
    [0.25, 0.25],
    [0.25, 0.25],
    [0.1, 0.1],
    [-1.00, -1.00],
]

num_segments = 5

sources, lengths = generate_line_segments(base_src, end_src, num_segments)

point_x = np.linspace(-2, 2, 21)
point_y = np.array([0])
point_z = np.linspace(-2, 2, 21)

def mkvc(arr):
    return np.reshape(arr, (-1,))

rec_x, rec_y, rec_z = np.meshgrid(point_x, point_y, point_z)
print(rec_x[0, :, 0], rec_y[:, 0, 0], rec_z[0, 0, :])
rec_loc = np.c_[mkvc(rec_x), mkvc(rec_y), mkvc(rec_z)]
np.savetxt("examples/time_domain/data_horn/rec.dat", rec_loc)

rec_azi = 0
rec_dip = 90
rec_list = [rec_x[0, :, 0], rec_y[:, 0, 0], rec_z[0, 0, :], rec_azi, rec_dip]

time_channel = np.logspace(-5, -2, 31)

np.savetxt("time.dat", time_channel)
depth = [0, np.infty]
# np.savetxt("res.dat", res)
res_high = 2e2
res_low = 1e2
res_list = [1e7, res_low]

parameters = {
    "src": base_src,
    "rec": [],
    "depth": depth,
    "res": res_list,
    "freqtime": time_channel,
    "srcpts": 31,
    "strength": 1,
    "signal": -1,
    "verb": 3,
    "mrec": False,
    "msrc": False,
}

import concurrent.futures

def calculate_bipole(rec_tuple):
    rec_x, rec_y, rec_z, rec_azi, rec_dip = rec_tuple
    parameters["rec"] = [rec_x, rec_y, rec_z, rec_azi, rec_dip]
    # print(parameters["rec"])
    bipole_result = empymod.bipole(**parameters) / 1000  # * mu_0
    rec_xyz = np.array([rec_x, rec_y, rec_z])
    rec_xyz = np.tile(rec_xyz, (31, 1))
    time_channel_reshaped = np.transpose(time_channel).reshape(-1, 1)
    bipole_result_sum = np.sum(bipole_result, axis=1)
    result_point = np.column_stack(
        (
            rec_xyz,
            time_channel_reshaped,
            bipole_result,
            bipole_result_sum,
        )
    )
    return result_point

i = 0
j = 0

with concurrent.futures.ThreadPoolExecutor(max_workers=18) as executor:
    for s in range(len(sources)):
        # 提交任务给线程池
        parameters["src"] = sources[s]
        max_length = max(max(length) for length in lengths)
        parameters["strength"] = max_length / lengths[s][0]
        futures = [
            executor.submit(
                calculate_bipole,
                (
                    rec_list[0][i],
                    rec_list[1],
                    rec_list[2][j],
                    rec_list[3],
                    rec_list[4],
                ),
            )
            for i in range(len(rec_list[0]))
            for j in range(len(rec_list[2]))
        ]

        results = [
            future.result() for future in concurrent.futures.as_completed(futures)
        ]
        merged_result = np.vstack([result for result in results])
        np.savetxt(
            f"examples/time_domain/data_horn/artical_vertical_trapezoid_Ez-_source{s+1}.dat",
            merged_result,
        )

etime = time.time()
print("正演计算时间为:", format(etime - stime, ".2f"), "s")

data_list = []

for i in range(1, len(sources) + 1):
    filename = (
        f"examples/time_domain/data_horn/artical_vertical_trapezoid_Ez-_source{i}.dat"
    )
    data = np.loadtxt(filename)
    data_list.append(data)
print(np.shape(data_list))

sorted_data_list = []
for data in data_list:
    sort_indices = np.lexsort((data[:, 2], data[:, 0]))
    sorted_data = data[sort_indices]
    sorted_data_list.append(sorted_data)

max_rows = max([data.shape[0] for data in sorted_data_list])

sum_last_column = np.zeros(max_rows)
for sorted_data in sorted_data_list:
    last_column = sorted_data[:, -1]
    sum_last_column[: last_column.shape[0]] += last_column

total_result = np.hstack((sorted_data[:, :3], sum_last_column.reshape(-1, 1)))

print(total_result.shape)
split_arrays = [total_result[i :: len(time_channel)] for i in range(len(time_channel))]
for i, arr in enumerate(split_arrays):
    np.savetxt(
        f"examples/time_domain/data_horn/chosen_time/artical_vertical_trapezoid_Ez-_source_time_channel_{i+1}.dat",
        arr,
    )

[prisae]

Hi @XueJingru - apologies, but I am lacking the time right now to look in detail into your long script. I do have the impression though that what you are trying to do could be achieved much simpler.

Could you maybe make a drawing or a sketch of the layout you are trying to model? That would help to understand your setup.

[XueJingru]

The sources I want to calculate is shown in the png. I just want to know whether these electric line sources have current direction or the result of empymod.bipole function is the just the intergration of the result of empymod.dipole along the line.

[prisae]

That looks super interesting! Although I fear it might be difficult to model with empymod, but not impossible. But you have to take care that the transforms do not incur any artefacts. Particularly, you use very early times (1e-5 to 1e-2 s), which requires high frequencies for the transforms. You might need a bandpass filter for it to properly work, see, e.g., https://empymod.emsig.xyz/en/stable/gallery/tdomain/tem_walktem.html

empymod.bipole is simply doing a Gaussian Quadrature over empymod.dipole.

However, even a dipole has a direction, and so das bipole. You can test that with the following script:

import empymod

model = {'rec': [1000, 0, 0, 0, 0], 'depth': [0, 100], 'res': [1, 100, 3], 'freqtime': 1, 'verb': 1}

bip1 = empymod.bipole(src=[0, 100, 0, 0, 0, 0], **model)
bip2 = empymod.bipole(src=[100, 0, 0, 0, 0, 0], **model)

print(bip1, bip2)

The result is

(2.0351810723580752e-10-4.254477006114425e-11j) (-2.0351810723580752e-10+4.254477006114425e-11j)

so you see the result differ by a factor -1, which comes from the direction.

[XueJingru]

Thanks very much for your guidance during your busy time.

[prisae]

I hope it helps. I close this issue, if there are still open points you would like to get help feel free to re-open it.