holterma
/
analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

initial import

master
Lars Umlauf 3 years ago
commit
12abaa520c
11 changed files with 964 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 8
      iow_fancypsd/README
  2. BIN
      iow_fancypsd/data.mat
  3. 165
      iow_fancypsd/iow_fancypsd.m
  4. 125
      iow_fancypsd/iow_fancypsd_test.m
  5. 44
      iow_fitellipse/iow_fitellipse.m
  6. 63
      iow_fitellipse/iow_fitellipse_test.m
  7. 120
      iow_fourierfilter/iow_fourierfilter.m
  8. 51
      iow_fourierfilter/iow_fourierfilter_test.m
  9. 7
      iow_iwmodes/README
  10. 182
      iow_iwmodes/iow_iwmodes.m
  11. 199
      iow_iwmodes/iow_iwmodes_test.m

8
iow_fancypsd/README
Unescape View File

@ -0,0 +1,8 @@
This is a MATLAB tool for spectral analysis including variable
confidence intervals. Complex stochastic veriables will be analyzed
with respect to clockwise and anti-clockwise rotary components.
Run "iow_fancypsd_test.m" for a quick test.
Info: lars.umlauf@io-warnemuende.de

BIN
iow_fancypsd/data.mat
View File

Binary file not shown.

165
iow_fancypsd/iow_fancypsd.m
Unescape View File

@ -0,0 +1,165 @@
function [F,P,Ctop,Cbot] = iow_fancypsd(x,NFFTmax,Ffac,Fmin,Fs,p,pPlot)
% This routine computes a series of spectral densities with increasingly
% smaller segment length, NFFT, according to the Welch algorithm with 50%
% overlap. The maximum segment length used in the first iteration is
% NFFTmax. For each iteration, the segment length is divided by two, and a
% new spectrum is computed. The spectrum with the largest segment length,
% NFFTmax, is used for the frequency windo from 0 to Fmin. The spectrum
% from the next interation is used for the frequency range from Fmin to
% Ffac*Fmin, and so on until the Nyquist frequency is reached. Due to the
% increasingly smaller segment length, the confidence intervals become
% stepwise smaller towards higher frequencies.
% Input arguments:
% x: Data vector. If x is complex, rotary spectra will be computed
% NFFTMax: Maximum segment length used in first iteration
% Ffac: Frequency window shifts by Ffac for every iteration
% Fmin: Lowest frequency window ranges from 0 to Fmin
% Fs: Sampling frequency
% p: Confidence interval (e.g. p=0.95)
% pPlot: Center confidence intervals around this value (for plotting)
% Output arguments:
% F: Frequency vector
% P: Power spectral density. If x is complex, also P is complex:
% real(P) contains the counter-clockwise rotary spectrum
% imag(P) contains the clockwise rotary spectrum
% Ctop: Upper confidence level (with respect to 'pPlot')
% Cbot: Lower confidence level (with respect to 'pPlot')
%
% Info: lars.umlauf@io-warnemuende.de
%% initialize
% Inititalize output arrays
F = [];
P = [];
Cbot = [];
Ctop = [];
% Inititalize spectrum parameters
NFFT = NFFTmax; % segment and window width
Fl = 0; % lower range for frequency window
Fu = Fmin; % upper range for frequency window
%% iterate until Nyquist frequency is exceeded
it = 0;
while( Fl < 0.5*Fs )
% count interation
it = it + 1;
% report current parameters
disp(['Iteration: ' num2str(it)] );
disp(['Freq. window: ' num2str(Fl,'%2.2e') ' - ' ...
num2str(Fu,'%2.2e')]);
disp(['NFFT: ' num2str(NFFT,'%6.0f')]);
disp(['Freq. resolution: ' num2str(Fs/NFFT,'%2.2e')]);
disp(' ');
% compute spectrum for current segment length
OverlapPercent = 50;
Hs = spectrum.welch('Hamming',NFFT,OverlapPercent);
if isreal(x)
hpsd = psd(Hs,x,'NFFT',NFFT,'Fs',Fs, ...
'ConfLevel',p,'SpectrumType','onesided');
isRotary = 0;
else
hpsd = psd(Hs,x,'NFFT',NFFT,'Fs',Fs, ...
'ConfLevel',p,'SpectrumType','twosided');
isRotary = 1;
end
% save current iteration
Fit = hpsd.Frequencies;
Pit = hpsd.Data;
Cit = hpsd.ConfInterval;
% construct rotary spectra from left and right hand sides of complex
% spectra
if isRotary
% Nup is the highest frequency that can be represented
Nup = NFFT/2;
if round(Nup) == Nup
% even
Nup = NFFT/2 + 1;
even = 1;
else
% odd
Nup = (NFFT+1)/2;
even = 0;
end
% interprete right hand side of spectrum as counter-clockwise
Fit = Fit(1:Nup);
Pplus = Pit(1:Nup);
% interprete left hand side of spectrum as clockwise
% Matlab spectra have be flipped
Pminus = zeros(Nup,1);
if even
Pminus(2:Nup-1) = Pit(NFFT:-1:Nup+1);
else
Pminus(2:Nup) = Pit(NFFT:-1:Nup+1);
end
% compose complex spectrum
Pit = Pplus + i*Pminus;
% just use the right hand side of confidence intervals
% (left hand side will be equal)
Cit = Cit(1:Nup,:);
end
% cut out the frequency window of this iteration
Ind = find(Fit > Fl & Fit < Fu);
Fit = Fit(Ind);
Pit = Pit(Ind);
Cit = Cit(Ind,:);
% convert Matlab confidence intervals to version that is constant on a
% log-log plot
[Cbotit,Ctopit] = confLabel(real(Pit),Cit,pPlot);
% append current frequency range to output
F = [F ; Fit];
P = [P ; Pit];
Cbot = [Cbot ; Cbotit];
Ctop = [Ctop ; Ctopit];
% update frequency range and segment length
Fl = Fu;
Fu = Ffac*Fu;
NFFT = floor(NFFT/2);
end
end
function [Cbot Ctop] = confLabel(Psd,ConfInterval,CentralValue)
% This function converts from the Matlab way to plot confidence intervals
% to confidence intervals that do not depend on frequency if plotted
% log-log scale.
DCtop = log10(ConfInterval(:,2)) - log10(Psd);
DCbot = log10(Psd) - log10(ConfInterval(:,1));
Ctop = 10.^(log10(CentralValue) + DCtop);
Cbot = 10.^(log10(CentralValue) - DCbot);
end

