Matlab code EMP system for anti UAV with range exceeding 5km but I am unable to get the intended graphs, I dont know what wrong I am doing

5 次查看（过去 30 天）
Taimour
Taimour 2024-12-31
回答： UDAYA PEDDIRAJU 2025-4-21
Ran in:
%%
% source MATLAB
%Pulse propagation:
%
% /\ /\
% E(x,t) = | dw | dx {E(kx,w) * exp[i*(w*t - kx*x - kz*z)]}
% \/ \/
%
% / / \ \
% = FT < FT | E(kx,w), (kx -> x) |, (w -> t-|k|z) >
% \ \ / /
% kz = sqrt(|k|^2 - kx^2) - |k|
% |k| = w/c
%
% Note this is in the time frame moving along the propagation axis.
%
% Angular dispersion is when the propagation direction is dependant on frequency:
% k -> k - alpha * w
%
% This is equivalent in the spatial domain to a linear phase:
% E(x, w) -> |E(x, w)| * exp(i * alpha * w * x)
%
% This is then equivalent to a spatially dependant shift in the time domain:
% E(x, t) -> E(x, t - alpha * x)
%
% This thus represents pulse front tilt.
% As the pulse propagates along axis (z direction), it will diffract. Thus the dominant
% phase term is the radius of curvature (almost quadratic) sqrt(|k| - kx^2)*z. Thus any
% STC is lost, except that the pulses is delayed the further from the propagation axis.
% Some useful functions
gauss = @(x, x0, FWHM) 2.^(-(2*(x-x0)./FWHM).^2);% Gaussian of full width half max (FWHM) centred at x0
chng_rng = @(x, c) 2*pi*c./x; % Changed between wavelength and angular frequency
% Initialise array sizes
Nx = 256; % # spatial points
Nk = Nx; % # k-vector points (FT of x)
Nw = 256; % # freq. points
Nt = Nw; % # temporal points (FT of w)
c = .3; % Speed of light um/fs
% Setup frequency data
wmin = chng_rng(1, c); % Min w value
wmax = chng_rng(.5, c); % Max w value
w0 = mean([wmin, wmax]); % Central wavelength (rad/fs)
wFWHM = w0*.01; % Bandwidth
w = (0:Nw-1)*(wmax-wmin)/(Nw-1) + wmin; % Freq. range rad/fs
t = ((0:Nt-1)-Nt/2)*2*pi/(wmax-wmin); % Time domain /fs
W = ones(Nx, 1)*w; % Freq. range for all kx
EW = gauss(W, w0, wFWHM); % Spectrum
% Setup spatial data & angular dispersion
x = ((2*(0:Nx-1)'/(Nx-1))-1)*1e4; % Spatial co-ord /um
k = fftshift(((0:Nk-1)'-Nk/2)*2*pi/(max(x)-min(x)));
% kx-values
K = k*ones(1, Nw); % kx values for all w.
alpha = 0.5; % Angular dispersion
xFWHM = 500; % Beam width /um
EX = (gauss(x, 0, xFWHM)*ones(1,Nw)) ... % Spatial profile
.* exp(i*alpha*x*(2*(0:Nw-1)/(Nw-1)-1)); % Angular dispersion
EK = abs(ifft(ifftshift(EX, 1), [], 1)); % k-vector components
% Propagation phase & Electric field
phi_kw = @(z) (sqrt((W/c).^2-K.^2)-W/c)*z; % Phase along z-axis
E_kw = @(z) EK.*EW.*exp(-i*phi_kw(z));
% Electric field for each k-w plane wave
% Initial Pulse
figure()
imagesc(t, x*1e-3, abs(fftshift(fft2(E_kw(0)))).^2);
axis([-50 50 -1 1]);
grid on;
% Pulse Propagation
delta = 5; % Used for plotting above & below axis profile
for z=0:.2e5:1e6
I_xt = abs(fftshift(fft2(E_kw(z)))).^2; % Pulse intensity in space & time at pos z
subplot(2,1,1)
imagesc(t, x*1e-3, I_xt);
axis([-350 350 -2 2])
grid on;
title(['z = ' num2str(z*1e-6, '%.3f') 'm']);
xlabel('Time /fs');
ylabel('Beam profile /mm');
subplot(2,1,2)
plot(t, I_xt([Nx/2-delta+1 Nx/2 Nx/2+delta+1], :));
title('Pulse profile along propagation axis');
xlabel('Time /fs');
ylabel('Intensity /arb.');
legend('Above axis', 'On axis', 'Below axis');
axis([-350 350 0 max(max(I_xt([Nx/2-delta Nx/2 Nx/2+delta], :), [], 2))]);
getframe;
end
%%
clear;
clc;
close all;
% Exponential Decay
% Define the exponential decay function
exponential_decay = @(y0, lambd, t) y0 * exp(-lambd * t);
% Parameters for Exponential Decay
y0 = 100; % Initial value
lambd = 0.1; % Decay constant
t_decay = linspace(0, 50, 500); % Time range from 0 to 50
% Compute the decay
y_decay = exponential_decay(y0, lambd, t_decay);
% EMP Pulse Propagation
% Simulation Parameters
c = 3e8; % Speed of light (m/s)
frequency = 1e9; % EMP pulse frequency (1 GHz)
pulse_duration = 1e-9; % Pulse duration (1 ns)
E_peak = 1e6; % Peak electric field strength (V/m)
domain_size = 100; % Domain size (meters)
grid_points = 200; % Number of spatial grid points
time_steps = 300; % Number of time steps
dx = domain_size / grid_points; % Spatial step (meters)
dt = dx / (2 * c); % Time step based on Courant condition
% UAV Target Parameters
target_position = 75; % Target UAV position (meters)
E_disruption = 5e5; % Electric field required to disrupt UAV electronics (V/m)
% Initialize Fields
E = zeros(1, grid_points); % Electric field
H = zeros(1, grid_points); % Magnetic field
% Gaussian Pulse Definition
t0 = pulse_duration * 3; % Pulse center time
pulse = @(t) E_peak * exp(-((t - t0) / (pulse_duration / 2))^2);
% Absorbing Boundary Conditions
attenuation = exp(-0.01 * (1:grid_points)); % Simple exponential attenuation at boundaries
% Set up figure for combined visualization
figure;
% Plot the Exponential Decay Curve
subplot(2, 1, 1);
plot(t_decay, y_decay, 'r', 'LineWidth', 1.5);
title('Exponential Decay');
xlabel('Time');
ylabel('y(t)');
grid on;
% Initialize EMP Visualization
subplot(2, 1, 2);
hold on;
% Time-Stepping Loop for EMP Pulse Propagation
for n = 1:time_steps
% Update Magnetic Field
H(1:end-1) = H(1:end-1) - (dt / dx) * diff(E);
% Update Electric Field
E(2:end) = E(2:end) - (dt / dx) * diff(H);
% Apply Absorbing Boundary Conditions
E = E .* attenuation;
% Inject Gaussian Pulse at Source
E(1) = pulse(n * dt);
% Plot EMP Field
subplot(2, 1, 2);
plot((1:grid_points) * dx, E, 'b', 'LineWidth', 1.5);
yline(E_disruption, 'r--', 'Disruption Threshold', 'LineWidth', 1.2);
scatter(target_position, 0, 100, 'k', 'filled', 'DisplayName', 'Target UAV');
title(['EMP Pulse Propagation - Time Step: ', num2str(n)]);
xlabel('Position (m)');
ylabel('Electric Field (V/m)');
grid on;
axis([0 domain_size -E_peak E_peak]);
pause(0.01);
% cla;
% Check if the UAV is disrupted
if abs(E(round(target_position / dx))) >= E_disruption
disp(['UAV disrupted due to EMP at time step: ', num2str(n)]);
break;
end
end
●
if abs(E(round(target_position / dx))) < E_disruption
disp('Target UAV not disrupted within simulation range.');
end
Target UAV not disrupted within simulation range.
  1 个评论
Walter Roberson
Walter Roberson 2025-1-1
We do not know what the intended graphs are like.
Note: you loop outputting the EMP Pulse Propagation. You would be able to monitor the progression "live" but only the final graph shows up here on MATLAB Answers.,

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

UDAYA PEDDIRAJU
UDAYA PEDDIRAJU 2025-4-21
Hi Taimour,
The issue likely stems from either incorrect parameter settings, insufficient simulation range or resolution, or errors in the visualization commands within the MATLAB code. Since the EMP anti-UAV system involves complex wave propagation and angular dispersion modeling, ensure that:
  • The spatial and frequency domain parameters (e.g., Nx, Nw, spatial coordinates, frequency range) are set to cover the physical scale exceeding 5 km appropriately.
  • The propagation phase and electric field calculations correctly reflect the EMP pulse behavior over long distances.
  • The plotting commands correctly display the results, with axis limits and scaling adjusted to visualize the expected phenomena.
  • Any necessary data transformations (e.g., FFT shifts, scaling) are applied properly before plotting.
Additionally, MATLAB's UAV and Simulink tools provide built-in functionalities for UAV system modeling, including flight control, obstacle avoidance, and sensor simulation, but may not directly support EMP system modeling beyond certain ranges without custom extensions.
Refer to the following documentation for more details:
1) https://www.mathworks.com/help/uav/ug/uav-obstacle-avoidance-in-simulink.html
2) https://www.mathworks.com/solutions/robotics/uav.html
3) https://www.mathworks.com/help/uav/ug/uav-scenario-tutorial.html

类别

在 Help Center 和 File Exchange 中查找有关 UAV 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by