%**************************************************************************
%This function performs an analysis related to undulator radiation.
%
%Input:
%param, parameter structure
%dat_file, output data file name
%plt_file, output plot file name
%field, field structure
%traj, trajectory structure
%phas, phase structure
%
%Output:
%rad, structure containing radiation information
%   omega_fund          %Fundamental angular frequency (1/s)
%	omega_har           %Array of fundamental and harmonic frequencies, _har quantities calculated at these frequencies in forward direction (1/s)
%	energy_kev_har      %Energy of fundamental and harmonics (keV)
%	ang_num_spect_har   %Number of photons per second for 0.1 A beam current per 0.1% bandwidth per square mrad solid angle (num/s/mrad^2/0.1%bw/0.1A) 
%	ang_num_spect_ideal_har     %Center period repeated with no phase errors (num/s/mrad^2/0.1%bw/0.1A)
%	efficiency_har      %Ratio of measured angular number spectrum to ideal
%	ex_fund             %Fourier transform Ex(omega_fund) (V/m)
%	ey_fund             %Fourier transform Ey(omega_fund) (V/m)
%	ex_mag_fund         %Magnitude of ex_fund (V/m)
%	ey_mag_fund         %Magnitude of ey_fund (V/m)
%	exy_ang_fund        %Phase angle_Ey - angle_Ex (rad)
%   s0, s1, s2, s3      %Stokes parameters S0 to S3 (V/m)^2
%	p1, p2, p3          %Polarization rates P1 to P3
%
%Zachary Wolf
%7/10/10

function [rad] = epu_anal_rad(param, dat_file, plt_file, field, traj, phas)

%Physical constants
global C_LIGHT;         %Speed of light (m/s)
global H_BAR;           %Planck constant / 2pi (Js)
global Q_ELEC;          %Electron charge (C)

%Extract required information
x = traj.x;
y = traj.y;
z = traj.z;
vx = traj.vx;
vy = traj.vy;
vz = traj.vz;
lambda_r = phas.lambda_r;
lambda_u = field.lambda_u;
%num_period = floor(param.num_poles / 2);
%num_period = floor((param.num_poles - (1/2) * 2 * param.num_end_poles) / 2);  %Assume end poles are 1/2 efficient at producing radiation
num_period = floor(max(length(field.bx_peak), length(field.by_peak)) / 2);
gamma = param.gamma;
a0_bx = field.a0_bx;
an_bx = field.an_bx;
bn_bx = field.bn_bx;
a0_by = field.a0_by;
an_by = field.an_by;
bn_by = field.bn_by;
    
%Initialize
rad.omega_fund = 0;
rad.omega_har = 0;
rad.energy_kev_har = 0;
rad.ang_num_spect_har = 0;
rad.ang_num_spect_ideal_har = 0;
rad.efficiency_har = 0;
rad.ex_fund = 0;
rad.ey_fund = 0;
rad.ex_mag_fund = 0;
rad.ey_mag_fund = 0;
rad.exy_ang_fund = 0;
rad.s0 = 0;
rad.s1 = 0;
rad.s2 = 0;
rad.s3 = 0;
rad.p1 = 0;
rad.p2 = 0;
rad.p3 = 0;

%Check
if lambda_u <= 0
    return;
end
if lambda_r <= 0
    return;
end

% %Use the average forward velocity in the undulator for vz for the radiation calculations
% k_eff = phas.k_eff;
% vz0 = C_LIGHT * (1 - (1 / (2 * gamma^2)) * (1 + (1/2) * k_eff^2));
% vz0 = C_LIGHT * lambda_u / (lambda_r + lambda_u);
% vz = vz0 * ones(1, length(z));

%Define n, only use forward direction for standard results
nx = 0;
ny = 0;
nz = 1;

%Define the fundamental angular frequency
omega_fund = C_LIGHT * 2 * pi / lambda_r;

%Define the fundamental and harmonic frequencies
omega_har = (1 : 15) * omega_fund;

%Calculate the energy of each harmonic
energy_kev_har = omega_har * H_BAR / abs(Q_ELEC) / 1000.;

