(*:Name: Chladni` *)

(*:Author:        Wilberd van der Kallen, 2006 *)

(*:Mathematica Version: 5.2 *)
(*:Package Version: 1.1 *)

(*:Summary:                                                               *)
(* We introduce functions that model some Chladni plate experiments.      *)

(*:Sources:                                                               *)
(* Wheatstone, On the Figures Obtained by Strewing Sand on Vibrating      *)
(* Surfaces, Commonly Called Acoustic Figures, Philosophical Transactions *)
(* of the Royal Society of London, vol 123, 1833.                         *)
(* Lord Rayleigh, Theory of Sound, Dover, New York 1945.                  *)
(* Mary Waller, Chladni Figures, A Study in Symmetry, Bell, London 1961.  *)
(* Rossing en Fletcher, Principles of Vibration and Sound, Springer,      *)
(* New York 1995.                                                         *)

(* :Discussion:                                                           *)
(* Consider a snapshot of an oscillating thin bending bar with free ends, *)
(* in an eigen mode.                                                      *)
(* The equations are                                                      *)
(* u''''= freq^2 u                                                        *)
(* with boundary conditions                                               *)
(*  u'''(1) = u'''(-1) = u''(1) = u''(-1) = 0.                            *)
(* Here freq is proportional to the frequency of the oscillation.         *)
(* The proportionality constant depends on things like the                *)
(* stiffness of the bar.                                                  *)
(*                                                                        *)
(* Next consider a snapshot for an oscillating thin                       *)
(* round plate of radius one with free boundary, in an eigen mode.        *)
(* We simplify the boundary conditions.                                   *)
(* If Delta denotes the Laplacian, then the equations are                 *)
(* Delta Delta u = freq^2 u                                                *)
(* with our simplified boundary conditions saying that both               *)
(* Delta u                                                                *)
(* and its gradient vanish on the boundary.                               *)
(* One gets a model that, given its simplicity, agrees reasonably well    *)
(* with reality, both in shape of Chladni pattern and in frequency.       *)
(*                                                                        *)
(* Next consider a snapshot of an oscillating thin bending square plate   *)
(* with free boundary. It is too hard to solve the equations              *)
(* symbolically. By way of approximation one could work with the pretend  *)
(* solution f(x)g(y) + (or -) f(y)g(x), where f and g are solutions       *)
(* for the bending bar. This pretend solution does not solve any          *)
(* equation. If f and g satisfy the free boundary conditions for the      *)
(* bending bar, then the product f(x)g(y) is no eigenfunction of          *)
(* Delta Delta.  But if one removes the exponentials from f and g,        *)
(* keeping only the sines and cosines, then one does get an eigenfunction *)
(* of Delta Delta, with a frequency that is the sum of the frequencies    *)
(* associated to f and g. Thus one gets a reasonable choice of frequency  *)
(* to associate with the pretend solution.                                *)
(* One gets a very rough model that shows some similarity with reality,   *)
(* both in shape of Chladni pattern and in frequency. This goes back to   *)
(* Wheatstone.                                                            *)

BeginPackage["Chladni`"]

bar::usage = "bar[m] gives a snapshot of an oscillating thin \
bending bar with free ends at -1 and 1."

freqbar::usage = 
     "freqbar[m] is proportional with the frequency associated with \
 bar[m]."
 
freqsquare::usage = 
     "freqsquare[m,n] = freqbar[m] + freqbar[n]"
 
disc::usage = "disc[m,n][r] gives the radial part of \
a snapshot of an oscillating thin disc of \
radius one with simplified free boundary. The angular part is Cos[m phi]."

freqdisc::usage = "freqdisc[m,n] is proportional with the \
frequency associated with disc[m,n]."

(*:Examples:

In[1]:=  <<chladni.m

In[2]:=  Plot[bar[3][t], {t, -1, 1},
  PlotLabel -> "disc(" <> ToString[3] <> ") freq=" <> ToString[freqbar[3]]] 

Out[2]= -Graphics-

In[3]:= ContourPlot[
    bar[7][x]bar[4][y] + bar[4][x]bar[7][y], {x, -1, 1}, {y, -1, 1},
    Contours -> {-0.01, 0.01}, PlotPoints -> 100,
    PlotLabel ->
      "(" <> ToString[7] <> "," <> ToString[4] <> ")plus(" <> ToString[4] <>
        "," <> ToString[7] <> ") freq=" <> ToString[freqsquare[7, 4]]] 

Out[3]= -ContourGraphics-

In[4]:= ContourPlot[(disc[3, 2][r] /. r -> Sqrt[u^2 + v^2] // Evaluate )Cos[
      3 Arg[u + I v]]If[u^2 + v^2 < 1, 1, 0], {u, -1, 1}, {v, -1, 1},
  Contours -> {-.01, .01}, PlotPoints -> 100,
  PlotLabel ->
    "(" <> ToString[3] <> "," <> ToString[2] <> ") freq=" <>
      ToString[freqdisc[3, 2]]]

Out[4]= -ContourGraphics-

*)

Begin["Chladni`Private`"] 

{a,m,b,n,x,c,nn,d,k,ka,y0,la,r};

Clear[a,m,b,n,x,c,nn,d,k,ka,y0,la,r];

bar[0] = Function[x, 1]
 
bar[1] = Function[x, x]
 
bar[m_Integer] := (bar[m] = If[OddQ[m], Function[x, 
        Evaluate[First[u[x, 0, b, 0, d, k] /. threecond] /. d -> 1 /. 
          k -> onevk[(m - 1)/2]]], Function[x, 
        Evaluate[First[u[x, a, 0, c, 0, k] /. threecondeven] /. c -> 1 /. 
          k -> evk[m/2]]]]) /; m > 1
 
u[x_, a_, b_, c_, d_, k_] = c*Cos[k*x] + a*Cosh[k*x] + d*Sin[k*x] + 
     b*Sinh[k*x]
 
threecond = {{a -> (d*Sin[k]*(-Sin[k] + Cos[k]*Tanh[k]))/
        (Cosh[k]*Sin[k] + Cos[k]*Sinh[k]), b -> d*Cos[k]*Sech[k], 
      c -> (d*Sinh[k]*(Sin[k] - Cos[k]*Tanh[k]))/(Cosh[k]*Sin[k] + 
         Cos[k]*Sinh[k])}}
 
onevk[0] = 0
 
onevk[m_Integer] := (onevk[m] = k /. FindRoot[trans[k] == 0, 
        {k, onevk[m - 1] + Pi}]) /; m >= 1
 
trans[k_] = -2*Sin[k] + 2*Cos[k]*Tanh[k]
 
threecondeven = {{a -> -(c*Csch[k]*Sin[k]), 
      b -> -((c*Cos[k]*(Cos[k] + Coth[k]*Sin[k]))/(-(Cosh[k]*Sin[k]) + 
          Cos[k]*Sinh[k])), d -> (c*Cosh[k]*(Cos[k] + Coth[k]*Sin[k]))/
        (Cosh[k]*Sin[k] - Cos[k]*Sinh[k])}}
 
evk[0] = -1
 
evk[m_Integer] := (evk[m] = k /. FindRoot[transeven[k] == 0, 
        {k, evk[m - 1] + Pi}]) /; m >= 1
 
transeven[k_] = -2*(Cos[k] + Coth[k]*Sin[k])
 
freqsquare[m_, n_] = freqbar[m] + freqbar[n]
 
freqbar[0] = 0
 
freqbar[1] = 0
 
freqbar[m_Integer] := If[OddQ[m], onevk[(m - 1)/2]^2, evk[m/2]^2] /; m > 1
 

freqdisc[m_Integer, 1] := (freqdisc[m, 1] = 
      (ka /. FindRoot[err[ka, m], {ka, (38*m/100 + 1)*Pi}])^2) /; 
     Inequality[0, LessEqual, m, Less, 11]
 
freqdisc[m_Integer, 1] := (freqdisc[m, 1] = 
      (ka /. FindRoot[err[ka, m], {ka, 1 + Sqrt[freqdisc[m - 1, 1]]}])^
       2) /; m > 10
 
freqdisc[m_Integer, 2] := (freqdisc[m, 2] = 
      (ka /. FindRoot[err[ka, m], {ka, 1 + Sqrt[freqdisc[m - 1, 2]]}])^
       2) /; m > 5
 
freqdisc[m_Integer, n_Integer] := (freqdisc[m, n] = 
      (ka /. FindRoot[err[ka, m], {ka, Sqrt[freqdisc[m, n - 1]] + Pi}])^
       2) /; m >= 0 && n > 0
 
err[ka_, 0] = ka*BesselJ[1, ka] - (ka^2*(BesselJ[0, ka] - BesselJ[2, ka]))/
      2 - (ka*BesselI[1, ka]*(ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - 
        ka*BesselJ[2, ka]))/(ka*BesselI[0, ka] + 2*BesselI[1, ka] + 
       ka*BesselI[2, ka]) + (ka^2*(BesselI[0, ka] + BesselI[2, ka])*
       (ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - ka*BesselJ[2, ka]))/
      (2*(ka*BesselI[0, ka] + 2*BesselI[1, ka] + ka*BesselI[2, ka])) + 
     (ka^2*(ka*BesselI[1, ka] + (ka*(BesselI[1, ka] + BesselI[3, ka]))/2)*
       (ka*BesselJ[0, ka] + 2*BesselJ[1, ka] - ka*BesselJ[2, ka]))/
      (2*(ka*BesselI[0, ka] + 2*BesselI[1, ka] + ka*BesselI[2, ka])) - 
     (ka^2*(-(ka*BesselJ[1, ka]) - (ka*(BesselJ[1, ka] - BesselJ[3, ka]))/2))/
      2
 
err[ka_, nn_] = -(ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 + 
     (ka*((ka*(BesselJ[-2 + nn, ka] - BesselJ[nn, ka]))/2 - 
        (ka*(BesselJ[nn, ka] - BesselJ[2 + nn, ka]))/2))/2 - 
     (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka])*
       (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 
        2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 
        2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/
      (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 
        2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 
        2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + 
     (ka*((ka*(BesselI[-2 + nn, ka] + BesselI[nn, ka]))/2 + 
        (ka*(BesselI[nn, ka] + BesselI[2 + nn, ka]))/2)*
       (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 
        2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 
        2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/
      (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 
        2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 
        2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + 
     (ka*((ka*((ka*(BesselI[-3 + nn, ka] + BesselI[-1 + nn, ka]))/2 + 
           (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka]))/2))/2 + 
        (ka*((ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka]))/2 + 
           (ka*(BesselI[1 + nn, ka] + BesselI[3 + nn, ka]))/2))/2)*
       (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 
        2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 
        2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/
      (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 
        2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 
        2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka])) + 
     2*nn^2*(BesselJ[nn, ka] + (BesselI[nn, ka]*(ka^2*BesselJ[-2 + nn, ka] + 
          2*ka*BesselJ[-1 + nn, ka] - 2*ka^2*BesselJ[nn, ka] - 
          4*nn^2*BesselJ[nn, ka] - 2*ka*BesselJ[1 + nn, ka] + 
          ka^2*BesselJ[2 + nn, ka]))/(-(ka^2*BesselI[-2 + nn, ka]) - 
         2*ka*BesselI[-1 + nn, ka] - 2*ka^2*BesselI[nn, ka] + 
         4*nn^2*BesselI[nn, ka] - 2*ka*BesselI[1 + nn, ka] - 
         ka^2*BesselI[2 + nn, ka])) - 
     nn^2*((ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 + 
       (ka*(BesselI[-1 + nn, ka] + BesselI[1 + nn, ka])*
         (ka^2*BesselJ[-2 + nn, ka] + 2*ka*BesselJ[-1 + nn, ka] - 
          2*ka^2*BesselJ[nn, ka] - 4*nn^2*BesselJ[nn, ka] - 
          2*ka*BesselJ[1 + nn, ka] + ka^2*BesselJ[2 + nn, ka]))/
        (2*(-(ka^2*BesselI[-2 + nn, ka]) - 2*ka*BesselI[-1 + nn, ka] - 
          2*ka^2*BesselI[nn, ka] + 4*nn^2*BesselI[nn, ka] - 
          2*ka*BesselI[1 + nn, ka] - ka^2*BesselI[2 + nn, ka]))) + 
     (ka*((ka*((ka*(BesselJ[-3 + nn, ka] - BesselJ[-1 + nn, ka]))/2 - 
           (ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2))/2 - 
        (ka*((ka*(BesselJ[-1 + nn, ka] - BesselJ[1 + nn, ka]))/2 - 
           (ka*(BesselJ[1 + nn, ka] - BesselJ[3 + nn, ka]))/2))/2))/2
 
disc[d_Integer, m_Integer] := (disc[d, m] = Function[x,(la = freqdisc[d, m]^2; 
       If[d == 0, y0[x], y[x] /. nn -> d])//Evaluate]) /; d >= 0 && m > 0
 
y0[r_] = BesselJ[0, la^(1/4)*r] + (BesselI[0, la^(1/4)*r]*
       (la^(1/4)*BesselJ[0, la^(1/4)] + 2*BesselJ[1, la^(1/4)] - 
        la^(1/4)*BesselJ[2, la^(1/4)]))/(la^(1/4)*BesselI[0, la^(1/4)] + 
       2*BesselI[1, la^(1/4)] + la^(1/4)*BesselI[2, la^(1/4)])
 
y[r_, 0, a_, b_, la_] = a*BesselI[0, la^(1/4)*r] + b*BesselJ[0, la^(1/4)*r]
 
y[r_] = BesselJ[nn, la^(1/4)*r] + (BesselI[nn, la^(1/4)*r]*
       (Sqrt[la]*BesselJ[-2 + nn, la^(1/4)] + 2*la^(1/4)*
         BesselJ[-1 + nn, la^(1/4)] - 2*Sqrt[la]*BesselJ[nn, la^(1/4)] - 
        4*nn^2*BesselJ[nn, la^(1/4)] - 2*la^(1/4)*BesselJ[1 + nn, la^(1/4)] + 
        Sqrt[la]*BesselJ[2 + nn, la^(1/4)]))/
      (-(Sqrt[la]*BesselI[-2 + nn, la^(1/4)]) - 
       2*la^(1/4)*BesselI[-1 + nn, la^(1/4)] - 
       2*Sqrt[la]*BesselI[nn, la^(1/4)] + 4*nn^2*BesselI[nn, la^(1/4)] - 
       2*la^(1/4)*BesselI[1 + nn, la^(1/4)] - 
       Sqrt[la]*BesselI[2 + nn, la^(1/4)])
 
y[r_, n_, a_, b_, c_, d_, la_] = a*BesselI[-n, la^(1/4)*r] + 
     c*BesselI[n, la^(1/4)*r] + b*BesselJ[-n, la^(1/4)*r] + 
     d*BesselJ[n, la^(1/4)*r]

End[]

EndPackage[]

Null
