%**************************************************************************
%This function generates the field from an EPU in the coordinate
%system where z is in the beam direction, y is up, and x is in the
%direction which makes a right handed system.  The EPU quadrants are at
%45 deg, 135 deg, 225 deg, and 315 deg.
%z=0 is at the first pole.
%
%Input:
%num_per, number of periods in the undulator, including ends
%lambda_u, undulator period (m)
%kr, parameter giving the radial exponential potential decrease (1/m)
%ks, parameter giving the azimuther potential roll off (1/m)
%phi_0, peak magnetic scalar potential (Tm)
%z_01, z-shift of quadrant 1 (m)
%z_02, z-shift of quadrant 2 (m)
%z_03, z-shift of quadrant 3 (m)
%z_04, z-shift of quadrant 4 (m)
%x, x-position to give the field at (m)
%y, y-position to give the field at (m)
%z, z-position to give the field at (m)
%
%Output:
%bx, Bx field at (x, y, z) (T)
%by, By field at (x, y, z) (T)
%bz, Bz field at (x, y, z) (T)
%
%Zachary Wolf
%4/20/15

function [bx, by, bz] = gen_epu_field(num_per, lambda_u, kr, ks, phi_0, z_01, z_02, z_03, z_04, x, y, z)
 
%Parameters
%d_taper = 0;  %(m)
%len_und = num_per * lambda_u;
%d_shift = 300e-6;  %(m)  Shift Q1 and Q4 up by d_shift, shift Q2 and Q3 down by d_shift
%d_shift = 0;  %(m)  Shift Q1 and Q4 up by d_shift, shift Q2 and Q3 down by d_shift

%Calculate the field at the given position from the quadrant fields-----

%Quadrant 1---

%Relative position to quadrant
r1 = (-x - y) / sqrt(2);
s1 = (x - y) / sqrt(2);
z1 = z - z_01;

% %Taper
% if (z >= 0) && (z <= len_und)
%     r1 = r1 + (d_taper / len_und) * z;
% end

% %Shift
% r1 = r1 + (d_shift / sqrt(2));
% s1 = s1 + (d_shift / sqrt(2));

%Calculate the fields from quadrant 1
[br1, bs1, bz1] = gen_epu_quadrant_field(num_per, lambda_u, kr, ks, phi_0, r1, s1, z1);
bx1 = (-br1 + bs1) / sqrt(2);
by1 = (-br1 - bs1) / sqrt(2);
bz1 = bz1;

% %Simulate sorting
% [mod_fact] = modulate_field(lambda_u, z1);
% bx1 = bx1 * mod_fact;
% by1 = by1 * mod_fact;
% bz1 = bz1 * mod_fact;

%Quadrant 2---

%Relative position to quadrant
r2 = (x - y) / sqrt(2);
s2 = (x + y) / sqrt(2);
z2 = z - z_02;

% %Taper
% if (z >= 0) && (z <= len_und)
%     r2 = r2 + (d_taper / len_und) * z;
% end

% %Shift
% r2 = r2 - (d_shift / sqrt(2));
% s2 = s2 + (d_shift / sqrt(2));

%Calculate the fields from quadrant 2
[br2, bs2, bz2] = gen_epu_quadrant_field(num_per, lambda_u, kr, ks, phi_0, r2, s2, z2);
bx2 = (br2 + bs2) / sqrt(2);
by2 = (-br2 + bs2) / sqrt(2);
bz2 = bz2;

% %Simulate sorting
% [mod_fact] = modulate_field(lambda_u, z2);
% bx2 = bx2 * mod_fact;
% by2 = by2 * mod_fact;
% bz2 = bz2 * mod_fact;

%Quadrant 3---

%Relative position to quadrant
r3 = (x + y) / sqrt(2);
s3 = (-x + y) / sqrt(2);
z3 = z - z_03;

% %Taper
% if (z >= 0) && (z <= len_und)
%     r3 = r3 + (d_taper / len_und) * z;
% end

% %Shift
% r3 = r3 + (d_shift / sqrt(2));
% s3 = s3 + (d_shift / sqrt(2));

%Calculate the fields from quadrant 3
[br3, bs3, bz3] = gen_epu_quadrant_field(num_per, lambda_u, kr, ks, -phi_0, r3, s3, z3);  %Opposite polarity blocks
bx3 = (br3 - bs3) / sqrt(2);
by3 = (br3 + bs3) / sqrt(2);
bz3 = bz3;

% %Simulate sorting
% [mod_fact] = modulate_field(lambda_u, z3);
% bx3 = bx3 * mod_fact;
% by3 = by3 * mod_fact;
% bz3 = bz3 * mod_fact;

%Quadrant 4---

%Relative position to quadrant
r4 = (-x + y) / sqrt(2);
s4 = (-x - y) / sqrt(2);
z4 = z - z_04;

% %Taper
% if (z >= 0) && (z <= len_und)
%     r4 = r4 + (d_taper / len_und) * z;
% end

% %Shift
% r4 = r4 - (d_shift / sqrt(2));
% s4 = s4 + (d_shift / sqrt(2));

%Calculate the fields from quadrant 4
[br4, bs4, bz4] = gen_epu_quadrant_field(num_per, lambda_u, kr, ks, -phi_0, r4, s4, z4);  %Opposite polarity blocks
bx4 = (-br4 - bs4) / sqrt(2);
by4 = (br4 - bs4) / sqrt(2);
bz4 = bz4;

% %Simulate sorting
% [mod_fact] = modulate_field(lambda_u, z4);
% bx4 = bx4 * mod_fact;
% by4 = by4 * mod_fact;
% bz4 = bz4 * mod_fact;

%-----

%Calculate the undulator fields by superposing the quadrant fields
bx = bx1 + bx2 + bx3 + bx4;
by = by1 + by2 + by3 + by4;
bz = bz1 + bz2 + bz3 + bz4;

%Done
end


%**************************************************************************
%Modulate the field to similate the sorting done to minimize phase errors.
%
%Zachary Wolf
%4/2/15

function [mod_fact] = modulate_field(lambda_u, z)

sin_mod = cos((2 * pi / (2 * lambda_u)) * z);

% if sin_mod >= 0
%     mod_fact = 1.01;
% else
%     mod_fact = .99;
% end

%mod_fact = 1 + .01 * sin_mod;
mod_fact = 1.;

end

%Done
