%**************************************************************************
%gen_bxyvsz
%This function generates the main field in a planar undulator
%and a small perpendicular error field.
%The convention is that the main field is in the Y direction,
%the perpendicular small error field is in the X direction,
%the beam is moving in the Z direction,
%The main 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_bxyvsz.dat';    %Output data file
PAR_FILE = 'c:\magdata\fake_data\fake_bxyvsz.par';    %Output parameter file

%Bx errors
BX_ERR_SIN = 0;                 %Switch, sinusoidal field variation through the undulator
BX_ERR_EARTH = 1;               %Switch, add constant field to all field values
BX_ERR_NOISE = 1;               %Switch, add random noise to each sample

%Bx error parameters
BX_ERR_SIN_AMP = 0.0001;        %Amplitude of the sinusoidal field variation
BX_ERR_SIN_HAR = 1;             %Harmonic number for the sinusoidal field variation
BX_ERR_EARTH_VAL = 0.00002;     %Earth's field (T)
BX_ERR_NOISE_VAL = 0.0001;      %Sigma of gaussian noise, mean zero (T)

%By errors
BY_ERR_ALL_POLE = 0;            %Switch, fractional error on all poles
BY_ERR_SGL_POLE = 0;            %Switch, fractional error on single pole
BY_ERR_GAUSS = 0;               %Switch, Gaussian error on all poles
BY_ERR_PHASE = 0;               %Switch, error on two poles making phase error
BY_ERR_BASE = 0;                %Switch, error over a range of poles simulating base shim error
BY_ERR_EARTH = 1;               %Switch, add constant field to all field values
BY_ERR_NOISE = 1;               %Switch, add random noise to each sample

%By error parameters
BY_ERR_ALL_POLE_FRAC = -.00054; %Fractional field change to all field values
BY_ERR_SGL_POLE_NUM = 51;       %Number of pole for single pole error
BY_ERR_SGL_POLE_FRAC = 0.01;    %Value of fractional error on single pole
BY_ERR_GAUSS_DIST_RMS = 0.0001; %RMS variation of gaussian error field added to peak pole field, zero mean (T)
BY_ERR_PHASE_POLE_NUM = 51;     %Number of first pole for a phase error
BY_ERR_PHASE_FRAC = 0.1;        %Fractional change in the strength of two neighboring poles
BY_ERR_BASE_POLE_NUM_START = 51; %Number of first pole for a base error
BY_ERR_BASE_POLE_NUM_END = 70;  %Number of lase pole for a base error
BY_ERR_BASE_FRAC = 0.01;        %Fractional change in the strength of the poles
BY_ERR_EARTH_VAL = -0.00002;    %Earth's field (T)
BY_ERR_NOISE_VAL = 0.0001;      %Sigma of gaussian noise, mean zero (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)

%Undulator Center -----

%Z ---
%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.;

%Bx ---
%Generate the ideal Bx
bx = zeros(1, length(z));

%By ---
%Initialize
B0 = B0_no_err;

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

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

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

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

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

%Undulator Center Errors -----

%Bx ---
%Generate a sinusoidal Bx field for the undulator
if BX_ERR_SIN ~= 0
    bx = bx + BX_ERR_SIN_AMP * sin(BX_ERR_SIN_HAR * pi * z / z(end));
end

%By ---
%Error, fractional change on all poles
if BY_ERR_ALL_POLE ~= 0
    by = by * (1 + BY_ERR_ALL_POLE_FRAC);
    B0 = B0 * (1 + BY_ERR_ALL_POLE_FRAC);
end

%Error, fractional change on a single pole
if BY_ERR_SGL_POLE ~= 0
    index = find(abs(z - z_pole(BY_ERR_SGL_POLE_NUM)) < LAMBDA_U / 4.);
    by(index) = by(index) * (1 + BY_ERR_SGL_POLE_FRAC);
end

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

%Error, phase error
if BY_ERR_PHASE ~= 0
    index = find(abs(z - z_pole(BY_ERR_PHASE_POLE_NUM)) < LAMBDA_U / 4.);
    by(index) = by(index) * (1 + BY_ERR_PHASE_FRAC);
    index = find(abs(z - z_pole(BY_ERR_PHASE_POLE_NUM + 1)) < LAMBDA_U / 4.);
    by(index) = by(index) * (1 + BY_ERR_PHASE_FRAC);
end

%Error, base error
if BY_ERR_BASE ~= 0
    for ii = BY_ERR_BASE_POLE_NUM_START : BY_ERR_BASE_POLE_NUM_END
        index = find(abs(z - z_pole(ii)) < LAMBDA_U / 4.);
        by(index) = by(index) * (1 + BY_ERR_BASE_FRAC);
    end
end

%Drift Regions -----

