%**************************************************************************
%This function calculates the field expansion coefficients in a planar
%undulator using the technique of Walker, NIM A237 (1985), 366.
%
%Input:
%h_m, H in the magnet blocks (A/m)
%lz_m, length of magnet blocks in z (m)
%lz_p, length of poles in z (m)
%g_u, undulator gap (m)
%num_bn, number of bn coeffieients to calculate
%plt_file, plot file, use to illustrate calculations
%
%Output:
%bn(1 to num_bn), field expansion coefficients
%
%Zachary Wolf, 11/11/09

function [bn] = field_calc_calc_bn(h_m, lz_m, lz_p, g_u, num_bn, plt_file)

%Constants
mu0 = 4 * pi * 10^-7;  %permeability of free space (Tm/A)

%Parameters
l_u = 2 * lz_m + 2 * lz_p;
k_u = 2 * pi / l_u;

%Set the range of z-values
z = - l_u / 2 : .000001 : l_u / 2;

%Field expansion
%by = Sum bn sin(nkz) cosh(nky)
%bz = Sum bn cos(nkz) sinh(nky)
%At y = g_u/2, assume rectangular shape to bz, bz=0 under poles, bz = mu0*h_m next to magnets
%(1/l_u) Int bz cos(nkz) = bn sinh(nkg/2) 1/2
%bn = (2 / sinh(nkg/2)) (1/l_u) Int bz cos(nkz)
%Format in z: mmppppmmmmppppmm
bz_expand = zeros(1, length(z));
index_array = find(z < - lz_m / 2 - lz_p);
bz_expand(index_array) = - mu0 * h_m;
index_array = find((z > - lz_m / 2) & (z < lz_m / 2));
bz_expand(index_array) = mu0 * h_m;
index_array = find(z > lz_m / 2 + lz_p);
bz_expand(index_array) = - mu0 * h_m;
for nn = 1 : num_bn
    integral = trapz(z, bz_expand .* cos(nn * k_u * z));
    bn(nn) = (2 ./ sinh(nn * k_u * g_u / 2)) * (1 / l_u) * integral;
    %Debug
    if nn <= 5
        %field_calc_plot_field_expansion(plt_file, nn, z, bz_expand, cos(nn * k_u * z));
    end
end

%Done
