%**************************************************************************
%gen_byvsz
%
%This function generates the main field in an undulator.
%The convention is that the field is in the Y direction, the beam is moving in
%the Z direction, and the trajectory is in the X direction.
%The data is generated for the beam axis with X = 0 and Y = 0.
%The point Z = 0 is at the undulator midpoint.
%The output is written to a file.
%
%Zachary Wolf  7/31/03, 1/3/06, 2/19/07

%Clear the workspace
clear;

%Set the working directory
cd 'c:\Analysis\Insertion Device Analysis\Generate Fake Data\';

%Parameters
B0_no_err = 1.249471;           %Peak field strength on the midplane (T)
LAMBDA_U = .03;                 %Undulator period (m)
NUM_POLE = 226;                 %Number of poles per jaw, even (for these ends)
SAMPLE_SPACING = .000200;       %Distance between magnetic field samples (m)
GAMMA = 26692.70;               %Electron Lorentz factor
DAT_FILE = 'c:\magdata\fake_data\fake_byvsz.dat';    %Output data file
PAR_FILE = 'c:\magdata\fake_data\fake_byvsz.par';    %Output parameter file

%Errors
B_ERR_ALL_POLE = 0;             %Switch, fractional error on all poles
B_ERR_SGL_POLE = 0;             %Switch, fractional error on single pole
B_ERR_GAUSS = 0;                %Switch, Gaussian error on all poles
B_ERR_PHASE = 0;                %Switch, error on two poles making phase error
B_ERR_BASE = 0;                 %Switch, error over a range of poles simulating base shim error
B_ERR_EARTH = 0;                %Switch, add constant field to all field values

%Error parameters
B_ERR_ALL_POLE_FRAC = -.00054;  %Fractional field change to all field values
B_ERR_SGL_POLE_NUM = 51;        %Number of pole for single pole error
B_ERR_SGL_POLE_FRAC = 0.01;     %Value of fractional error on single pole
B_ERR_GAUSS_DIST_RMS = 0.0001;  %RMS variation of gaussian error field added to peak pole field, zero mean (T)
B_ERR_PHASE_POLE_NUM = 51;      %Number of first pole for a phase error
B_ERR_PHASE_FRAC = 0.1;         %Fractional change in the strength of two neighboring poles
B_ERR_BASE_POLE_NUM_START = 51; %Number of first pole for a base error
B_ERR_BASE_POLE_NUM_END = 70;   %Number of lase pole for a base error
B_ERR_BASE_FRAC = 0.01;         %Fractional change in the strength of the poles
B_ERR_EARTH_VAL = -0.00003;     %Earth's field (T)

%Physical constants
Q_ELEC = -1.60217733e-19;       %Electron charge (C)
M_ELEC = 9.1093897e-31;         %Electron mass (kg)
C_LIGHT = 2.99792458e8;         %Speed of light (m/s)

%Initialize
B0 = B0_no_err;

%Generate z positions in the undulator starting from the end
z = 0 : SAMPLE_SPACING : NUM_POLE * LAMBDA_U / 2.;

%Calculate the pole positions
z_pole = [1 : NUM_POLE] * LAMBDA_U / 2. - LAMBDA_U / 4.;

%Generate a sinusoidal field for the undulator
b = B0 * sin(2 * pi * z / LAMBDA_U);

%Make the first pole quarter strength
index = find(abs(z - z_pole(1)) < LAMBDA_U / 4.);
b(index) = b(index) * .25;

%Make the second pole three quarter strength
index = find(abs(z - z_pole(2)) < LAMBDA_U / 4.);
b(index) = b(index) * .75;

%Make the NUM_POLE-1'th pole three quarter strength
index = find(abs(z - z_pole(NUM_POLE - 1)) < LAMBDA_U / 4.);
b(index) = b(index) * .75;

%Make the NUM_POLE'th pole quarter strength
index = find(abs(z - z_pole(NUM_POLE)) < LAMBDA_U / 4.);
b(index) = b(index) * .25;

%Error, fractional change on all poles
if B_ERR_ALL_POLE ~= 0
    b = b * (1 + B_ERR_ALL_POLE_FRAC);
    B0 = B0 * (1 + B_ERR_ALL_POLE_FRAC);
end

%Error, fractional change on a single pole
if B_ERR_SGL_POLE ~= 0
    index = find(abs(z - z_pole(B_ERR_SGL_POLE_NUM)) < LAMBDA_U / 4.);
    b(index) = b(index) * (1 + B_ERR_SGL_POLE_FRAC);
end

%Error, Gaussian distribution added to all poles, random noise
if B_ERR_GAUSS ~= 0
    b = b + B_ERR_GAUSS_DIST_RMS * randn(1, length(b));
end

%Error, phase error
if B_ERR_PHASE ~= 0
    index = find(abs(z - z_pole(B_ERR_PHASE_POLE_NUM)) < LAMBDA_U / 4.);
    b(index) = b(index) * (1 + B_ERR_PHASE_FRAC);
    index = find(abs(z - z_pole(B_ERR_PHASE_POLE_NUM + 1)) < LAMBDA_U / 4.);
    b(index) = b(index) * (1 + B_ERR_PHASE_FRAC);
end

%Error, base error
if B_ERR_BASE ~= 0
    for ii = B_ERR_BASE_POLE_NUM_START : B_ERR_BASE_POLE_NUM_END
        index = find(abs(z - z_pole(ii)) < LAMBDA_U / 4.);
        b(index) = b(index) * (1 + B_ERR_BASE_FRAC);
    end