%Z ---
%Generate positions outside the undulator
z_end = 0 : SAMPLE_SPACING : (30 + 1) * LAMBDA_U / 2.;

%Add the ends to the undulator positions
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];

%Bx ---
%Generate zeros outside the undulator
bx_end = zeros(1, length(z_end));

%Add the ends to the Bx undulator field
bx_gen = [bx_end, bx, bx_end];

%By ---
%Generate zeros outside the undulator
by_end = zeros(1, length(z_end));

%Add the ends to the By undulator field
by_gen = [by_end, by, by_end];

%Measured Values -----

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

%Bx -----
%Switch from generated to measured coordinates
bx_meas = bx_gen;

%By -----
%Switch from generated to measured coordinates
by_meas = by_gen;

%Measured Value Errors -----

%Bx ---
%Error, add the Earth's field to all samples
if BX_ERR_EARTH ~= 0
    bx_meas = bx_meas + BX_ERR_EARTH_VAL;
end

%Error, random noise added to all the samples
if BX_ERR_NOISE ~= 0
    bx_meas = bx_meas + BX_ERR_NOISE_VAL * randn(1, length(bx_meas));
end

%By -----
%Error, add the Earth's field to all samples
if BY_ERR_EARTH ~= 0
    by_meas = by_meas + BY_ERR_EARTH_VAL;
end

%Error, random noise added to all the samples
if BY_ERR_NOISE ~= 0
    by_meas = by_meas + BY_ERR_NOISE_VAL * randn(1, length(by_meas));
end

%Output -----

%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           Vy           Vx           By           Bx     \r\n');
fprintf(fid, '    (m)          (V)          (V)          (T)          (T)    \r\n');
fprintf(fid, '-----------  -----------  -----------  -----------  -----------\r\n');
fprintf(fid, '%11.6f  %11.6f  %11.6f  %11.6f  %11.6f\r\n', [z_meas; by_meas; bx_meas; by_meas; bx_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, 'Bx Errors:\r\n');
if BX_ERR_SIN ~= 0
    fprintf(fid, 'BX_ERR_SIN_HAR = %4i\r\n', BX_ERR_SIN_HAR);
    fprintf(fid, 'BX_ERR_SIN_AMP = %10.6f\r\n', BX_ERR_SIN_AMP);
end
if BX_ERR_EARTH ~= 0
    fprintf(fid, 'BX_ERR_EARTH_VAL = %11.7f T\r\n', BX_ERR_EARTH_VAL);
end
if BX_ERR_NOISE ~= 0
    fprintf(fid, 'BX_ERR_NOISE_VAL = %11.7f T\r\n', BX_ERR_NOISE_VAL);
end
fprintf(fid, 'End of Bx error list.\r\n');
fprintf(fid, '\r\n');
fprintf(fid, 'By Errors:\r\n');
if BY_ERR_ALL_POLE ~= 0
    fprintf(fid, 'BY_ERR_ALL_POLE_FRAC = %10.6f\r\n', BY_ERR_ALL_POLE_FRAC);
end
if BY_ERR_SGL_POLE ~= 0
    fprintf(fid, 'BY_ERR_SGL_POLE_NUM = %4i\r\n', BY_ERR_SGL_POLE_NUM);
    fprintf(fid, 'BY_ERR_SGL_POLE_FRAC = %10.6f\r\n', BY_ERR_SGL_POLE_FRAC);
end
if BY_ERR_GAUSS ~= 0
    fprintf(fid, 'BY_ERR_GAUSS_DIST_RMS = %10.6f T\r\n', BY_ERR_GAUSS_DIST_RMS);
end
if BY_ERR_PHASE ~= 0
    fprintf(fid, 'BY_ERR_PHASE_POLE_NUMS = %4i, %4i\r\n', BY_ERR_PHASE_POLE_NUM, BY_ERR_PHASE_POLE_NUM + 1);
    fprintf(fid, 'BY_ERR_PHASE_FRAC = %10.6f\r\n', BY_ERR_PHASE_FRAC);
end
if BY_ERR_BASE ~= 0
    fprintf(fid, 'BY_ERR_BASE_POLE_NUMS = %4i to %4i\r\n', BY_ERR_BASE_POLE_NUM_START, BY_ERR_BASE_POLE_NUM_END);
    fprintf(fid, 'BY_ERR_BASE_FRAC = %10.6f\r\n', BY_ERR_BASE_FRAC);
end
if BY_ERR_EARTH ~= 0
    fprintf(fid, 'BY_ERR_EARTH_VAL = %11.7f T\r\n', BY_ERR_EARTH_VAL);
end
if BY_ERR_NOISE ~= 0
    fprintf(fid, 'BY_ERR_NOISE_VAL = %11.7f T\r\n', BY_ERR_NOISE_VAL);
end
fprintf(fid, 'End of By 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
