Hello Djaborebbi,
Let’s consider a basic approach where the lightning channel is divided into small dipoles, and the contribution from each dipole to the electromagnetic field is calculated and summed. You can modify this according to your requirements
The code below follows this. Here I have represented the return stroke current with the variable “I”
n = 100; % Number of segments
dz = H / n; % Length of each segment
z = linspace(0, H, n); % Positions of the segments along the channel
I0 = 1; % Peak current at the channel base (amplitude of Heidler function)
v = 1.5e8; % Return stroke speed (m/s), adjust as needed
epsilon0 = 8.854e-12; % Permittivity of free space (F/m)
mu0 = 4 * pi * 1e-7; % Permeability of free space (H/m)
% Observation point
r = 5000; % Radial distance from the channel (m)
theta = pi / 2; % Observation angle (horizontal plane)
% Initialize fields
E_r = 0;
E_theta = 0;
H_phi = 0;
% Calculate the electromagnetic field contributions from each dipole
for k = 1:n
% Position of the k-th dipole
zk = z(k);
% Current at the k-th dipole
Ik = I(k);
% Distance from dipole to observation point (approximation)
R = sqrt(r^2 + zk^2);
% Calculate the contribution to the fields
dE_r = (mu0 / (4 * pi)) * (Ik * dz * zk / R^3) * (3 * cos(theta)^2 - 1);
dE_theta = (mu0 / (4 * pi)) * (Ik * dz * r / R^3) * (3 * cos(theta) * sin(theta));
dH_phi = (Ik * dz * sin(theta)) / (4 * pi * R^2);
% Accumulate the contributions
E_r = E_r + dE_r;
E_theta = E_theta + dE_theta;
H_phi = H_phi + dH_phi;
end
% Display results
fprintf('Electric Field (Radial Component): %e V/m\n', E_r);
fprintf('Electric Field (Theta Component): %e V/m\n', E_theta);
fprintf('Magnetic Field (Phi Component): %e A/m\n', H_phi);
