function t = stab220_BW(theta,mu)
t = stab210_BW(theta,mu) + 3*mu*mu*(1 - mu).*exp(-i.*theta).*(exp(i.*theta)-1).*(1 - exp(-i.*theta)).^2;

