input_case nil$

n:=2$

% Finite depth case
% -----------------
operator b; weight a=1$ wtlevel n+1$
for all x let cosh(x)^2=1+sinh(x)^2;

% general form of the velocity potential
F:=a*cos(xi)*cosh(h+z)+
for i:=2:n+1 sum a^i*
if evenp(i)
then for j:=2 step 2 until i sum b(i,j)*sin(j*xi)*cosh(j*(h+z))
else for j:=3 step 2 until i sum b(i,j)*cos(j*xi)*cosh(j*(h+z))$

fz:=df(F,z)$

% Taylor expansion
for all x such that not freeof(x,z) let cosh(x)=
begin scalar x0,x1,s,u,i; x0:=sub(z=0,x); x1:=x-x0;
    s:=cosh(x0); u:=i:=1; weight z=1;
    while (u:=u*x1/i) neq 0 do
    <<s:=s+u*(if evenp(i) then cosh(x0) else sinh(x0));i:=i+1>>;
    clear z; return s
end;

for all x such that not freeof(x,z) let sinh(x)=
begin scalar x0,x1,s,u,i; x0:=sub(z=0,x); x1:=x-x0;
    s:=sinh(x0); u:=i:=1; weight z=1;
    while (u:=u*x1/i) neq 0 do
    <<s:=s+u*(if evenp(i) then sinh(x0) else cosh(x0));i:=i+1>>;
    clear z; return s
end;

f:=F$ fz:=fz$
for all x such that not freeof(x,z) clear cosh(x),sinh(x);

% cosh and sinh of multiple arguments
for all n,x such that fixp(n) and n>1 let
    cosh(n*x)=cosh((n-1)*x)*cosh(x)+sinh((n-1)*x)*sinh(x),
    sinh(n*x)=cosh((n-1)*x)*sinh(x)+sinh((n-1)*x)*cosh(x);

f:=f$ fz:=fz$
for all n,x such that fixp(n) and n>1 clear cosh(n*x),sinh(n*x);

weight z=1$
depend xi,x,t; let df(xi,x)=1,df(xi,t)=-v; fx:=df(f,x)$

for all x let cos(x)^2=(1+cos(2*x))/2,
              sin(x)^2=(1-cos(2*x))/2;
for all x,y let cos(x)*cos(y)=(cos(x-y)+cos(x+y))/2,
                sin(x)*sin(y)=(cos(x-y)-cos(x+y))/2,
                sin(x)*cos(y)=(sin(x-y)+sin(x+y))/2;

% water level from the boundary condition
Z:=-df(f,t)-fz^2/2-fx^2/2$

repeat <<clear z;Z1:=sub(z=z1,Z);weight z=1;Z:=sub(z1=Z,Z1)>>
until freeof(Z,z);

clear z;

% this expression should vanish due to the boundary condition
q:=df(Z,t)-sub(z=Z,fz)+sub(z=Z,fx)*df(Z,x)$

clear f,fz,fx; clear df(xi,x),df(xi,t);
for all x clear cos(x)^2,sin(x)^2;
for all x,y clear cos(x)*cos(y),sin(x)*sin(y),sin(x)*cos(y);

weight u=2$
let v^2=sinh(h)/cosh(h)*(1+u); q:=q$ Z:=Z$ clear v^2;

let v=sqrt(sinh(h)/cosh(h))*
begin scalar a,b,j,i0; a:=1; b:=3/2; j:=0;
    i0:=if evenp(n) then n/2 else (n-1)/2;
    return 1+for i:=1:i0 sum <<b:=b-1;j:=j+1;a:=a*b/j*u>>
end;

q:=q$ Z:=Z$ clear v;

% resonant term
let z*cos(xi)=1; r:=q*z$ clear z*cos(xi);
r:=sub(z=0,r)$ q:=q-r*cos(xi)$ U:=u-r/(a*sinh(h))$

repeat <<clear u;U1:=sub(u=u1,U);weight u=2;U:=sub(u1=U,U1)>>
until freeof(U,u);

clear u; u:=U$ clear U; q:=q$ Z:=Z$
for all j let cos(j*xi)=z^j,sin(j*xi)=z^j; q:=q$
for all j clear cos(j*xi),sin(j*xi); clear a;
qa:=coeff(q,a)$ clear q; qa:=rest(rest(qa))$

for i:=2:n+1 do
<<  qi:=first(qa); qa:=rest(qa);
    j0:=if evenp(i) then 2 else 3;
    qi:=sub(z=sqrt(z),qi/z^j0); qz:=coeff(qi,z); clear qi;
    for j:=j0 step 2 until i do
    <<  qj:=first(qz); qz:=rest(qz);
        b(i,j):=-sub(b(i,j)=0,qj)/df(qj,b(i,j));
        clear qj;
    >>
>>;

factor a; on revpri,rat; % printing results
u:=u$ F:=F$ Z:=Z$
for all x clear cosh(x)^2;

% simplification of expressions with cosh and sinh
procedure simh1(x);
if x=0 then 0 else
begin scalar y,z; y:=x; z:=1;
    while remainder(y,sinh(h))=0 do <<y:=y/sinh(h);z:=z*sinh(h)>>;
    let sinh(h)^2=cosh(h)^2-1; y:=y; clear sinh(h)^2; y:=y*z;
    return y
end$

procedure simh(x);
begin scalar n,d,l,y,z; n:=num(x); d:=den(x);
    d:=simh1(d); l:=coeff(n,a); z:=1; n:=0;
    while l neq {} do
    <<y:=first(l);l:=rest(l);n:=n+simh1(y)*z;z:=z*a>>;
    return n/d
end$

1+simh(u); % omega^2
simh(F);   % velocity potential
simh(Z);   % water level

showtime;

end;