%**************************************************************************
%This function calculates the fields when the EPU is in circular polarization
%left hand mode.
%
%Input:
%lambda_u, undulator period (m)
%num_per, number of periods in the center of the undulator, even
%phi0, overall scale factor for the fields (Tm)
%kx, parameter giving the x-dependence (1/m)
%ky, parameter giving the y-dependence (1/m)
%kz, parameter giving the z-dependence, kz = ku (1/m)
%xc, x-position of the magnetic center (m)
%yc, y-position of the magnetic center (m)
%zc, z-position of the magnetic center (m)
%x, x-position to calculate the field at (m)
%y, y-position to calculate the field at (m)
%z, z-position to calculate the field at (m)
%
%Output:
%bx, x-component of the field strength (T)
%by, y-component of the field strength (T)
%bz, z-component of the field strength (T)
%
%Zachary Wolf
%4/2/13

function [bx, by, bz] = cp_lh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z)

%Check the constraint on the k parameters
if (kx^2 + ky^2 < (0.98 * kz)^2) || (kx^2 + ky^2 > (1.02 * kz)^2)
    error('The k parameters do not satisfy the constraint');
end

%Check that num_per is even
if mod(num_per, 2) > 0
    error('The number of periods must be even');
end

%Calculate the fields
bx_cos = (1 / (2 * sqrt(2))) * phi0 * kz * cosh((1 / sqrt(2)) * kz * ((x - xc) + (y - yc))) .* cos(kz * (z - zc));
bx_sin = (1 / (2 * sqrt(2))) * phi0 * kz * cosh((1 / sqrt(2)) * kz * ((x - xc) - (y - yc))) .* sin(kz * (z - zc));
by_cos = (1 / (2 * sqrt(2))) * phi0 * kz * cosh((1 / sqrt(2)) * kz * ((x - xc) + (y - yc))) .* cos(kz * (z - zc));
by_sin = - (1 / (2 * sqrt(2))) * phi0 * kz * cosh((1 / sqrt(2)) * kz * ((x - xc) - (y - yc))) .* sin(kz * (z - zc));
bz_sin = -(1 / 2) * phi0 * kz * sinh((1 / sqrt(2)) * kz * ((x - xc) + (y - yc))) .* sin(kz * (z - zc));
bz_cos = (1 / 2) * phi0 * kz * sinh((1 / sqrt(2)) * kz * ((x - xc) - (y - yc))) .* cos(kz * (z - zc));

%Sine terms-----
%Make the fields zero outside the undulator
z_end_us = zc + (num_per / 2) * lambda_u + lambda_u;
z_end_ds = zc - (num_per / 2) * lambda_u - lambda_u;
index_array = find((z < z_end_ds) | (z > z_end_us));
bx_sin(index_array) = 0;
by_sin(index_array) = 0;
bz_sin(index_array) = 0;

%Make the first half period quarter strength
z_hp1_us = z_end_us - lambda_u / 2;
z_hp1_ds = z_end_ds + lambda_u / 2;
index_array = find(((z > z_end_ds) & (z < z_hp1_ds)) | ((z < z_end_us) & (z > z_hp1_us)));
bx_sin(index_array) = bx_sin(index_array) * 1/4;
by_sin(index_array) = by_sin(index_array) * 1/4;
bz_sin(index_array) = bz_sin(index_array) * 1/4;

%Make the second half period three-quarter strength
z_hp2_us = z_hp1_us - lambda_u / 2;
z_hp2_ds = z_hp1_ds + lambda_u / 2;
index_array = find(((z > z_hp1_ds) & (z < z_hp2_ds)) | ((z < z_hp1_us) & (z > z_hp2_us)));
bx_sin(index_array) = bx_sin(index_array) * 3/4;
by_sin(index_array) = by_sin(index_array) * 3/4;
bz_sin(index_array) = bz_sin(index_array) * 3/4;

%Cosine terms-----
%Make the fields zero outside the undulator
z_end_us = zc + (num_per / 2) * lambda_u + lambda_u + lambda_u / 4;
z_end_ds = zc - (num_per / 2) * lambda_u - lambda_u + lambda_u / 4;
index_array = find((z < z_end_ds) | (z > z_end_us));
bx_cos(index_array) = 0;
by_cos(index_array) = 0;
bz_cos(index_array) = 0;

%Make the first half period quarter strength
z_hp1_us = z_end_us - lambda_u / 2;
z_hp1_ds = z_end_ds + lambda_u / 2;
index_array = find(((z > z_end_ds) & (z < z_hp1_ds)) | ((z < z_end_us) & (z > z_hp1_us)));
bx_cos(index_array) = bx_cos(index_array) * 1/4;
by_cos(index_array) = by_cos(index_array) * 1/4;
bz_cos(index_array) = bz_cos(index_array) * 1/4;

%Make the second half period three-quarter strength
z_hp2_us = z_hp1_us - lambda_u / 2;
z_hp2_ds = z_hp1_ds + lambda_u / 2;
index_array = find(((z > z_hp1_ds) & (z < z_hp2_ds)) | ((z < z_hp1_us) & (z > z_hp2_us)));
bx_cos(index_array) = bx_cos(index_array) * 3/4;
by_cos(index_array) = by_cos(index_array) * 3/4;
bz_cos(index_array) = bz_cos(index_array) * 3/4;

%Combine fields
bx = bx_sin + bx_cos;
by = by_sin + by_cos;
bz = bz_sin + bz_cos;

%Done
