(* Dispersion relation and wave diagram tools.                            *)

(*$Log: disperse.m%v $
# Revision 1.3  1994/11/20  19:50:16  bill
# added duckwaves in parametric form
#
# Revision 1.2  1994/11/19  14:49:58  bill
# removed examples to comments
#
# Revision 1.1  1994/11/18  16:40:38  bill
# Initial revision
#
#*)

(**************************************************************************)
(*                                                                        *)
(*  To calculate the group velocity for deep water waves:                 *)
(*                                                                        *)
(*  In[] :=    Vgroup[ w^2 - k, {k,w}]                                    *)
(*    output in the same order as the wavenumbers                         *)
(*                                                                        *)
(*  Some standard dispersion relations are predefined below.              *)
(*                                                                        *)
(**************************************************************************)
(*                                                                        *)
(*  Subs[ function, datalist, args] returns a list of the same length     *)
(*    as the datalist, each element is the function with specified        *)
(*    args substituted according to the data:                             *)
(*                                                                        *)
(*  Subs[ a+b, {1,2}, a]  ->gt  { 1+b, 2+b }                                *)
(*                                                                        *)
(*    datalist must be a list of type args                                *)
(*      perhaps having only a single element:                             *)
(*                                                                        *)
(*  Subs[ a+b, {{1,2}}, {a,b}] ->gt {3}                                     *)
(*                                                                        *)
(*  Subs[ a+b, {1,2}, {a,b}] ->gt INCORRECT ANSWER!                         *)
(*                                                                        *)
(*    function can itself be a list:                                      *)
(*                                                                        *)
(*  Subs[ {2 a, b}, {{1,2}}, {a,b}] ->gt {{2, 2}}                           *)
(*                                                                        *)
(**************************************************************************)


Print["Dispersive Wave Diagram Tools -- William L. Burke"]
Print["$Header: d:/math/bb/waves/disperse.m%v 1.3 1994/11/20 19:50:16 bill Exp bill $"]

Vgroup[ disprel_, pts_List:{k,w} ] := Module[ {vraw},
  vraw = Map[ D[ disprel, #]&, pts ];
  vraw / Apply[Plus, vraw*pts]
  ]

Subs[ g_, d_, pt_ ] := (g /. Thread[Rule[pt, #]])& /@ d

SetOptions[ListPlot, PlotJoined ->gt True];

deepwater = w^2 - k;

tanhwater = w^2 - k Tanh[k L];

capillary = w^2 - k^3;

(**************************************************************************)
(*  Waves with gravity and surface tension:                               *)
(*    one space dimension, one time dispersion relation: duckwaves        *)
(*    two spatial dimensions:  duckxy                                     *)
(*    parametric solution for {k,l,w}: duckpar                            *)
(*    waves in comoving frame, velocity unity: coduck                     *)
(*    wavediagram using parametric form: coduckvgp                        *)
(*    waves in comoving frame, velocity unity: coduck                     *)
(*  To generate a wave diagram:                                           *)
(*    T=0.01; ListPlot[Subs[coduckvgp, Table[i/100, {i, 110, 400}], lm]]  *)

duckwaves = w^2 - k (1 + T k^2);
duckxy    = duckwaves /. k ->gt Sqrt[k^2 + l^2]
duckpar   = { Sqrt[ lm(1 + T lm^2) ], 
              Sqrt[ lm^2 - lm (1 + T lm^2) ], 
              Sqrt[ lm(1 + T lm^2) ]}
coduck    = duckxy /. w ->gt k
coduckvgp =
  {((1 - 2*lm + 3*lm^2*T)*(lm + lm^3*T)^(1/2))/(-lm^2 + lm^4*T), 
    (lm^(1/2)*(-1 + lm - lm^2*T)^(1/2)*(1 + 3*lm^2*T))/(-lm^2 + lm^4*T)}

(*                                                                        *)
(**************************************************************************)

totalwater = w^2 - k (1 + T k^2) Tanh[k L];

(* Uses Plot to choose spacings for plottig. Expects rule ss for py *
 *   and vgroup defined                                             *
 *   Uses full range of px >gt 0                                      *)

AutoPlot := (
  tmp = Plot[(1./py /. N[ss]) /. px->gtu/(1-u),
    {u, 0, 1}, DisplayFunction ->gt Identity];
  udata = Nest[First, tmp, 4] // Transpose // First;
  pxdata = #/(1.-#)& /@ udata;
  pxydata = Subs[{px, (py /. N[ss])}, pxdata, {px}];
  gpdata = Subs[ N[vgroup], pxydata, {px, py}];
  gpxdata = Select[gpdata, Part[#, 2] gt $DisplayFunction];
)