end

%Generate 15 periods of zeros for the ends
z_end = 0 : SAMPLE_SPACING : (30 + 1) * LAMBDA_U / 2.;
b_end = zeros(1, length(z_end));

%Add the ends to the undulator field
z_gen = z_end;
z_gen = [z_gen, z_gen(end)+ SAMPLE_SPACING + z];
z_gen = [z_gen, z_gen(end)+ SAMPLE_SPACING + z_end];
b_gen = [b_end, b, b_end];

%Make z = 0 the center of the undulator
z_meas = z_gen - (z_gen(1) + z_gen(end)) / 2.;
b_meas = b_gen;

%Error, add the Earth's field to all samples
if B_ERR_EARTH ~= 0
    b_meas = b_meas + B_ERR_EARTH_VAL;
end

%Write the data to a file
fid = fopen(DAT_FILE, 'w');
fprintf(fid, 'SLAC Magnetic Measurements\r\n');
fprintf(fid, 'Date: mm-dd-yyyy\r\n');
fprintf(fid, 'Time: hh:mm:ss\r\n');
fprintf(fid, '\r\n');
fprintf(fid, 'Project: Test\r\n');
fprintf(fid, 'Magnet Type: Undulator\r\n');
fprintf(fid, 'Magnet Name: Test\r\n');
fprintf(fid, 'Serial Number:\r\n');
fprintf(fid, 'Run Number: 1\r\n');
fprintf(fid, 'Measurement Device: Hall probe\r\n');
fprintf(fid, 'Operator: zw\r\n');
fprintf(fid, 'Comment: Test\r\n');
fprintf(fid, '\r\n');
fprintf(fid, '     Z           By     \r\n');
fprintf(fid, '    (m)          (T)    \r\n');
fprintf(fid, '-----------  -----------\r\n');
fprintf(fid, '%11.6f  %11.6f\r\n', [z_meas; b_meas]);
fclose(fid);

%Write the parameters to a file
fid = fopen(PAR_FILE, 'w');
fprintf(fid, 'SLAC Magnetic Measurements\r\n');
fprintf(fid, 'Simulated data to test undulator analysis software.\r\n');
fprintf(fid, '\r\n');
fprintf(fid, 'Simulated Data Parameters:\r\n');
fprintf(fid, 'B0_no_err = %10.5f T\r\n', B0_no_err);
fprintf(fid, 'GAMMA = %10.2f\r\n', GAMMA);
fprintf(fid, 'LAMBDA_U = %10.5f m\r\n', LAMBDA_U);
fprintf(fid, 'NUM_POLE = %5i\r\n', NUM_POLE);
fprintf(fid, 'SAMPLE_SPACING = %10.7f m\r\n', SAMPLE_SPACING);
fprintf(fid, 'N_SAMPLES = %7i\r\n', length(z_meas));
fprintf(fid, 'Magnetic Field Length = %10.6f m\r\n', z(end) - z(1));
fprintf(fid, 'Measurement Length = %10.6f m\r\n', z_meas(end) - z_meas(1));
fprintf(fid, '\r\n');
fprintf(fid, 'Errors:\r\n');
if B_ERR_ALL_POLE ~= 0
    fprintf(fid, 'B_ERR_ALL_POLE_FRAC = %10.6f\r\n', B_ERR_ALL_POLE_FRAC);
end
if B_ERR_SGL_POLE ~= 0
    fprintf(fid, 'B_ERR_SGL_POLE_NUM = %4i\r\n', B_ERR_SGL_POLE_NUM);
    fprintf(fid, 'B_ERR_SGL_POLE_FRAC = %10.6f\r\n', B_ERR_SGL_POLE_FRAC);
end
if B_ERR_GAUSS ~= 0
    fprintf(fid, 'B_ERR_GAUSS_DIST_RMS = %10.6f T\r\n', B_ERR_GAUSS_DIST_RMS);
end
if B_ERR_PHASE ~= 0
    fprintf(fid, 'B_ERR_PHASE_POLE_NUMS = %4i, %4i\r\n', B_ERR_PHASE_POLE_NUM, B_ERR_PHASE_POLE_NUM + 1);
    fprintf(fid, 'B_ERR_PHASE_FRAC = %10.6f\r\n', B_ERR_PHASE_FRAC);
end
if B_ERR_BASE ~= 0
    fprintf(fid, 'B_ERR_BASE_POLE_NUMS = %4i to %4i\r\n', B_ERR_BASE_POLE_NUM_START, B_ERR_BASE_POLE_NUM_END);
    fprintf(fid, 'B_ERR_BASE_FRAC = %10.6f\r\n', B_ERR_BASE_FRAC);
end
if B_ERR_EARTH ~= 0
    fprintf(fid, 'B_ERR_EARTH_VAL = %11.7f T\r\n', B_ERR_EARTH_VAL);
end
fprintf(fid, 'End of error list.\r\n');
fprintf(fid, '\r\n');
fprintf(fid, 'Computed Values:\r\n');
K = LAMBDA_U * abs(Q_ELEC) * B0 / (2 * pi * M_ELEC * C_LIGHT);
fprintf(fid, 'K = %10.6f\r\n', K);
lambda_r = (LAMBDA_U / (2 * GAMMA^2)) * (1 + K^2 / 2);
fprintf(fid, 'lambda_r = %15.6e m\r\n', lambda_r);
fclose(fid);

%Done
