%**************************************************************************
%This function calculates the angular spectrum of the energy from 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:
%ang_spect, energy per unit frequency per unit solid angle (J s)
%
%Zachary Wolf
%7/9/10

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

%Method
%See Jackson, p. 670, eq. 14.65

%Physical constants
global Q_ELEC;          %Electron charge (C)
global M_ELEC;          %Electron mass (kg)
global C_LIGHT;         %Speed of light (m/s)
global EPS_0;           %Permittivity of free space (F/m)
global H_BAR;           %Planck constant / 2pi (Js)
  
%Check
if (length(z) ~= length(vx)) || (length(z) ~= length(vy)) || (length(z) ~= length(vz))
    msgbox('rad_calc_ang_spect: arrays have different lengths', 'Error');
    return;
end
if length(z) < 10
    msgbox('rad_calc_ang_spect: problem with array sizes', 'Error');
    return;
end

%Calculate the acceleration
[ax, ay, az] = rad_calc_acc(z, vx, vy, vz);

%Put components into vectors
r = [x; y; z];
v = [vx; vy; vz];
a = [ax; ay; az];
n = [nx * ones(1, length(z)); ny * ones(1, length(z)); nz * ones(1, length(z))];
beta = v / C_LIGHT;
beta_dot = a / 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));

%Calculate the numerator in the integral
num = cross(n, cross((n - beta), beta_dot));

%Calculate the denominator in the integral
denom = (1 - dot(n, beta)).^2;

%Loop over omega's
ang_spect = zeros(1, length(omega));
for ii = 1 : length(omega)

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

%Calculate the integral
%dt = (1 / vz) dz
vals = int_calc_int_y_dx(z, (num(1, :) ./ denom) .* expon .* (1 ./ vz));
integ(1) = vals(end);
vals = int_calc_int_y_dx(z, (num(2, :) ./ denom) .* expon .* (1 ./ vz));
integ(2) = vals(end);
vals = int_calc_int_y_dx(z, (num(3, :) ./ denom) .* expon .* (1 ./ vz));
integ(3) = vals(end);

%Calculate the magnitude squared of the integral
mag_integ_sq = integ * integ';  %' is hermetian conjugate

%Calculate the angular spectrum
ang_spect(ii) = (Q_ELEC^2 / (16 * pi^3 * EPS_0 * C_LIGHT)) * mag_integ_sq;

%End loop over omega's
end

%Done
