%**************************************************************************
%This function calculates the fourier transform of the electric field of
%the radiation from a charge in the undulator.
%
%Input:
%x, x position of the trajectory (m)
%y, y position of the trajectory (m)
%z, z position along the undulator (m)
%vx, x component of the velocity (m/s)
%vy, y component of the velocity (m/s)
%vz, z component of the velocity (m/s)
%nx, x component of unit vector in direction of interest
%ny, y component of unit vector in direction of interest
%nz, z component of unit vector in direction of interest
%omega, angular frequency of interest (1/s)
%
%Output:
%ft_er, fourier transform of the electric field times the distance to the observation point as a function of z (V/m 1/(rad/sec) m)
%
%Zachary Wolf
%7/9/10

function [ft_er] = rad_calc_ft_er(x, y, z, vx, vy, vz, nx, ny, nz, omega)

%Method
%See Hofmann, p. 36, eq. 2.38
%E(t) = (q / (4 pi eps_0 c)) [(n x ((n - beta) x beta_dot)) / (r (1 - n . beta)^3)]_ret
%B(t) = (1 / c) (n x E)
%E(omega) = (i omega q / (4 pi sqrt(2 pi) eps_0 c r)) Int((n x (n x beta)) exp(-i omega (t + r(t)/c)) dt)
%B(omega) = (1 / c) (n x E(omega))
%vz = dz / dt
%t = Int((1 / vz) dz)
%dt = (1 / vz) dz
%f(t) = f(z)

%Debug
%tic;

%Physical constants
global Q_ELEC;          %Electron charge (C)
global C_LIGHT;         %Speed of light (m/s)
global EPS_0;           %Permittivity of free space (F/m)

%Initialize
ft_er = 0;

%Check
if (length(z) ~= length(x)) || (length(z) ~= length(y))
    msgbox('rad_calc_ft_er: arrays have different lengths', 'Error');
    return;
end
if (length(z) ~= length(vx)) || (length(z) ~= length(vy)) || (length(z) ~= length(vz))
    msgbox('rad_calc_ft_er: arrays have different lengths', 'Error');
    return;
end
if length(z) < 10
    msgbox('rad_calc_ft_er: problem with array sizes', 'Error');
    return;
end
if length(omega) ~= 1
    msgbox('rad_calc_ft_er: problem with omega', 'Error');
    return;
end

%Put components into vectors
r = [x; y; z];
v = [vx; vy; vz];
n = [nx * ones(1, length(z)); ny * ones(1, length(z)); nz * ones(1, length(z))];
beta = v / C_LIGHT;

%Calculate the time from the z position
%vz = dz / dt
%dt = dz / vz
%t = int (1 / vz) dz
%t = int_calc_int_y_dx(z, (1 ./ vz));
t = cumtrapz(z, (1 ./ vz));

%Calculate the triple cross product in the integral
trip_cross_prod = cross(n, cross(n, beta));

%Debug
%tic;

%Calculate the exponential in the integral
expon = exp(-1i * omega * (t - dot(n, r) / C_LIGHT));

%Debug
% figure;
% hold on;
% plot(z, trip_cross_prod(1, :) / max(trip_cross_prod(1, :)), 'b-', 'linewidth', 2.);
% plot(z, real(expon), 'r-', 'linewidth', 1.);
% plot(z, imag(expon), 'g-', 'linewidth', 1.);
% title('Intermediate Steps In ft-er Calculation');
% ylabel('Calculated Quantities');
% xlabel('Z (m)');
% grid on;
% zoom on;
% hold off;

%Calculate the integral
%dt = (1 / vz) dz
% integ(1, :) = int_calc_int_y_dx(z, trip_cross_prod(1, :) .* expon .* (1 ./ vz));
% integ(2, :) = int_calc_int_y_dx(z, trip_cross_prod(2, :) .* expon .* (1 ./ vz));
% integ(3, :) = int_calc_int_y_dx(z, trip_cross_prod(3, :) .* expon .* (1 ./ vz));
integ(1, :) = cumtrapz(z, trip_cross_prod(1, :) .* expon .* (1 ./ vz));
integ(2, :) = cumtrapz(z, trip_cross_prod(2, :) .* expon .* (1 ./ vz));
integ(3, :) = cumtrapz(z, trip_cross_prod(3, :) .* expon .* (1 ./ vz));

%Calculate the fourier transform
ft_er = (1i * omega * Q_ELEC / (4 * pi * sqrt(2 * pi) * EPS_0 * C_LIGHT)) * integ;

%Debug
%dt = toc;
%disp(['dt = ', num2str(dt), ' sec']);

%Done