%Calculate the angular number spectrum at each harmonic
[ang_num_spect_har] = rad_calc_ang_num_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega_har);

%Calculate the ideal angular number spectrum at each harmonic
%Average period repeated with no phase errors
[ang_num_spect_ideal_har] = rad_calc_ang_num_spect_ideal(x, y, z, vx, vy, vz, nx, ny, nz, omega_har, omega_fund, num_period, lambda_u);
%Fourier coefficients of field
%[ang_num_spect_ideal_har] = rad_calc_ang_num_spect_har(a0_bx, an_bx, bn_bx, a0_by, an_by, bn_by, nx, ny, nz, gamma, omega_har, omega_fund, num_period, lambda_u);

%Calculate the efficiency of producing light at each harmonic
efficiency_har = ang_num_spect_har ./ ang_num_spect_ideal_har;

%Calculate the fourier transform of the electric field at the fundamental
[ft_er] = rad_calc_ft_er(x, y, z, vx, vy, vz, nx, ny, nz, omega_fund);
ft_er_und = ft_er(:, end);
ex_fund = ft_er_und(1);
ey_fund = ft_er_und(2);
ex_mag_fund = abs(ex_fund);
ey_mag_fund = abs(ey_fund);
exy_ang_fund = angle(ey_fund * conj(ex_fund));

%Calculate the Stokes parameters
[s0, s1, s2, s3, p1, p2, p3] = rad_calc_stokes(ex_fund, ey_fund);

%Plot Fourier Transform of E-----

%Calculate the fourier transform of the electric field at the fundamental
[ft_er] = rad_calc_ft_er(x, y, z, vx, vy, vz, nx, ny, nz, omega_fund);
[ft_er_ideal] = rad_calc_ft_er_ideal(ft_er, z, lambda_u);

%Plot the fourier transform of the electric field at the fundamental
%rad_plot_ft_er_x(plt_file, ft_er);
%rad_plot_ft_er_y(plt_file, ft_er);
rad_plot_ft_er_x_ideal(plt_file, ft_er, ft_er_ideal);
rad_plot_ft_er_y_ideal(plt_file, ft_er, ft_er_ideal);

%Scan Many Harmonics-----

%Define a frequency scan range over many harmonics
omega = (1:150) / 10 * omega_fund;
%omega = (1:10000) / 1000 * omega_fund;

% %Calculate the angular energy spectrum
% [ang_ener_spect] = rad_calc_ang_ener_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
% 
% %Plot the angular energy spectrum
% rad_plot_ang_ener_spect(plt_file, omega_fund, omega, ang_ener_spect);

%Calculate the angular number spectrum
[ang_num_spect] = rad_calc_ang_num_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);

%Plot the angular number spectrum
rad_plot_ang_num_spect(plt_file, omega_fund, omega, ang_num_spect);

%Scan First Harmonic-----

%Define a frequency scan range over the first harmonic
omega = (1 - 0.04) * omega_fund : 0.001 * omega_fund : (1 + 0.04) * omega_fund;

%Debug
% for ii = 1 : 4 : length(omega)
%     [ft_er] = rad_calc_ft_er(x, y, z, vx, vy, vz, nx, ny, nz, omega(ii));
%     [ft_er_ideal] = rad_calc_ft_er_ideal(ft_er, z, lambda_u);
%     rad_plot_ft_er_x_ideal(plt_file, ft_er, ft_er_ideal);
%     %rad_plot_ft_er_y_ideal(plt_file, ft_er, ft_er_ideal);
% end

% %Calculate the angular energy spectrum
% [ang_ener_spect] = rad_calc_ang_ener_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
% 
% %Plot the angular energy spectrum
% rad_plot_ang_ener_spect(plt_file, omega_fund, omega, ang_ener_spect);

%Calculate the angular number spectrum
[ang_num_spect] = rad_calc_ang_num_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
%(Average period repeated with no phase errors)
[ang_num_spect_ideal] = rad_calc_ang_num_spect_ideal(x, y, z, vx, vy, vz, nx, ny, nz, omega, omega_fund, num_period, lambda_u);
%(Fourier coefficients of field)
%[ang_num_spect_ideal] = rad_calc_ang_num_spect_har(a0_bx, an_bx, bn_bx, a0_by, an_by, bn_by, nx, ny, nz, gamma, omega, omega_fund, num_period, lambda_u);

