#
# Study of Membranes Vibration using Maple - version 1, Oct 1997.
#
# D. Frenkel[frenkel@venus.rdc.puc-rio.br] , 
# L. Golebiowski[lgolebi@netrio.com.br] , 
# R. Portugal[portugal@cat.cbpf.br]
# 
# Circular Membrane
#

with(plots):

#
# protect reserved variables
#
protect('t,r,theta,m,n');

#
# Zeros de BesselJ
#
read `/u1/rportugal/zeros/zerosj`: read `/u1/rportugal/zeros/ZerosAiry`:
alias(j=BesselJZeros):

# 
#
circmemb[omega] := (m,n) -> k(m,n)*vel_c:
protect('omega');
#  
circmemb[k] := (m,n) -> BesselJZeros(m,n)/radius_a:
protect('k');
# 


circmemb[normal_mode11] := proc (m, n)  
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
  BesselJ(m,k(m,n)*r)*cos(k(m,n)*vel_c*t)*cos(m*theta) 
end:

circmemb[normal_mode12] := proc (m, n)  
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
  BesselJ(m,k(m,n)*r)*sin(k(m,n)*vel_c*t)*cos(m*theta) 
end:

circmemb[normal_mode21] := proc (m, n)  
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
  BesselJ(m,k(m,n)*r)*cos(k(m,n)*vel_c*t)*sin(m*theta) 
end:

circmemb[normal_mode22] := proc (m, n)  
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
  BesselJ(m,k(m,n)*r)*sin(k(m,n)*vel_c*t)*sin(m*theta) 
end:
# 
# 

# 
# Normal mode animation
# 
circmemb[animate_mode] := proc(m0::name=nonnegint,n0::name=posint)
local FRAM,Options,m1,n1;
global radius_a,vel_c,r,theta,t;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if not type([radius_a,vel_c],list('numeric')) 
    then ERROR(`The constants radius_a and vel_c must be numerical.`) 
fi;
m1:=subs(m0,m);
n1:=subs(n0,n);
if not type(m1,'nonnegint') 
    then ERROR(`wrong argument:`,m0) 
fi;
if not type(n1,'posint') 
    then ERROR(`wrong argument:`,n0) 
fi;

if nargs=1
    then Options:='frames'=8,'style=patch'; FRAM:=8;
    else Options:=args[3..nargs];
         if has({Options},'frames')
           then FRAM:=subs({Options},'frames');
           else Options:=Options,'frames'=8; FRAM:=8
         fi;
         if not has({Options},'style')
           then Options:=Options,'style=patch'
         fi;
         if not has({Options},'constrained')
           then Options:=Options,'scaling=constrained'
         fi;
fi; 
plots[animate3d]([r,theta,evalf(circmemb[normal_mode11](m1,n1))],r=0..radius_a,theta=0..2*Pi,t=0..
evalf(2*Pi*(FRAM-1)/(circmemb[omega](m1,n1)*FRAM)),coords=cylindrical,Options) 
end:

#
# NODAL LINES
#
circmemb[nodal_lines] := proc(m0::name=nonnegint,n0::name=posint)
local ci,sl,p,m1,n1;
global radius_a,r;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if not assigned(BesselJZeros) 
    then ERROR(`The procedure BesselJZeros must be loaded.`)
fi;
if not type(radius_a,'numeric') 
    then ERROR(`The constant radius_a must be numerical.`) 
fi;
m1:=subs(m0,m);
n1:=subs(n0,n);
if not type(m1,'nonnegint') 
    then ERROR(`wrong argument:`,m0) 
fi;
if not type(n1,'posint') 
    then ERROR(`wrong argument:`,n0) 
fi;
ci[n1]:=plots[polarplot](radius_a,thickness=2);
for p to n1-1 do ci[p] := plots[polarplot](radius_a*BesselJZeros(0,p)/BesselJZeros(0,n1)) od;
for p to m1 do sl[p] := plots[polarplot]([r,(2*p-1)*Pi/(2*m1),r=-radius_a..radius_a])
od;
plots[display]({seq(ci[p],p=1..n1),seq(sl[p],p=1..m1)},axes=none,
scaling=constrained,title=`Nodal lines of mode m = `.m1.` and n = `.n1)
end:

# 
# Animation with initial conditions
# 
circmemb[A1] :=proc(m,n) global u0,v0,radius_a,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`A1/A1uv`(m,n,u0,radius_a,vel_c,eval(AB_simp)) 
end:
`A1/A1uv` := proc(m,n,u0,radius_a,vel_c,AB_simp)
local s;
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if m=0 then s:=2 else s:=1 fi;
AB_simp(int(r*BesselJ(m,BesselJZeros(m,n)/radius_a*r)*int(u0*
cos(m*theta),theta=0..2*Pi),r=0..radius_a)/
(s*Pi*radius_a^2*BesselJ(m+1,BesselJZeros(m,n))^2));
end:

