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


#
# protect reserved variables
# 
protect('x,y,t,m,n');
 
rectmemb[omega] := proc(m,n) 
(m^2/side_a^2+n^2/side_b^2)^(1/2)*Pi*vel_c 
end:

rectmemb[normal_mode1] := proc(m,n) 
sin(Pi*m/side_a*x)*sin(Pi*n/side_b*y)*cos(rectmemb[omega](m,n)*t)
end:

rectmemb[normal_mode2] := proc(m,n) 
sin(Pi*m/side_a*x)*sin(Pi*n/side_b*y)*sin(rectmemb[omega](m,n)*t)
end:

#
# Animation of the Normal Modes
#  
rectmemb[animate_mode] := proc(m0::name=posint,n0::name=posint)
local FRAM,Options,m1,n1;
global side_a,side_b,x,y,t;
if not type([side_a,side_b,vel_c],list(numeric))
     then ERROR(`The constants side_a, side_b and vel_c must be numerical.`) 
elif not assigned(rectmemb[omega])
     then ERROR(`The function omega must be defined.`) 
elif not assigned(normal_mode1)
    then ERROR(`The function normal_mode1 must be defined.`)
fi;
m1:=subs(m0,m);
n1:=subs(n0,n);
if not type(m1,'posint') 
    then ERROR(`wrong argument:`,m0) 
fi;
if not type(n1,'posint') 
    then ERROR(`wrong argument:`,n0) 
fi;

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},'scaling')
 then Options:=Options,'scaling=constrained'
fi;

plots[animate3d](evalf(rectmemb[normal_mode1](m1,n1)), x=0..side_a, y=0..side_b,
t=0..evalf(2*Pi*FRAM/rectmemb[omega](m1,n1)/(FRAM+1)), Options)
end:


# 
# Degeneracy Analysis
# 
# The following procedure does the animation of the linear combination 
# c[1]*normal_mode(m,n)+c[2]*normal_mode(n,m)
# The procedure must be called in the form animate_lc( mode(m,n), c1, c2):
#
rectmemb[animate_lc] := proc(m0::name=posint,n0::name=posint,c1::numeric,c2::numeric)
local FRAM,Options,m1,n1;
global side_a,side_b,x,y,t;
if not type([side_a,side_b,vel_c],list('numeric')) 
  then ERROR(`The constants side_a, b and vel_c must be numerical.`) 
elif not assigned(normal_mode1) 
  then ERROR(`The function normal_mode1 must be defined.`) 
elif side_a<>side_b 
  then ERROR(`The membrane must be square.`)
elif nargs<4 
  then ERROR(`The first and second arguments are of the form m=m0 and n=0. \
The second and the third are the coefficients of the linear combination (numbers).`)
fi;

m1:=subs(m0,'m');
n1:=subs(n0,'n');
if not type(m1,'posint') 
  then ERROR(`wrong argument:`,m0) 
fi;
if not type(n1,'posint') 
  then ERROR(`wrong argument:`,n0) 
fi;

Options:=args[5..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},'scaling')
 then Options:=Options,'scaling=constrained'
fi;

plots[animate3d](c1*rectmemb[normal_mode1](m1,n1)+c2*rectmemb[normal_mode1](n1,m1),
   x=0..side_a,
   y=0..side_b,
   t=0..evalf(2*Pi*FRAM/rectmemb[omega](m1,n1)/(FRAM+1)), Options)
end:


#
# The procedure nodal_line plots the nodal line. 
# The arguments of  this procedure are the same as the procedure animate_lc. 
#
rectmemb[nodal_lines] := proc(m0::name=posint,n0::name=posint,c1::numeric,c2::numeric)
local m1,n1,Options,s;
global side_a,side_b,x,y,t;
if not type([side_a,side_b],list('numeric')) 
  then ERROR(`The constants side_a, side_b and vel_c must be numerical.`) 
elif not assigned(normal_mode1) 
  then ERROR(`The function normal_mode1 must be defined.`) 
elif side_a<>side_b 
  then ERROR(`The membrane must be square.`)
elif nargs<4 
  then ERROR(`The first and second arguments are of the form m=m0 and n=0. \
The second and the third are the coefficients of the linear combination (numbers).`)
fi;
m1:=subs(m0,'m');
n1:=subs(n0,'n');
if not type(m1,'posint') 
  then ERROR(`wrong argument:`,m0) 
fi;
if not type(n1,'posint') 
  then ERROR(`wrong argument:`,n0) 
fi;
Options := args[5..nargs];
if not has({Options},'color')
  then Options:=Options,'color=black'
fi;
s:=subs(t=0,c1*rectmemb[normal_mode1](m1,n1)+c2*rectmemb[normal_mode1](n1,m1));
plots[implicitplot](s,x=0..side_a,y=0..side_b,'scaling=constrained','xtickmarks'
=[0,side_a],'ytickmarks'=[0,side_b],Options);
end:

# 
# 
# Animation with Initial Conditions
# 
rectmemb[u] := proc(trunc_m::{name,posint},trunc_n::{name,posint})
global side_a,side_b,x,y,t;
options `Copyright 1997 by D.Frenkel, L.Golebiowski and R.Portugal`;
sum(sum(sin(Pi*m/side_a*x)*sin(Pi*n/side_b*y)*
(cos((m^2*side_b^2+n^2*side_a^2)^(1/2)*Pi*vel_c/side_b/side_a*t)*A(m,n)+
sin((m^2*side_b^2+n^2*side_a^2)^(1/2)*Pi*vel_c/side_b/side_a*t)*B(m,n)),
m = 1 .. trunc_m),n = 1 .. trunc_n)
end:
protect('u'):


# 
rectmemb[A] := proc(m,n) 
global u0,v0,side_a,side_b;
`A/Auv`(m,n,u0,v0,side_a,side_b,eval(AB_simp)) 
end:
protect(A):

`A/Auv` := proc(m,n,u0,v0,side_a,side_b,AB_simp)
global x,y;
options remember;
AB_simp(4/side_a/side_b*int(int(u0*sin(m*Pi*x/side_a)*
sin(n*Pi*y/side_b),x=0..side_a),y=0..side_b));
end:

rectmemb[B] := proc(m,n) 
global u0,v0,side_a,side_b;
`B/Buv`(m,n,u0,v0,side_a,side_b,eval(AB_simp)) 
end:
protect(B):

`B/Buv` := proc(m,n,u0,v0,side_a,side_b,AB_simp)
global x,y;
options remember;
AB_simp(4/side_a/side_b/rectmemb[omega](m,n)*int(int(v0*sin(m*Pi*x/side_a)*
sin(n*Pi*y/side_b),x=0..side_a),y=0..side_b)) 
end:

AB_simp := simplify:

# 
rectmemb[animate_ic] := proc(m0::name=posint,n0::name=posint,t0::name=0..numeric) 
local FRAM, Options,p1,p2,tf; 
global side_a, side_b,x,y,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([side_a,side_b,vel_c],list(numeric))
     then ERROR(`The constants side_a, side_b and vel_c must be numerical.`) 
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](evalf(rectmemb[u](p1,p2)), x=0..side_a, y=0..side_b, t=0..tf, Options) 
end:
# 
# 
