FFT of gaussian if radial domain
2 views (last 30 days)
Show older comments
Hi
I defind a gaussian with cartesian dimain as written below.
How do I make the same but with radial? without the angle (theta).
In other words I want to transfer instead of x,y to 'r' :
exp(-(r(i_index)^2 )/w0^2)
and keep everything.
% Constants and Parameters
c = 299792458; % Speed of light [m/s]
twidth = 500e-15; % Time window [s]
%twidth = 5e-12;
n = 2^10; % Number of points
lambda0 = 795e-9; % Center wavelength [m]
W0 = 2*pi*c/lambda0; % Center frequency [rad/s]
lambda0_um = 1e6*lambda0;
FWHM = 34e-15; % Initial pulse duration [s]
Epsilon_0 = 8.8541878128*1e-12;
% Time domain
T = linspace(-twidth/2, twidth/2, n); % Time grid
dt = T(2) - T(1);
Dnu = 1/dt;
delta_T = twidth/n; % Time interval [s]
% Gaussian pulse in time domain
y = exp(1i*W0*T); % Carrier
sigma = FWHM/((2*log(2))^0.5);
gaussian = exp(-T.^2/sigma^2);
N = sqrt(sum(gaussian.^2)*dt);
gaussian = gaussian/N;
fin_t = y.*gaussian; % Electric field in time domain
% Spatial domain
x = linspace(-10e-3, 10e-3, 256); % at the focus, beam will be ~300um, so step size should be there approx 30um, not including self focusing
dx = x(2) - x(1);
w0 = 13e-3/2;
gaussian_space = zeros(length(x),length(x));
for i_index = 1:length(x)
for j =1:length(x)
gaussian_space(i_index,j) = exp(-(x(i_index).^2 + x(j).^2)/w0^2);
end
end
N = sqrt(sum(sum(gaussian_space.^2))*(dx)^2);
gaussian_space = gaussian_space / N;
Average_power = 14; % watt
Laser_operation_freq = 5e3; % Laser operation freq
Pulse_energy = Average_power / Laser_operation_freq;
Total_initial_intensity = 0.5*c*Epsilon_0*(sum(sum(gaussian_space.^2))*(dx)^2) * (sum(gaussian.^2)*dt);
Normalization = sqrt(Total_initial_intensity/Pulse_energy);
% Electric field in spatial and time domain
E_x_y_t = (1/Normalization)*gaussian_space.*reshape(fin_t,1,1,[]);
% Create frequency grids
frequencies = linspace(-Dnu/2, Dnu/2, n); % Frequency grid (time dimension)
n_x = length(x);
dfx = 1/(x(end)-x(1));
dfy = 1/(x(end)-x(1));
fx = linspace(-dfx*(n_x/2), dfx*(n_x/2), n_x); % Frequency grid (x dimension)
fy = linspace(-dfy*(n_x/2), dfy*(n_x/2), n_x); % Frequency grid (y dimension)
wavelength_freq_um = 1e6*c./frequencies; %wavelength in um
%% First Lens
first_focal_length = 4; %focal length in [m]
[i_index, j, k] = ndgrid(1:length(x), 1:length(x), 1:length(frequencies));
phase_first_focal_length = (2.*pi*frequencies(k)./c).*(x(i_index).^2+x(j).^2)/(2*first_focal_length);
%% now let's add the lens phase and return to the time domain
E_f = fftshift(fft(E_x_y_t, [], 3),3); % Frequency domain signal
E_f = E_f.*exp(-1i*phase_first_focal_length);
E_x_y_t = ifft(ifftshift(E_f,3), [], 3);
later i also add a different phase which correspond to the fourier of x,y:
phase_shift = sqrt(k_w_air(k).^2 - (2*pi*fx(i_index)).^2 - (2*pi*fy(j)).^2) * L_first_air_path;
delta_f = frequencies(k) - c/lambda0;
E_fx_fy_f_shifted = E_fx_fy_f_shifted .* exp(1i * phase_shift) .* exp(-1i*2*pi*delta_f *L_first_air_path/Vg_air);
Where fx , fy are the corresponding fourier of x, y so idealy should change to fr
0 Comments
Answers (1)
Balaji
on 8 Sep 2023
Hi Elis,
I understand that you want to perform ‘fft’ in radial domain.
Firstly you can change the Electric Field to radial domain as follows:
r = 0:dr:r_max;
n_r = length(r);
frequencies = linspace(-Dnu/2, Dnu/2, n);
% Electric field in radial and time domain
E_r_t = zeros(n_r, n);
for i = 1:n_r
radial_profile = exp(-r(i)^2/w0^2);
E_r_t(i, :) = radial_profile * fin_t;
end
Secondly for implementing fft in radial domain please refer to the following MATLAB answer:
Hope this helps!
Thanks
Balaji
0 Comments
See Also
Categories
Find more on Parametric Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!