mathsam

classdef SpectrumAnalysis < handle
    properties (SetAccess = public, GetAccess = public)
        EperWavenumber
        wavenumber
        energyContainWavenumber
        t_start
        t_end
        p_lowb
        p_upperb
        numFourier

mathsam

%-----calculate the pressure at the tropopause
% do lat=ncread, pfull=ncread first!
pfull_range=20; % levels of pressure in modelling
lat_point=49; % 
T_p=zonal_temp(:,lat_point);
d_T=T_p(1:pfull_range-2)-T_p(3:pfull_range);
d_pfull=pfull(1:pfull_range-2)-pfull(3:pfull_range);
dT_dz=-9.8/287*pfull(2:pfull_range-1).*d_T./d_pfull./T_p(2:pfull_range-1);
subplot(1,2,1)
plot(pfull(2:pfull_range-1),1000*dT_dz,'linewidth',2)

mathsam

classdef DynamicsAnalysis < handle
    properties (SetAccess = public, GetAccess = public)
        filename
        t_start = 300;
        t_end
        numLat
        latDeg % latitudes in degree
        pfull % pressure in mb
        refP % (1:2) don't use surface or top wind
        refLat % reference latitude for f and beta

mathsam

syms c epsilon k N gamma

A = [-(-c+0.5*(3+epsilon))*(1+k^2)+k*(1+gamma), -c+0.5*(3+epsilon), 0; ...
    N*(-c+0.5*(1+epsilon)), -(-c+0.5*(1+epsilon))*(k^2+2*N)+k*(gamma-N*(0.5-0.5*epsilon)), N*(-c+0.5*(1+epsilon)); ...
    0, -c*N/epsilon, c*(k^2+N/epsilon)+k*(gamma-0.5*N*(1+epsilon)/epsilon)];

det_A = vpa(det(A));

%%
clear

mathsam

syms c epsilon1 epsilon2 k N1 N2 gamma

A = [-(-c+0.5*(2*epsilon1+epsilon2+3))*(1+k^2)+1+gamma, -c+0.5*(2*epsilon1+epsilon2+3), 0, 0;...
    -c+0.5*(2*epsilon1+epsilon2+1), -(-c+0.5*(2*epsilon1+epsilon2+1))*(k^2+1+N1)+gamma-1+0.5*N1*(1+epsilon1), N1*(-c+0.5*(2*epsilon1+epsilon2+1)),0; ...
    0, N1/epsilon1*(-c+0.5*(epsilon1+epsilon2)), -(-c+0.5*(epsilon1+epsilon2))*(k^2+N1/epsilon1)-N2/epsilon1, N2/epsilon1;...
    0, 0, -c*N2/epsilon2, c*(k^2+N2/epsilon2)+gamma-0.5*N2/epsilon2*(epsilon1+epsilon2)];

det_A = vpa(det(A));

%%

mathsam

function ave_Y = plot_ensemble(gamma,Y)
% input: gamma(N); Y(m,N)
% plot average shading for the N experiments
% return the mean

alpha_confidence = 0.1;              % significance level
if length(gamma) ~= size(Y,2)
    error('Y has the wrong size');
    return
end

mathsam

%PRINT_PDF  Prints cropped figures to pdf with fonts embedded
%
% Examples:
%   print_pdf filename
%   print_pdf(filename, fig_handle)
%
% This function saves a figure as a pdf nicely, without the need to specify
% multiple options. It improves on MATLAB's print command (using default
% options) in several ways:
%   - The figure borders are cropped

mathsam

function field2d = zonalMean2d(field3d)
% input is a filed3d(log,lat,p)
% return average field2d(lat,p)
numLog  = size(field3d,1);
field2d = reshape(sum(field3d,1)/numLog,size(field3d,2),size(field3d,3));
return

mathsam

function flux2dIsen = calIsenFlux(vcomp3d,PT3d,pfull,isenCoord,deltaP)
% calculate mass flux in isentropic coordinate
%----INPUT----
% vcomp3d:  meridonal velocity at one time
% PT3d:     potential temperature
% pfull:    temperature coordinate, unit mb
% isenCoord:isentropic coordinate that converts to, increases with index,
% should have regular spacing
% deltaP:   usually 10 mb

mathsam

% driver script for callsenFlux
temp  = ncread('Symmetric_HS=60K.nc','temp');
vcomp = ncread('Symmetric_HS=60K.nc','vcomp');
pfull = ncread('Symmetric_HS=60K.nc','pfull');

tStart    = 1000;
tEnd      = 1100;
isenCoord = 240:5:360;
deltaP    = 10;