125
iow_fancypsd/iow_fancypsd_test.m
Unescape View File

@ -0,0 +1,125 @@
%IOW_FANCYPSD_TEST Test IOW spectral analysis tool
% This routine is supposed to help understanding how iow_fancypsd() works,
% and what kind of arguments it expects.
%
% Info: lars.umlauf@io-warnemuende.de
%% load sample data
% This contains velocities from the Baltic Sea over a couple of months
% The inertial frequency is about f=1.2E-4 (about 0.069 cph).
data = load('data');
% date and time
date = data.date; % MATLAB date
days = (date - date(1));
hours = days*24;
% velocity components
u = data.u/100; % convert to m/s
v = data.v/100;
w = u + i*v; % complex velocity
% sampling parameters
N = length(date);
dt = mean(diff(hours)); % sampling period (in hours)
Fs = 1/dt; % sampling frequency (in cph)
% some plotting stuff
XULim = [days(1) days(end)];
ULim = [-0.15 0.15];
XPLim = [5e-4 1e0];
PLim = [2e-6 2e0];
%% do spectral analysis
% detrend
u = detrend(u,'constant');
v = detrend(v,'constant');
% set spectrum parameters
NFFTmax = floor(N/4); % set largest segment length
% for low-frequency part of the spectrum
Fmin = 5e-3; % lowest frequency range from 0 to Fmin
Ffac = 5; % increase frequency range by this factor for
% each interation
p = 0.95; % set confidence interval
pPlot = 1e-5; % plot confidence intervals around this value
% compute u and v spectra
[F,Pu,Ctop,Cbot] = iow_fancypsd(u,NFFTmax,Ffac,Fmin,Fs,p,pPlot);
[F,Pv,Ctop,Cbot] = iow_fancypsd(v,NFFTmax,Ffac,Fmin,Fs,p,pPlot);
% compute rotary spectra (since w is complex this is done automatically
% by iow_fancypsd
[F,Pw,Ctop,Cbot] = iow_fancypsd(w,NFFTmax,Ffac,Fmin,Fs,p,pPlot);
%% plot velocity time series
figure(1);
clf;
subplot(2,1,1)
plot(days,u,'r')
set(gca,'FontSize',14,'XLim',XULim,'YLim',ULim,'XTickLabel',{});
ylabel('u (m s^{-1})');
title('Velocity Timeseries')
box on;
subplot(2,1,2)
plot(days,v,'b')
set(gca,'FontSize',14,'XLim',XULim,'YLim',ULim);
xlabel('t (days)');
ylabel('v (m s^{-1})');
box on;
%% plot spectra
figure(2);
clf;
% first the ordinary spectra of u and v
subplot(2,1,1)
line(F,Pu,'Color','r','LineWidth',2);
line(F,Pv,'Color','b','LineWidth',2);
line(F,Ctop,'Color','k');
line(F,Cbot,'Color','k');
legend('u','v',[num2str(100*p) '%'])
set(gca,'FontSize',14,'XScale','Log','YScale','Log')
set(gca,'XLim',XPLim,'YLim',PLim,'XTickLabel',{})
ylabel('P (m^2 s^{-2}) cph^{-1}');
title('Power Spectral Densities')
box on;
% second the rotary components
subplot(2,1,2)
line(F,imag(Pw),'Color','r','LineWidth',2);
line(F,real(Pw),'Color','b','LineWidth',2);
line(F,Ctop,'Color','k');
line(F,Cbot,'Color','k');
legend('clockwise','anti-clockwise',[num2str(100*p) '%'])
set(gca,'FontSize',14,'XScale','Log','YScale','Log')
set(gca,'XLim',XPLim,'YLim',PLim)
xlabel('f (cph)')
ylabel('P (m^2 s^{-2}) cph^{-1}');
box on;

44
iow_fitellipse/iow_fitellipse.m
Unescape View File

@ -0,0 +1,44 @@
function [uAmp,vAmp,alpha,wFit] = iow_fitellipse(t,w,omega)
% IOW_FITELLIPSE Computes elliptical fit and fit parameters from the
% complex input signal w=u+i*v (corresponding, for example, to two velocity
% components. The theory is described in the section about rotatary
% spectra in the notes.
%
% Input arguments:
% t time vector of size(N,1)
% w complex signal of size(N,1)
% omega frequency of fitted siginal (rad/s)
%
% Output arguments:
% uAmp Major axis amplitude
% vAmp Minor axis amplitude
% alpha Rotation angle of major axis with respect to u
% wFit Fitted signal
%% get Fourier coefficients and compute parameters
N = length(t);
wMean = mean(w);
wHatPlus = 1/N* sum(w .* exp(-1i*omega*t) );
wHatMinus = 1/N* sum(w .* exp( 1i*omega*t) );
APlus = abs(wHatPlus);
PhiPlus = angle(wHatPlus);
AMinus = abs(wHatMinus);
PhiMinus = angle(wHatMinus);
alpha = (PhiPlus + PhiMinus)/2;
uAmp = APlus + AMinus;
vAmp = APlus - AMinus;
%% reconstruct signal
wFit = wHatPlus * exp(1i*omega*t) ...
+ wHatMinus * exp(-1i*omega*t) + wMean;

63
iow_fitellipse/iow_fitellipse_test.m
Unescape View File

@ -0,0 +1,63 @@
% This routine demonstrates how current data (e.g. from tidal flows or
% near-inertial wave motions can by analyzed with respect to their
% elliptical character (major axis, rotation).
%% Generate tidal data with some noise
T = 4; % duration of signal
dt = 0.001; % time resolution
omega = 2*pi; % frequency of signal (rad/s)
phi = pi/4; % rotation of major axis
e = 1.2; % ellipticty (ratio of major to minor axis)
% construct data
t = [0:dt:T]'; % time vector
N = length(t);
u = e*cos(omega*t)+(rand(N,1)-0.5)*0.3;
v = sin(omega*t)+(rand(N,1)-0.5)*0.3;
% rotate current ellipse
w = u+1i*v;
w = abs(w).*exp(1i*(angle(w)+phi));
u = real(w);
v = imag(w);
%% get fit parameters and report
[uAmp,vAmp,alpha,wFit] = iow_fitellipse(t,w,omega);
disp('Fit parameters')
disp(['Ellipticity = ' num2str(uAmp/vAmp)]);
disp(['Rotation (rad) = ' num2str(alpha)]);
disp(' ')
%% do the plotting
figure(1);
clf;
line(t,real(w),'LineWidth',2,'Color','k')
line(t,real(wFit),'LineWidth',2,'Color','r')
line(t,imag(w),'LineWidth',1,'Color','k')
line(t,imag(wFit),'LineWidth',1,'Color','r')
legend('u','uFit','v','vFit')
set(gca,'FontSize',14)
xlabel('t (s)')
ylabel('u (m/s)')
box on;
figure(2);
clf;
line(real(w),imag(w),'LineWidth',2,'Color','k')
line(real(wFit),imag(wFit),'LineWidth',2,'Color','r')
axis equal;
box on;
legend('w','wFit')
set(gca,'FontSize',14)
xlabel('u (m/s)')
ylabel('v (m/s)')

120
iow_fourierfilter/iow_fourierfilter.m
Unescape View File

@ -0,0 +1,120 @@
function uFilt = iow_fourierfilter(u,NDelta,FilterType,DoPlot)
% FOURIERFILTER Filters the input signal u according to a given
% filter window with NDelta points. For argument FilterType, there is
% currently the choice between 'box' and 'bartlett' available. For
% DoPlot=1 some plotting of the signal and filter window (FFT and
% physical space) is done.
% For filtering, the Fourier coefficients are multiplied by a variable
% transfer function, and the signal is transformed back to physical space.
%% check input signal
[N,M] = size(u);
if M~=1
error('Input signal u is not of size(N,1)')
end
%% compute number of represented frequencies
% index of hightest representable frequency in two-sided spectra
% This depends on the way fft() deals with even and odd record numbers
Nup = N/2;
if round(Nup) == Nup
% even
Nup = N/2 + 1;
m = [0:Nup-1]';
m = [-flipud(m) ; m(2:end-1)];
even = 1;
else
% odd
Nup = (N+1)/2;
m = [0:(Nup-1)]';
m = [-flipud(m) ; m(2:end)];
even = 0;
end
%% create filtering window
wD = 2*pi*m/N*NDelta;
switch lower(FilterType)
case 'box'
gHat = sin(wD/2) ./ (wD/2);
case 'bartlett'
gHat = ( sin(wD/4) ./ (wD/4) ).^2;
otherwise
error('unkown FilterType')
end
gHat(Nup) = 1;
%% filtering
% everything rescaled by 1/N and shifted to the notation used in the notes
% compute FFT of the signal
uHat = 1/N*fft(u);
uHat = fftshift(uHat);
% spectrum is filtered by multiplying with FFT of window function
uFiltHat = uHat .* gHat;
% compute inverse FFT to get the filtered signal in physical space
g = ifft(N*ifftshift(gHat));
uFilt = ifft(N*ifftshift(uFiltHat));
%% plotting
if DoPlot
FontSize = 14;
nFig = 99;
% Fourier coefficients
nFig = nFig + 1;
figure(nFig);
clf;
subplot(4,1,1)
plot(m,real(uHat),'o-k',m,real(uFiltHat),'o-r')
grid on;
set(gca,'FontSize',FontSize,'XTickLabel',{})
ylabel('$\hat{u}_m$','Interpreter','latex')
title('FFT of signal (real part)');
legend('raw','filtered')
subplot(4,1,2)
plot(m,imag(uHat),'o-k',m,imag(uFiltHat),'o-r');
grid on;
set(gca,'FontSize',FontSize,'XTickLabel',{})
ylabel('$\hat{u}_m$','Interpreter','latex')
title('FFT of signal (imaginary part)')
legend('raw','filtered')
subplot(4,1,3)
plot(m,gHat,'o-k');
grid on;
set(gca,'FontSize',FontSize)
ylabel('$\hat{g}_m$','Interpreter','latex')
xlabel('$m$','Interpreter','Latex')
title('FFT of filter window')
subplot(4,1,4)
plot([0:N-1]',g,'o-k');
grid on;
set(gca,'FontSize',FontSize)
ylabel('$g_i$','Interpreter','latex')
xlabel('$i$','Interpreter','Latex')
title('Filter window')
end

51
iow_fourierfilter/iow_fourierfilter_test.m
Unescape View File

@ -0,0 +1,51 @@
% illustrates different kinds of Fourier filters
% using the routine iow_fourierfilter
%% define parameters
N = 20; % number of sampling intervals
n = 2; % mode number of input signal
% filter specification
FilterType = 'bartlett'; % 'box' or 'bartlett'
NDelta = 8; % width (points) of filtering window
% other stuff
DoPlot = 1; % do plotting of spectra etc.
FontSize = 14;
%% create some sinusoidal data
% compose signal index vector
% (indexing follows notes with 0 as lowest index)
t = [0:N-1]';
% generate signal
u = cos(2*pi*n/N*t);
% compute filtered signal
uFilt = iow_fourierfilter(u,NDelta,FilterType,DoPlot);
%% plotting
nFig = 0;
% time series
nFig = nFig + 1;
figure(nFig);
clf;
plot(t,u,'o-k',t,uFilt,'o-r')
grid on;
legend('raw','filtered')
set(gca,'FontSize',FontSize)
xlabel('$i$','Interpreter','Latex')
ylabel('$u_i$','Interpreter','Latex')
title('Raw and filtered signals')

7
iow_iwmodes/README
Unescape View File

@ -0,0 +1,7 @@
This is a MATLAB script computing the vertical eigenmodes of internal
waves in stratified basins with flat bottom.
Run "iow_iwmodes_test.m" to see an example.
Info: lars.umlauf@io-warnemuende.de

182
iow_iwmodes/iow_iwmodes.m
Unescape View File

@ -0,0 +1,182 @@
function [W,Psi,dPsidz,hW,c,nm] = iow_iwmodes(NN,h,omega,f,m,isRot,isNH)
% IW_MODES computes the vertical structure and speed of an internal wave
% with frequency omega from the equation W'' + N^2 /c^2 F W = 0. Here, W
% is the vertical strucuture of the vertical velocity, c the propagation
% speed, N the buoyancy frequency and F = (1 - rN) / (1 - rf), where
% rN = (omega/N)^2 and rf = (f/omega^2) with f denoting the rotation rate.
%
% The returned eigenvectors are checked for physical plausibility. They may
% though still contain numerical modes that have no physical significance.
%
% INPUT ARGUMENTS:
% NN: square of buoyancy frequency given at n vertical levels
% h: grid spacing between vertical levels (n-1 values)
% (note that the grid spacing may be variable)
% omega: frequency of internal wave
% (only used if isRot=1 or isNH=1, see below)
% f: Coriolis parameter
% m: return the m largest eigenvalues and eigenvectors
% isRot: False, if rotation effects are ignored (rf=0)
% is NH: False, if non-hydrostatic effects are ignored (rN=0)
%
% OUTPUT ARGUMENTS:
% W: vertial eigenvectors of w, size(n,m)
% Psi: vertial eigenvectors of u,v,etc, size(n-1,m)
% dPsidz: derivative of Psi, size(n,m)
% hW: distance between Psi-levels (at W-levels), size(n,m)
% c: propagation speeds (vector of size(1,m)
% nm: number of valid eigenvalues found (smaller than m)
%
% Info: lars.umlauf@io-warnemuende.de
%% do some consistency checks
n = length(NN);
if ~isvector(NN) && ~isvector(h)
error('Input arguments NN and h have to be vectors')
end
if length(NN) ~= length(h)+1
error('Size of input arguments NN and h inconsistent');
end
%% add effects of rotation
if isRot
rf = (f/omega)^2;
if (rf < 0) || (rf > 1)
error('Omega must be in the range 0...f');
end
else
rf = 0;
end
%% add non-hydroststic effects
if isNH
rN = omega^2 ./ NN;
else
rN = zeros(size(NN));
end
% compose F
F = (1-rN) ./ (1-rf);
%% compose discritization matrices
% allocate
A = zeros(n,n);
B = zeros(n,n);
% this is a simple trick to yield W(1)=W(n)=0
A(1,1) = 1;
A(n,n) = 1;
% populate A matrix
for i=2:n-1
a = 1/( h(i-1) + h(i) );
b = 1/h(i-1) + 1/h(i);
A(i,i-1) = 2*a/h(i-1);
A(i,i ) = -2*a*b;
A(i,i+1) = 2*a/h(i);
end
% populate B-matrix
for i=2:n-1
B(i,i) = F(i)*NN(i);
end
%% compute eigenvales and vectors
[E,L] = eig(A,B);
l = diag(L)';
%% check and sort in descending order
[W,c,nm] = iw_check_modes(E,l,m);
% compute Eigenfunction for horizontal velocities
Psi = nan(n-1,nm);
for i=1:n-1
Psi(i,:) = ( W(i+1,:) - W(i,:) ) / h(i) ;
end
% compute Eigenfunction for vertical shear
hW = nan(n,1); % distance between Psi-levels (at W-levels)
dPsidz = nan(n,nm); % shear at W-levels
hW(1) = 0; % thicknesses at boundaries undefined
hW(end) = 0;
dPsidz(1,:) = 0; % no shear at boundaries
dPsidz(end,:) = 0;
for i=2:n-1
hW(i) = 0.5*( h(i) + h(i-1) );
dPsidz(i,:) = ( Psi(i,:) - Psi(i-1,:) ) / hW(i) ;
end
end
function [W,c,nm] = iw_check_modes(E,l,m)
% This routine does some plausibility checks, and rejects eigenvalues and
% eigenvectors that fail.
%
% INPUT ARGUMENTS:
% E: matrix of size(n,n) with the eigenvectors in columns
% l: vector containing corresponding eigenvalues
% m: number of desired eigenvectors to be returned
%
% OUTPUT ARGUMENTS:
% W: matrix of size(n,m) with m returned eigenvectors
% c: vector containing corresponding phase speeds
% (largest speeds come first)
% nm: number of physically (probably) plausible
% eigenvectors
% check if all eigenvalues are real
if sum(~isreal(l)) ~= 0
error('Complex Eigenvalues detected');
end
% check if all eigenvalues are non-NaN
if sum(isnan(l)) ~= 0
error('NaN Eigenvalues detected');
end
% remove positive (unphysical) eigenvalues
Ind = find(l >= 0);
l( Ind) = [];
E(:,Ind) = [];
% sort with smallest eigenvalue first
[l,Ind] = sort(l,'descend');
E = E(:,Ind);
% compute phase speeds (largest first)
cRaw = 1 ./ sqrt(-l);
% return fields
c = nan(1,m);
W = nan(size(E,1),m);
nm = min(m,length(cRaw));
c(1:nm) = cRaw(1:nm);
W(:,1:nm) = E(:,1:nm);
end

199
iow_iwmodes/iow_iwmodes_test.m
Unescape View File

@ -0,0 +1,199 @@
% IOW_IW_MODES_TEST
% This function provides information about how to use the iow_iwmodes()
% computing internal waves modes in a stratified fluid.
%
% Info: lars.umlauf@io-warnemuende.de
%% set general parameters
isRot = 0; % set to 1 to include rotation
isNH = 1; % set to 1 to include non-hydrostatic effects
f = 1.0e-4; % set Coriolis parameter
H = 50; % water depth
n = 200; % number of vertical grid levels
m = 5; % number of modes to analyze
%% set grid parameters
GridType = 'const'; % vertical grid spacing is computed from:
% 'const': constant layer spacing
% 'sine': sinusoidal layer spacing
% 'adapt': spacing is finer at locations
% with large N
hfac = 0.1; % mimum grid size is "hfac" times maximum size
%% set stratifiction parameters
StratType = 'pycno'; % type of vertical stratification:
% 'const': constant N
% 'pycno': pycnocline a zP with
% half-width wP
% buoyancy frequency
NN1 = 1.0e-4; % used for constant stratification
NNP = 1.0e-3; % maximum value in pycnocline
% geometry
zP = -15; % z-position of thermocline
wP = 2; % half-width of thermocline
%% construct 'measured' stratification
nM = 1000; % number of 'measured' data points
NNM = nan(nM,1); % buoyancy frequency squared
zM = linspace(-H,0,nM)'; % vertical positions
switch lower(StratType)
case 'const'
NNM = NN1*ones(nM,1);
case 'pycno'
% Gaussian bump in NN for thermocline
NNM = NNP * exp( - 0.5*( (zM-zP)/wP ) .^2 );
otherwise
error(['StratType "' StratType '" is unsupported'])
end
%% construct grid
h = linspace(0,1,n-1)';
switch lower(GridType)
case 'const'
h = ones(n-1,1);
case 'sine'
% zooming towards surface and bottom
h = linspace(0,1,n-1)';
h = sin(pi*h);
h = max(h,hfac*max(h));
otherwise
error(['GridType "' GridType '" is unsupported'])
end
% stretch to fill complete water column
h = H/sum(h) * h;
% construct z-axis
z = nan(n,1);
z(1) = -H;
for i=2:n
z(i) = z(i-1)+h(i-1);
end
zp(1) = -H + 0.5*h(1);
for i=2:n-1
zp(i) = zp(i-1) + 0.5*( h(i-1) + h(i) );
end
% interpolated 'measured' values to grid
NN = interp1(zM,NNM,z);
% implement no-flux condition at surface and bottom (if applies)
NN(1) = 0;
NN(end) = 0;
%% plot stratification
figure(1);
clf;
plot(NN,z,'ko-')
set(gca,'FontSize',14,'XLim',[0 1.1*max(NN)],'YLim',[z(1) z(n)]);
xlabel('N^2 (s^{-2})');
ylabel('z (m)');
title(['Grid type is: ' GridType '. Stratification type is: ' StratType]);
%% get analytical solution
% pre-allocate
ca = nan(1,m); % propgation speed
Wa = nan(n,m); % vertical velocity structure of modes
% compute solution
haveAnalyt = 1;
switch StratType
case 'const'
N = sqrt(NN1);
for i=1:m
ca(i) = N*H/i/pi;
Wa(:,i)= (-1)^(i+1)*sin(i*pi*z/H);
end
otherwise
warning(['No analytical solution for StratType=' StratType]);
haveAnalyt = 0;
end
%% compute numerical solution
% put something reasonable here for isRot=1 or isNH=1
omega = 0.7*sqrt(NNP);
% get solution
[W,Psi,dPsidz,hW,c,nm] = iow_iwmodes(NN,h,omega,f,m,isRot,isNH);
%% plot results
% phase velocities
figure(2)
clf;
plot(c,'ok-')
set(gca,'FontSize',14);
xlabel('Mode');
ylabel('c (m s^{-1})');
title('Phase Speeds');
box on;
% mode structure
figure(3);
clf;
for i=1:nm
h1 = subplot(1,2,1);
cla(h1);
line( W(:,i),z,'Marker','o','Color','k','Parent',h1);
line(Wa(:,i),z,'Marker','o','Color','r','Parent',h1);
if haveAnalyt
legend('numerical','analytical')
end
set(h1,'FontSize',14,'YLim',[z(1) z(n)]);
title(['Mode-Number: ' num2str(i)]);
xlabel('W');
ylabel('z (m)');
box on;
h2 = subplot(1,2,2);
cla(h2);
line( Psi(:,i),zp,'Marker','o','Color','k','Parent',h2);
set(h2,'FontSize',14,'YTickLabel',{},'YLim',[z(1) z(n)])
xlabel('\Psi');
title(['Press key to proceed.']);
box on;
pause
end
Write Preview
Loading…
Cancel
Save