#
# 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]
# 
# Ring Membrane
#

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

#
# Zeros de BesselJ
#
read `/u1/rportugal/zeros/zerosj`: read `/u1/rportugal/zeros/ZerosAiry`:

ringmemb[BesselJY] := proc(m,x,radius_a,radius_b) 
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
BesselJ(m,x*radius_a)*BesselY(m,x*radius_b) - 
BesselY(m,x*radius_a)*BesselJ(m,x*radius_b)
end:


ringmemb[Rmn] := proc(m,n)
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
global radius_a,radius_b;
- BesselY(m,BesselJYZeros(m,n)*radius_a)*BesselJ(m,BesselJYZeros(m,n)*r) + 
BesselJ(m,BesselJYZeros(m,n)*radius_a)*BesselY(m,BesselJYZeros(m,n)*r)
end:


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

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

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

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

# Animagco dos Modos Normais
ringmemb[omega] := proc(m,n) 
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
BesselJYZeros(m,n)*vel_c
end:
protect('omega'):



ringmemb[animate_mode] := proc(m0::name=nonnegint,n0::name=posint)
local FRAM,Options,m1,n1;
global radius_a,radius_b,vel_c,r,theta,t;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if not type([radius_a,radius_b,vel_c],list('numeric')) 
    then ERROR(`The constants radius_a, radius_b 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; 
#RETURN(plots[animate](evalf(ringmemb[normal_mode11](m1,n1)),
#r=0..radius_b, 
#t=0..evalf(2*Pi*(FRAM-1)/(ringmemb[omega](m1,n1)*FRAM)),
#Options) );
plots[animate3d]([r,theta,evalf(ringmemb[normal_mode11](m1,n1))],
r=radius_a..radius_b, theta=0..2*Pi, 
t=0..evalf(2*Pi*(FRAM-1)/(ringmemb[omega](m1,n1)*FRAM)),
coords=cylindrical,
Options) 
end:

# 

`ringmemb/zeros` := proc(p::posint)
local z,A,B,R0,i;
global radius_a,radius_b,r,t,m,n;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if p=1 then RETURN(table([])) fi;
A := radius_a;
B := evalf(Pi/BesselJYZeros(0,p));
R0 := eval(ringmemb[Rmn](0,p));
for i to p-1 do
  z[i]:= fsolve(R0, r=A..A+B);
  A := z[i]
od;
eval(z)
end:

#   
#      As linhas nodais sco mostradas pelo seguinte procedimento:
ringmemb[nodal_lines] := proc(m0::name=nonnegint,n0::name=posint)
local ci,sl,p,radius,m1,n1;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
if not type([radius_a,radius_b],list('numeric')) 
  then ERROR(`radius_a and radius_b must be numeric.`) 
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[0] := plots[polarplot](radius_a,thickness=2);
ci[n1] := plots[polarplot](radius_b,thickness=2);
radius:=`ringmemb/zeros`(n1);
for p to n1-1 do ci[p] := plots[polarplot](radius[p]) od;
sl[0]:=NULL;
for p to m1 do 
  sl[p] := plots[polarplot]([r,(2*p-1)*Pi/(2*m1),r=radius_a..radius_b]);
  sl[-p] := plots[polarplot]([r,(2*p-1)*Pi/(2*m1),r=-radius_b..-radius_a]);
od;
plots[display]({seq(ci[p],p=0..n1),seq(sl[p],p=-m1..m1)},axes=none,
scaling=constrained,title=`Nodal lines of mode m = `.m1.` and n = `.n1)
end:
# 


# Animagco com Condigues Iniciais
# 
ringmemb[A1] := proc(m,n) 
global u0,v0,radius_a,radius_b,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`A1/A1uv`(m,n,eval(u0),radius_a,radius_b,vel_c,eval(AB_simp)) 
end:

`A1/A1uv` := proc(m,n,u0,radius_a,radius_b,vel_c,AB_simp)
local R,s;
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;;
if m=0 then s:=2 else s:=1 fi;
R:=ringmemb[Rmn](m,n);
AB_simp(-1/s/Pi*int(r*R*int(u0(r,theta)*cos(m*theta),theta=0..2*Pi),r=
radius_a..radius_b)/int(r*R^2,r =radius_a..radius_b));
end:

ringmemb[A2] := proc(m,n) 
global u0,v0,radius_a,radius_b,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`A2/A2uv`(m,n,eval(v0),radius_a,radius_b,vel_c,eval(AB_simp)) 
end:

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

ringmemb[B1] := proc(m,n) 
global u0,v0,radius_a,radius_b,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`B1/B1uv`(m,n,eval(u0),radius_a,radius_b,vel_c,eval(AB_simp)) 
end:

`B1/B1uv` := proc(m,n,u0,radius_a,radius_b,vel_c,AB_simp)
local R;
options remember,`Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;;
R:=ringmemb[Rmn](m,n);
AB_simp(-1/Pi*int(r*R*int(u0(r,theta)*sin(m*theta),theta=0..2*Pi),
r=radius_a..radius_b)/int(r*R^2,r =radius_a..radius_b));
end:

ringmemb[B2] := proc(m,n) 
global u0,v0,radius_a,radius_b,vel_c;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
`B2/B2uv`(m,n,eval(v0),radius_a,radius_b,vel_c,eval(AB_simp)) 
end:

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

# 
AB_simp := x->x: #x -> normal(simplify(x),expanded):


ringmemb[u] := proc(p1, p2) 
global m,n,r,theta,t;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
sum(-Rmn(0,n)*(
ringmemb[A1](0,n)*cos(BesselJYZeros(0,n)*vel_c*t)+
ringmemb[A2](0,n)*sin(BesselJYZeros(0,n)*vel_c*t))+
sum(-Rmn(m,n)*(
ringmemb[A1](m,n)*cos(m*theta)*cos(BesselJYZeros(m,n)*vel_c*t)+
ringmemb[A2](m,n)*cos(m*theta)*sin(BesselJYZeros(m,n)*vel_c*t)+
ringmemb[B1](m,n)*sin(m*theta)*cos(BesselJYZeros(m,n)*vel_c*t)+
ringmemb[B2](m,n)*sin(m*theta)*sin(BesselJYZeros(m,n)*vel_c*t)),
m = 1 .. p1),n = 1 .. p2) 
end:


#
# Procedure to animate with initial conditions
#
ringmemb[animate_ic] := proc(m0::name=nonnegint,n0::name=posint,t0::name=0..numeric) 
local FRAM,Options,p1,p2,tf; 
global radius_a, radius_b, 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,radius_b,vel_c],list(numeric))
     then ERROR(`The constants radius_a, radius_b 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(ringmemb[u](p1,p2))],r=radius_a..radius_b,
theta=0..2*Pi,t=0..tf,coords=cylindrical,Options)

end:
# 



############### ZEROS OF BESSELJY #######################

alias(k=BesselJYZeros):


# Zeros of BesselJY
_Envfdig := NULL:
_Envfdig1 := NULL:
_EnvBesselJYZeros := 'false':
BesselJYZeros:=proc(m::algebraic,n::algebraic)
global radius_a,radius_b; 
if nargs>3 or nargs<2 
  then ERROR(`wrong number of arguments`) fi;
if n=infinity or n=-infinity or not type(n,realcons) 
  then RETURN('procname(m,n)') fi;
if not type(n,nonnegint) 
  then ERROR(`2nd argument must be a positive integer`) fi;
if m=infinity or m=-infinity or not type(m,realcons) 
  then RETURN('procname(m,n)') fi;
if evalf(m)<0 then ERROR(`negative order not implemented`) fi;
if not assigned(BesselJZeros) 
  then ERROR(`The function BesselJZeros must be loaded`) fi;
if not type([radius_a,radius_b],list('numeric')) 
  then ERROR(`radius_a and radius_b must be numeric.`) fi;
_Envfdig1 := NULL;
if nargs=3
  then if args[3]='fulldigits' 
       then _Envfdig1 := 'fulldigits'
       else ERROR(`third argument must be fulldigits`)
       fi
fi;
_EnvBesselJYZeros := 'false'; 
'BesselJYZeros'(m,n)
end:

`evalf/BesselJYZeros` := proc(mm,n)
local m,M,N,i,x,q;
global radius_a,radius_b;
m:=evalf(mm);
if n=infinity or n=-infinity or not type(n,realcons) 
  then RETURN('BesselJYZeros'(m,n)) fi;
if not type(n,nonnegint) 
  then ERROR(`2nd argument must be a positive integer`) fi;
if m=infinity or m=-infinity or not type(m,realcons) 
  then RETURN('BesselJYZeros'(m,n)) fi;
if not type([radius_a,radius_b],list('numeric')) 
  then ERROR(`radius_a and radius_b must be numeric.`) fi;
if _EnvBesselJYZeros 
  then RETURN(fsolve(ringmemb[BesselJY](m,x,radius_a,radius_b),x,
                 _Envq1.._Envq2,_Envfdig1))
fi;
N:=1;
M:=1;
for i while N<=n do 
  _Envq1 := evalf(BesselJZeros(m,i))/radius_b;
  _Envq2 := evalf(BesselJZeros(m,i+1))/radius_b;
  q := evalf(BesselJZeros(m,M))/radius_a;
  if _Envq1-q<0 and q-_Envq2<0 
  then M:=M+1
  else _EnvBesselJYZeros := 'true';
       evalf('BesselJYZeros'(m,i));
       N:=N+1
  fi;
od;
evalf('BesselJYZeros'(m,n)) 
end:
protect(BesselJYZeros):

