Internal Gravity Waves 
in an Atmosphere

This notebook was used to find the trapped wave solution shown in Figure 1 of "Anelastic Equations for the Atmosphere".  The equation we solve is:

((d^2)/(dz^2))(OverHat[w]) -(1/H)(d/dz)(OverHat[w]) -(k^2)(1-((N^2)/(ω^2)))(OverHat[w])=0

This notebook could benefit from more comments and explanations.  Basically, we are solving an ODE with boundary conditions of (OverHat[w])=0 at z=0 and z=lid , but solving it as an initial value problem.  This means we keep "guessing" or improving the eigenvalue ω until the shooting scheme hits the target of (OverHat[w])=0 at z=lid.

Functions for solving the ODE

soln[ω_,k_]:=NDSolve[{
			w''[z]-w'[z]/H-k^2*(1-BuoyFreq^2/ω^2)*w[z]==0,
			w[0]==0,
			w'[0]==1},
		w,{z,0,lid}]

wAtlid[ω_,k_]:=w[lid]/.soln[ω,k][[1]]

Some parameters we don't change much

H=10.;

lid=20;

Try it with constant BuoyFreq

BuoyFreq=.01;

trysoln=soln[.006,1]

ParametricPlot[{w[z]/.trysoln[[1]],z},{z,0,lid},AspectRatio->2]

Try shooting for w=0 at the top

The roots  in the following plot "hit the target" of w=0 at z=lid:

Plot[wAtlid[ω,1],{ω,0.006,0.007}]

Let's find one eigenvalue ω exactly:

aroot=FindRoot[wAtlid[ω,1],{ω,.006,.0064}]

frustration=Interpolation[Table[{ω,wAtlid[ω,1]},{ω,.006,.007,.00001}]]

Plot[frustration[ω],{ω,.006,.007}]

aroot=FindRoot[frustration[ω],{ω,.006,.0064}]

hitsoln=soln[ω/.aroot[[1]],1]

ParametricPlot[{w[z]/.hitsoln[[1]],z},{z,0,lid},AspectRatio->2]

Trapped Waves

BuoyFreq:=0.01*Sech[8*(z-lid/2)/lid]^2

buoyfreqplot=ParametricPlot[{BuoyFreq,z},{z,0,lid},AspectRatio->2,Ticks->{{0.005`,0.01`},Automatic},AxesLabel->{"N","z"}]

Export["buoyfreqplot.eps",buoyfreqplot,"EPS"];

wAtlid[.001,1]

Plot[wAtlid[ω,1],{ω,0.002`,0.008`}]

pain=Interpolation[Table[{ω,wAtlid[ω,1]},{ω,.003,.0035,.00001}]]

aroot=FindRoot[pain[ω],{ω,.003,.0035}]

hit=ω/.aroot

hitsoln=soln[hit,1]

wmax=Max[Table[w[z]/.hitsoln[[1]],{z,0,lid,.001*lid}]]

wplot=ParametricPlot[{((w[z])/wmax)/.hitsoln[[1]],z},{z,0,lid},AspectRatio->2,AxesLabel->{"w","z"}]

Export["wplot.eps",wplot,"EPS"];