# drumhead problem, General case, Problem 4.3-3
# u_{tt} = c^2(laplacian u), u=0 on clamped edge
# circle radius=a, c=wave speed, 
# position = f(r,theta) = (a^2-r^2)*r*sin(theta),
# velocity = g(r,theta) = 1 
a:=2;c:=1;N:=4;
f:=(r,theta)->(a^2-r^2)*r*sin(theta);
g:=(r,theta)->1;
# Because g=1, then 4.3 equations (18) and (19) have value zero [integrate on theta].
# This implies a-star_mn=0 except for m=0. And b-star_mn=0.
# Modify EXAMPLE 2, Section 4.3, to define u1(r,theta,t)
# Then a_mn=0 for all subscripts. And b_mn=0 except for m=1.
# Compute b_1n from (14) section 4.3 [we use function b1(n)]
# Computing u1(r,theta,t)
alpha1:=evalf([BesselJZeros(1,1..N)]);
b1:=n->2*2*(a^5)*4/(alpha1[n]^3*BesselJ(2,alpha1[n]));
u1:=(r,theta,t)->sin(theta)*sum(b1(n)*BesselJ(1,alpha1[n]*r)*cos(alpha1[n]*c*t),n=1..N);
# Computing u2(r,theta,t)
alpha0:=evalf([BesselJZeros(0,1..N)]);
astar0:=n->(2*Pi/(Pi*c*alpha0[n]*a*BesselJ(1,alpha0[n])^2))*int(BesselJ(0,alpha0[n]*r/a)*r, r=0..a);
u2:=(r,theta,t)->cos(0*theta)*sum(astar0(n)*BesselJ(0,alpha0[n]*r/a)*sin(alpha0[n]*c*t/a),n=1..N);
# Define u=u1+u2
u:=(r,theta,t)->u1(r,theta,t)+u2(r,theta,t);
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plot3d(f(r,theta),r=0..a,theta=0..2*Pi,coords=z_cylindrical,orientation=[13,81]);
plot3d(g(r,theta),r=0..a,theta=0..2*Pi,coords=z_cylindrical);
plots[animate](plot3d,
[u(r,theta,t),r=0..a,theta=0..2*Pi,coords=z_cylindrical,orientation=[13,81]],t=0..1.5);