Gaussian vertical profile eddy


In [1]:
%matplotlib inline
%load_ext secnum
%secnum

import numpy as np
import scipy as sp

import sys
sys.path.append('/media/data/Work/tools/python/dc/')
from dcutils import *

from IPython import display

# matplotlib section
import matplotlib as mpl
import matplotlib.pyplot as plt
# big figures
mpl.rcParams['savefig.dpi'] = 2 * mpl.rcParams['savefig.dpi']
# better color cycling
# import brewer2mpl
# colors =  brewer2mpl.get_map('Set2', 'qualitative', 8).mpl_colors
# mpl.rcParams['axes.color_cycle'] = colors


Horizontal Structure


In [16]:
x = np.linspace(-5,5,50)
y = np.linspace(-5,5,50)

r = np.hypot(x,y)
th = np.arctan2(y,x)

R = 1;
v0 = -1;

v = v0 * (r/R) * np.exp(-(r/R)**2) * np.cos(th)

plt.figure;
plt.plot(x,v)
plt.hold(True)
plt.plot(x, v/np.cos(th))


Out[16]:
[<matplotlib.lines.Line2D at 0x7f6d668f0d90>]

Vertical Structure

The modes for $(u,v,p)$ are $$ F_k = A_k \cos\left[ \frac{k \pi z}{Z}\right]$$

Thus, we need the cosine transform of the Gaussian profile = Gaussian again.

$$ \Sigma_k A_k \cos\left[ \frac{k \pi z}{Z}\right] = U_0 e^{-(z/H)^2}$$

In [121]:
# define vertical profile
zvec = np.linspace(-1200,0,200)
Z = zvec.min()
Hscale = 500
zprof = np.exp((zvec/Hscale)**1)

# baroclinic modes for (u,v,p)
mode1 = np.cos(np.pi/Z * zvec)
mode2 = np.cos(2*np.pi/Z * zvec)

plt.plot(zprof, zvec)
plt.hold(True)
plt.plot(mode1, zvec)
plt.plot(mode2, zvec)
plt.legend(['Gaussian', 'mode1', 'mode2'], 'upper left')

# determine amplitudes
modek = xrange(10)

A = np.zeros([len(modek), 1])
B = np.zeros([len(modek), 1])
modesum = np.zeros_like(zprof)

for k in modek:
    index = k;
    cmode = np.cos(k * np.pi/Z * zvec)
    smode = np.sin(k * np.pi/Z * zvec)
    A[index] = 2/Z * np.trapz(zvec, zprof * cmode)
    B[index] = 0 #2/Z * np.trapz(zvec, zprof * smode)
    
    if index == 0:
        A[0] = A[0]/2
        
    modesum += A[index] * cmode + B[index] * smode
    
norm = 1 #np.sum(zprof) / np.sum(modesum)
plt.figure()
plt.plot(modek, A)
plt.hold(True)
plt.plot(modek, B)

plt.figure()
plt.hold(True)
plt.plot(zprof, zvec)
plt.plot(norm*modesum, zvec)


Out[121]:
[<matplotlib.lines.Line2D at 0x1568de90>]

In [39]:



Out[39]:
array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
       -0., -0., -0., -0.,  0.])