%**************************************************************************
%This function calculates the H field in the magnet blocks.  It uses Ross
%Schlueter's method of finding the charge on the pole and setting it equal
%to the flux out of the pole.
%
%Input:
%b_r, remnant field in magnet blocks (T)
%mu_r, magnet block relative permeability
%lx_m, width of magnet blocks in x (m)
%ly_m, height of magnet blocks in y (m)
%lz_m, length of magnet blocks in z (m)
%lx_p, width of poles in x (m)
%ly_p, height of poles in y (m)
%lz_p, length of poles in z (m)
%g_u, undulator gap (m)
%
%Output:
%h_m, H in the magnet blocks (A/m)
%
%Zachary Wolf, 11/11/09

function [h_m] = field_calc_calc_hm(b_r, mu_r, lx_m, ly_m, lz_m, lx_p, ly_p, lz_p, g_u)

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

%Magnetic charge into half pole
%B = mu0 H + Br
%div B = 0 = mu0 div H + div Br
%div H = - div Br / mu0 == rho
%surf int H dA = vol int rho dV = Q
q = (b_r / (mu_r * mu0)) * lx_p * ly_p;

%Flux out of pole = charge on pole-----

%Flux through magnet at magnet center (V = 0 surface)
%phi_m = h_m * area of pole
%(phi / h_m == p_h)
p_h_m = lx_p * ly_p;

%Flux through gap at y=0 (V = 0 surface)
%hy = Sum hn sin(nkz) cosh(nky)
%hz = Sum hn cos(nkz) sinh(nky)
%phi_y = lx_p * int h_y dz
%phi_y = lx_p * int h0 sin(k z) dz
%phi_y = -lx_p (h0 / k) cos(k z)|_0 ^(l_u/4)
%phi_y = lx_p h0 / k
%hm = h0 sinh(k g/2)  %Hz continuous at magnet surface
%phi_y = (lx_p / k) hm / sinh(k g/2)
p_h_y = (lx_p / k_u) * (1 / sinh(k_u * g_u / 2));

%Flux through gap at z=0 (V = 0 surface)
%by = Sum bn sin(nkz) cosh(nky)
%bz = Sum bn cos(nkz) sinh(nky)
%phi_z = lx_p * int h_z dy
%phi_z = lx_p * int h0 sinh(k y) dy
%phi_z = lx_p h0 (1/k) cosh(k y)|_0 ^ g/2
%phi_z = lx_p h0 (1/k) [cosh(k g/2) - 1]
%hm = ho sinh(k g/2)
%phi_z = lx_p (1/k) [cosh(k g/2) - 1] hm (1 / sinh(k g/2))
p_h_z = (lx_p / k_u) * (cosh(k_u * g_u / 2) - 1) / sinh(k_u * g_u / 2);

%Total flux out of pole divided by h_m
p_h = p_h_m + p_h_y + p_h_z;

%Total flux = charge
%p_h * h_m = q
h_m = q / p_h;

%Done
