%**************************************************************************
%This function expands the measured field as N terms in the Fourier series.
%
%Input:
%z, position the field was measured at (m)
%b, measured field strength (T)
%z_start, start of region to analyze (m)
%z_end, end of region to analyze (m)
%lambda_u, undulator period (m)
%num_har, number of harmonics to calculate
%
%Output:
%a0, constant term (T)
%an, cosine terms (T)
%bn, sine terms (T)
%
%Zack Wolf  8/18/10

function [a0, an, bn] = field_calc_har(z, b, z_start, z_end, lambda_u, num_har)

%Method
%The region to analyze must be an integer number of periods.
%Expand the field in this region as
%b(z) = a0 + Sum{ an cos(n 2 pi z / lambda_u) + bn sin(n 2 pi z / lambda_u) }
%Calculate the coefficients
%a0 = (1 / L) Int{ b(z) dz }
%an = (2 / L) Int{ b(z) cos(n 2 pi z / lambda_u) }
%bn = (2 / L) Int{ b(z) sin(n 2 pi z / lambda_u) }

%Initialize
a0 = 0;
an = zeros(1, num_har);
bn = zeros(1, num_har);
  
%Check
if length(z) ~= length(b)
    msgbox('field_calc_har: z and b have different lengths', 'Error');
    return;
end
if (length(z) < 10) || (length(b) < 10)
    msgbox('field_calc_har: problem with array sizes', 'Error');
    return;
end
if (z_start < z(1)) || (z_end > z(end))
    %msgbox('field_calc_har: problem with parameters', 'Error');
    return;
end
if z_end <= z_start
    %msgbox('field_calc_har: problem with parameters', 'Error');
    return;
end

%Determine the range of integration
z1 = z_start;
num_per = floor((z_end - z_start) / lambda_u);
z2 = z1 + num_per * lambda_u;

%Isolate the data to integrate
index_array = find((z >= z1) & (z <= z2));
z_int = z(index_array);
b_int = b(index_array);
l_int = z2 - z1;

%Calculate a0
%integ = int_calc_int(z_int, b_int);
%integ = integ(end);
integ = trapz(z_int, b_int);
a0 = (1 / l_int) * integ;

%Calculate an and bn
for ii = 1 : num_har
    %integ = int_calc_int(z_int, b_int .* cos(ii * 2 * pi * z_int / lambda_u));
    %integ = integ(end);
    integ = trapz(z_int, b_int .* cos(ii * 2 * pi * z_int / lambda_u));
    an(ii) = (2 / l_int) * integ;
    %integ = int_calc_int(z_int, b_int .* sin(ii * 2 * pi * z_int / lambda_u));
    %integ = integ(end);
    integ = trapz(z_int, b_int .* sin(ii * 2 * pi * z_int / lambda_u));
    bn(ii) = (2 / l_int) * integ;
end

%Done
end