%Plot the angular number spectrum
%rad_plot_ang_num_spect(plt_file, omega_fund, omega, ang_num_spect);
rad_plot_ang_num_spect_ideal(plt_file, omega_fund, omega, ang_num_spect, ang_num_spect_ideal);

%Scan Third Harmonic-----

%Define a frequency scan range over the third harmonic
omega = (3 - 0.04) * omega_fund : 0.001 * omega_fund : (3 + 0.04) * omega_fund;

% %Calculate the angular energy spectrum
% [ang_ener_spect] = rad_calc_ang_ener_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
% 
% %Plot the angular energy spectrum
% rad_plot_ang_ener_spect(plt_file, omega_fund, omega, ang_ener_spect);

%Calculate the angular number spectrum
[ang_num_spect] = rad_calc_ang_num_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
%(Average period repeated with no phase errors)
[ang_num_spect_ideal] = rad_calc_ang_num_spect_ideal(x, y, z, vx, vy, vz, nx, ny, nz, omega, omega_fund, num_period, lambda_u);
%(Fourier coefficients of field)
%[ang_num_spect_ideal] = rad_calc_ang_num_spect_har(a0_bx, an_bx, bn_bx, a0_by, an_by, bn_by, nx, ny, nz, gamma, omega, omega_fund, num_period, lambda_u);

%Plot the angular number spectrum
%rad_plot_ang_num_spect(plt_file, omega_fund, omega, ang_num_spect);
rad_plot_ang_num_spect_ideal(plt_file, omega_fund, omega, ang_num_spect, ang_num_spect_ideal);

%Scan Fifth Harmonic-----

%Define a frequency scan range over the fifth harmonic
omega = (5 - 0.04) * omega_fund : 0.001 * omega_fund : (5 + 0.04) * omega_fund;

% %Calculate the angular energy spectrum
% [ang_ener_spect] = rad_calc_ang_ener_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
% 
% %Plot the angular energy spectrum
% rad_plot_ang_ener_spect(plt_file, omega_fund, omega, ang_ener_spect);

%Calculate the angular number spectrum
[ang_num_spect] = rad_calc_ang_num_spect(x, y, z, vx, vy, vz, nx, ny, nz, omega);
%(Average period repeated with no phase errors)
[ang_num_spect_ideal] = rad_calc_ang_num_spect_ideal(x, y, z, vx, vy, vz, nx, ny, nz, omega, omega_fund, num_period, lambda_u);
%(Fourier coefficients of field)
%[ang_num_spect_ideal] = rad_calc_ang_num_spect_har(a0_bx, an_bx, bn_bx, a0_by, an_by, bn_by, nx, ny, nz, gamma, omega, omega_fund, num_period, lambda_u);

%Plot the angular number spectrum
%rad_plot_ang_num_spect(plt_file, omega_fund, omega, ang_num_spect);
rad_plot_ang_num_spect_ideal(plt_file, omega_fund, omega, ang_num_spect, ang_num_spect_ideal);

%Results-----

%Save results
rad.omega_fund = omega_fund;
rad.omega_har = omega_har;
rad.energy_kev_har = energy_kev_har;
rad.ang_num_spect_har = ang_num_spect_har;
rad.ang_num_spect_ideal_har = ang_num_spect_ideal_har;
rad.efficiency_har = efficiency_har;
rad.ex_fund = ex_fund;
rad.ey_fund = ey_fund;
rad.ex_mag_fund = ex_mag_fund;
rad.ey_mag_fund = ey_mag_fund;
rad.exy_ang_fund = exy_ang_fund;
rad.s0 = s0;
rad.s1 = s1;
rad.s2 = s2;
rad.s3 = s3;
rad.p1 = p1;
rad.p2 = p2;
rad.p3 = p3;

%Write results to the data file
rad_dat_rad_summary(dat_file, param, rad);

%Done
