function [a, b, delay] = gethosewave(phase, ga, w0, length, terms);

xt= [40.70 8.731 0.9235 0.1553 129.2/2.5, 0.94;
    26.45 6.798 0.6998 0.1564 92.33/2.5, 0.915;
	15.40 8.546 1.032 0.1551 92.85/2.5, 0.90;
	9.427 6.322 0.7744 0.1557 66.05/2.5, 0.86;
	7.413 4.633 0.5349 0.1578 48.48/2.5, 0.82;
	5.262 4.180 0.5384 0.1576 42.95/2.5, 0.79;
	4.833 3.251 0.3726 0.1612 33.52/2.5, 0.75;
	4.219 2.385 0.2890 0.1755 26.42/2.5, 0.70];

pt = [1 1.5 2 3 4 5 6 8];
x = interp1(pt, xt, phase, 'linear', 'extrap');

if terms < 3
    terms = 3;
end
a(2) = x(1);
a(3) = x(2);
a(4) = x(3);
if terms >=5
    for i=5:terms
        a(i)=a(i-1)*x(4);
    end
end
a = a(1:terms);
a(1)=-sum(a(2:end));
b(2)=x(5);
b(1)=2.5*b(2);
if terms >= 3
    for i=3:terms
        b(i)=b(i-1)/2.5;
    end
end
b = b(1:terms);
a=a/sum(a./b);

c=1/ga;
ph=phase*pi/180;
p=1/w0*(w0*length/c)^(1/(1-(2*ph/pi)));
a=a/p;
b=b/p;
delay=x(6)*p;