circmemb[A2] :=proc(m,n) global u0,v0,radius_a,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`A2/A2uv`(m,n,v0,radius_a,vel_c,eval(AB_simp)) 
end:
`A2/A2uv` := proc(m,n,v0,radius_a,vel_c,AB_simp)
local s;
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if m=0 then s:=2 else s:=1 fi;
AB_simp(int(r*BesselJ(m,k(m,n)*r)*int(v0*cos(m*theta),theta=0..2*Pi
),r=0..radius_a)/
(s*circmemb[omega](m,n)*Pi*radius_a^2*(BesselJ(m+1,k(m,n)*radius_a))^2));
end:

circmemb[B1] :=proc(m,n) global u0,v0,radius_a,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`B1/B1uv`(m,n,u0,radius_a,vel_c,eval(AB_simp)) 
end:
`B1/B1uv` := proc(m,n,u0,radius_a,vel_c,AB_simp)
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
AB_simp(int(r*BesselJ(m,BesselJZeros(m,n)/radius_a*r)*int(u0*
sin(m*theta),theta=0..2*Pi),r=0..radius_a)/(Pi*radius_a^2*BesselJ(m+1,BesselJZeros(m,n))^2));
end:

circmemb[B2] :=proc(m,n) global u0,v0,radius_a,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`B2/B2uv`(m,n,v0,radius_a,vel_c,eval(AB_simp)) 
end:
`B2/B2uv` := proc(m,n,v0,radius_a,vel_c,AB_simp)
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
AB_simp(int(r*BesselJ(m,k(m,n)*r)*int(v0*sin(m*theta),theta=0..2*Pi
),r=0..radius_a)/
(circmemb[omega](m,n)*Pi*radius_a^2*(BesselJ(m+1,k(m,n)*radius_a))^2));
end:
protect('A1,A2,B1,B2');

AB_simp := simplify:



circmemb[u] := proc(p1, p2)  
global m,n;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
2*sum(BesselJ(0,BesselJZeros(0,n)/radius_a*r)*
(cos(BesselJZeros(0,n)/radius_a*vel_c*t)*A1(0,n)+A2(0,n)*
sin(BesselJZeros(0,n)/radius_a*vel_c*t))+
sum(BesselJ(m,BesselJZeros(m,n)/radius_a*r)*(A1(m,n)*
cos(BesselJZeros(m,n)/radius_a*vel_c*t)*cos(m*theta)+
A2(m,n)*sin(BesselJZeros(m,n)/radius_a*vel_c*t)*cos(m*theta)+
B1(m,n)*cos(BesselJZeros(m,n)/radius_a*vel_c*t)*sin(m*theta)+
sin(m*theta)*sin(BesselJZeros(m,n)/radius_a*vel_c*t)*B2(m,n)),
m = 1 .. p1),n = 1 .. p2) 
end:
protect('u');

#
# INITIAL CONDITIONS
# 
circmemb[animate_ic] := proc(m0::name=nonnegint,n0::name=posint,t0::name=0..numeric) 
local FRAM,Options,p1,p2,tf; 
global radius_a,vel_c,r,theta,t; 
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if nargs<3 
then ERROR(`This procedure has at least 3 arguments. The first one is \
trunc_m=truncation_value, the second is trunc_n=truncation_value and \
the third is the time range t=0..final_value.`) 
fi;
if not type([radius_a,vel_c],list(numeric))
     then ERROR(`The constants radius_a and vel_c must be numeric.`) 
elif not assigned(BesselJZeros)
     then ERROR(`The procedure BesselJZeros must be loaded.`) 
fi;

p1:=subs(m0,'trunc_m');
p2:=subs(n0,'trunc_n');
if not type(p1,nonnegint) 
    then ERROR(`wrong argument:`,m0) 
fi;
if not type(p2,posint) 
    then ERROR(`wrong argument:`,n0) 
fi;

tf := op(2,op(2,t0));
if not(type(tf,numeric) and tf>0)
    then ERROR(`wrong argument:`,t0) 
fi;

Options := args[4 .. nargs];
 
if not has({Options},'style') 
  then Options:=Options,'style=patch'
fi;
if not has({Options},'constrained')
  then Options:=Options,'scaling=constrained'
fi;

plots[animate3d]([r,theta,evalf(circmemb[u](p1,p2))],r=0..radius_a,theta=0..2*Pi,
t=0..tf,coords=cylindrical,Options)

end:
# 
# 
