API

radial_operators(nr, nℓ;
    r_in_frac = 0.6,
    r_out_frac = 0.985,
    ν = 5e11,
    trackingrate = :hanson2020,
    scalings = DefaultScalings,
    superadiabaticityparams = (;))

Compute the Chebyshev-basis representation of radial operators such as derivatives and factors that appear in the differential equations. The positional arguments nr and nℓ specify the number of radial and latitudinal coefficients in the stream functions.

Keyword arguments

• r_in_frac: inner boundary of the domain as a fraction of the solar radius
• r_out_frac: outer boundary of the domain as a fraction of the solar radius
• ν: coefficient of viscosity, in CGS units
• trackingrate: the tracking rate Ω0. This may either by a number, or one of :hanson2020 (corresponding to 453.1 nHz), :carrington (corresponding to 456 nHz), :cutoff (corresponding to the rotation rate at the equator at r = r_out_frac * Rsun), or :surface (corresponding to the rotation rate at the equator at r = Rsun).
• superadiabaticityparams: NamedTuple of parameters used to choose the superadiabaticity model. Values specified here are passed on to superadiabaticity, and these supersede the defaults in that function. This is empty by default.
• scalings: scalings applied to the equations to improve the balancing of the matrices to be diagonalized. Default values are (; Wscaling = 1e1, Sscaling = 1e3, Weqglobalscaling = 1e-4, Seqglobalscaling = 1.0), which should typically be good enough. The interpretation is that the fields that feature in the eigenvalue problem are (V, Wscaling * W, Sscaling * S) and the equations for W and S are multiplied by Weqglobalscaling and Seqglobalscaling respectively. Such scalings leave the eigenvalues unchanged.

RotMatrix(::Val{:matrix}, V_symmetric::Bool, rotation_profile::Symbol; operators, kw...)
RotMatrix(::Val{:spectrum}, V_symmetric::Bool, rotation_profile::Symbol; operators, kw...)

Struct containing the spectral decomposition of the rotation profile and the vorticity profile associated with differential rotation in the Sun, along with information about the parity of the solutions that are sought and the function that is to be used to generate the matrices.

A RotMatrix R may be called with arguments as R(args...; kwargs....). If R was created using Val(:matrix), this will return the matrix representation of the operator A that features in the eigenvalue problem $A\mathbf{v}=(ω/Ω_0)B\mathbf{v}$. These matrices are computed using the rotation profile that R is constructed with. On the other hand, if R was created using Val(:spectrum) as the first argument, this will return the filtered spectrum and eigenvectors.

mass_matrix(m::Integer;
    operators,
    V_symmetric::Bool)

Compute the mass matrix $B$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.

mass_matrix!(B, m::Integer;
    operators,
    V_symmetric::Bool)

Compute the mass matrix $B$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m. This operates in-place on the pre-allocated matrix B, and overwrites it with the result. The matrix B may be allocated using allocate_mass_matrix.

See mass_matrix for further details.

uniform_rotation_matrix(m::Integer; operators, V_symmetric::Bool)

Compute the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m, assuming that the Sun is rotating like a solid body. The matrix would therefore not contain any terms corresponding to differential rotation.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.

uniform_rotation_matrix!(A, m::Integer; operators, V_symmetric::Bool)

Compute the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m, assuming that the Sun is rotating like a solid body. This operates in-place on the pre-allocated matrix A, and overwrites it with the output. The matrix A may be allocated using allocate_operator_matrix.

See uniform_rotation_matrix for further details.

differential_rotation_matrix(m::Integer;
    operators, rotation_profile::Symbol, V_symmetric::Bool, kw...)

Return the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specified azimuthal order m, assuming a differentially rotating model of the Sun. The profile of rotation may be specified using the keyword argument rotation_profile.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• rotation_profile: label to select the rotation profile used to compute the spectrum. See solar_differential_rotation_profile_derivatives_grid for the list of possible options.

Additional keyword arguments will be passed on to the specific functions populating the matrix elements, such as solar_differential_rotation_terms!.

differential_rotation_matrix!(A, m::Integer;
    operators, rotation_profile::Symbol, V_symmetric::Bool, kw...)

Compute the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specified azimuthal order m assuming a differentially rotating model of the Sun, and store the result in A (which will be overwritten).

See differential_rotation_matrix for further details.

uniform_rotation_spectrum(m::Integer; operators, V_symmetric::Bool)

Compute the inertial-mode spectrum for the specified azimuthal order m, assuming that the Sun is rotating like a solid body, and it is being tracked from a frame that is rotating at the same rate as the Sun. This solves the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ and returns (ω_over_Ω0, v, (A,B)), where ω_over_Ω0 is a vector of eigenvalues, and v is a matrix whose columns are the eigenvectors. A is computed internally by calling uniform_rotation_matrix!, and B is computed by calling mass_matrix!.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• rotation_profile: label to select the rotation profile used to compute the spectrum. See solar_differential_rotation_profile_derivatives_grid for the list of possible options.

uniform_rotation_spectrum!((A, B), m::Integer; operators, V_symmetric::Bool)

Compute the inertial-mode spectrum for the specified azimuthal order m, assuming that the Sun is rotating like a solid body, and it is being tracked from a frame that is rotating at the same rate as the Sun. This overwrites the pre-allocated pencil matrices (A, B).

See uniform_rotation_spectrum for further details.

differential_rotation_spectrum(m::Integer; operators,
    V_symmetric::Bool, rotation_profile::Symbol, kw...)

Compute the inertial-mode spectrum for the specified azimuthal order m, assuming that the Sun is rotating with a solar-like rotation profile (the exact model may be specified using rotation_profile), and it is being tracked in a frame rotating at a rate Ω0 = operators.constants[:Ω0].

This solves the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ and returns (ω_over_Ω0, v, (A,B)), where ω_over_Ω0 is a vector of eigenvalues, and v is a matrix whose columns are the eigenvectors. A is computed internally by calling differential_rotation_matrix!, and B is computed by calling mass_matrix!.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• rotation_profile: label to select the rotation profile used to compute the spectrum. See solar_differential_rotation_profile_derivatives_grid for the list of possible options.

Additional keyword arguments are forwarded to differential_rotation_spectrum!.

differential_rotation_spectrum!((A, B), m::Integer; operators,
    V_symmetric::Bool, rotation_profile::Symbol, kw...)

Compute the inertial-mode spectrum for the specified azimuthal order m, assuming that the Sun is rotating with a solar-like rotation profile (the exact model may be specified using rotation_profile), and it is being tracked in a frame rotating at a rate Ω0 = operators.constants[:Ω0]. This overwrites the pre-allocated pencil matrices (A, B).

See differential_rotation_spectrum for further details.

filter_eigenvalues(ω_over_Ω0::AbstractVector{<:Number},
    v::AbstractMatrix{<:Number},
    (A,B), m::Integer;
    operators,
    V_symmetric::Bool,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators),
    filtercache = RossbyWaveSpectrum.allocate_filter_caches(m; operators, constraints),
    scale_eigenvectors::Bool = false,
    filterflags = RossbyWaveSpectrum.DefaultFilter,
    filterparams...)

Filter the set of eigenvalue-eigenvector pairs (ω_over_Ω0, v) for a specified azimuthal order m to remove the potentially spurious ones. This returns a filtered set (ω_over_Ω0f, vf). The eigenvalue pencil (A,B) should be provided if Filters.EIGEN in filterflags, and will be used in the filtering process. Otherwise, the parameters (A,B) are ignored.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• constraints: boundary condition constraints on the spectral coefficients, optional
• filtercache: pre-allocated workspace used in the filtering process, optional
• scale_eigenvectors: flag to indices whether to compute the unscaled stream functions from the scaled ones that are used in the eigenvalue problem.
• filterflags: the flags that specify which filters are used, optional. See Filters for a list of possible flags.
• filterparams: additional filter parameters, passed on to filterfn. See RossbyWaveSpectrum.DefaultFilterParams for the full list of parameters that may be specified, and their default values. The parameters specified in filterparams override the default values.

filter_eigenvalues(ω_over_Ω0s::AbstractVector{<:AbstractVector{<:Number}},
    vs::AbstractVector{<:AbstractMatrix{<:Number}},
    mr::AbstractVector{<:Integer};
    matrixfn!,
    operators,
    V_symmetric::Bool,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators),
    filtercache = RossbyWaveSpectrum.allocate_filter_caches(m; operators, constraints),
    filterflags = RossbyWaveSpectrum.DefaultFilter,
    filterparams...)

Filter the sets of eigenvalue-eigenvector pairs (ω_over_Ω0, v) for each azimuthal order m in the range mr, taken from the vectors ω_over_Ω0s and vs. The filtering operation for each m is independent of the others, and these are therefore carried out in parallel if multiple threads are available. This function returns a filtered set (ω_over_Ω0sf, vsf), where ω_over_Ω0sf[i] and vsf[i] represents the filtered set of eigenvalues and eigenvectors corresponding to m = mr[i].

Keyword arguments

• matrixfn!: The function used to compute the matrix A in the eigenvalue pencil (A,B). This may be one of uniform_rotation_matrix! or differential_rotation_matrix! or a RotMatrix, and will be called as matrixfn!(A, m; kw...) to populate the matrix A.
• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• constraints: Boundary condition constraints on the spectral coefficients, optional
• filtercache: Pre-allocated workspace used in the filtering process, optional
• filterflags: The flags that specify which filters are used, optional. See Filters for a list of possible flags.
• filterparams: additional filter parameters, passed on to filterfn. See RossbyWaveSpectrum.DefaultFilterParams for the full list of parameters that may be specified.

filter_eigenvalues(spectrumfn!,
    mr::AbstractVector;
    operators,
    V_symmetric,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators),
    filtercache = RossbyWaveSpectrum.allocate_filter_caches(m; operators, constraints),
    filterflags = RossbyWaveSpectrum.DefaultFilter,
    filterparams...)

Compute the spectrum of inertial waves for the specified range of azimuthal orders mr, and filter the solutions to remove potentially spurious solutions. This function returns (ω_over_Ω0s, vs), where ω_over_Ω0s[i] and vs[i] represent the filtered set of eigenvalues and eigenvectors corresponding to m = mr[i].

The argument spectrumfn! must be one of uniform_rotation_spectrum! or differential_rotation_spectrum!.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• constraints: Boundary condition constraints on the spectral coefficients, optional
• filtercache: Pre-allocated workspace used in the filtering process, optional
• filterflags: The flags that specify which filters are used, optional. See Filters for a list of possible flags.
• filterparams: additional filter parameters, passed on to filterfn. See RossbyWaveSpectrum.DefaultFilterParams for the full list of parameters that may be specified, and their default values. The parameters specified in filterparams override the default values.

save_eigenvalues(spectrumfn!, mr; operators, V_symmetric, rotation_profile, kw...)
save_eigenvalues(spectrumfn!::RotMatrix, mr; operators, kw...)

Compute the spectra for all azimuthal orders m in mr using the function spectrumfn!, and save the results to disk, with a filename given by rossbyeigenfilename Additional keyword arguments are passed on to filter_eigenvalues and rossbyeigenfilename.

If spectrumfn! isa RotMatrix, the keyword arguments V_symmetric and rotation_profile are automatically inferred from it.

FilteredEigen

Struct to hold the filtered eigenvalues and eigenfunctions for a range of m. Given a FilteredEigen F, the eigenfunctions for a specific m may be obtained using F[m], which returns a FilteredEigenSingleOrder.

FilteredEigen(filename::AbstractString)

Load the filtered eigenvalues and eigenvectors from filename.

datadir(filename)

Preprend the data directory to the filename. This is a convenience function to make loading files less cumbersome.

Filters

Flags that may be passed to the filterflags keyword argument in filter_eigenvalues to specify which conditions to filter solutions by.

Possible options are:

• Filters.NONE
• Filters.EIGEN
• Filters.EIGVAL
• Filters.BC
• Filters.SPATIAL_EQUATOR
• Filters.SPATIAL_RADIAL
• Filters.NODES

The individual filters may be combined using |, e.g. to include the Filters.EIGEN and the Filters.SPATIAL_EQUATOR flags, use Filters.EIGEN | Filters.SPATIAL_EQUATOR.

Filters.BC

Filter for eigenfunctions that satify the boundary conditions $C\mathbf{v} = 0$. Typically, all solutions should satisfy this condition, as this is imposed when solving the eigenvalue problem. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.DefaultFilter

The default set of filter flags that are used to constrain the set of solutions. This is defined as

julia> Filters.DefaultFilter
(NODES | SPATIAL_RADIAL | SPATIAL_EQUATOR | BC | EIGVEC | EIGVAL | EIGEN)::RossbyWaveSpectrum.Filters.FilterFlag = 0x00df

Filters.EIGEN

Filter for solutions (ω, v) that satisfy the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.EIGVAL

Filter for solutions that satisfy constraints on the eigenvalue, e.g. pick out damped solutions and ignore growing ones. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.EIGVAL

Filter for smooth eigenfunctions, by applying spectral cutoffs in the Chebyshev order n as well as the spherical-harmonic degree ℓ. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.NODES

Filter for eigenfunctions that have fewer than a specified number of radial nodes. A heuristic criteria is applied to ignore zero-crossings due to numerical or resolution limitations. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.NONE

Do not apply any filter. This retains all the solutions.

Filters.SPATIAL_EQUATOR

Filter for eigenfunctions that latitudinally peak about the solar equator. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

Filters.SPATIAL_RADIAL

Filter for eigenfunctions where the power is not concentrated very close to the top and the bottom boundaries of the domain. See RossbyWaveSpectrum.DefaultFilterParams for the parameters associated with this filter.

eigensystem_satisfy_filter(ω_over_Ω0::Number,
    v::StructVector{<:Complex},
    (A,B)::NTuple{2,AbstractMatrix{<:Number}},
    MVcache = RossbyWaveSpectrum.allocate_MVcache(size(A, 1));
    rtol = RossbyWaveSpectrum.DefaultFilterParams[:eigen_rtol])

Return whether the eigenvalue ω_over_Ω0 and eigenvector v satisfy $A\mathbf{v}=(\omega/\Omega_0) B\mathbf{v}$ to within the relative tolerance rtol.

eigvec_spectrum_filter(v::AbstractVector{<:Number}, m::Integer;
    operators,
    n_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:n_cutoff],
    Δl_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:Δl_cutoff],
    eigvec_spectrum_power_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:eigvec_spectrum_power_cutoff],
    filterfieldpowercutoff = RossbyWaveSpectrum.DefaultFilterParams[:filterfieldpowercutoff],
    low_n_power_lowercutoff = RossbyWaveSpectrum.DefaultFilterParams[:eigvec_spectrum_low_n_power_fraction_cutoff],
    kw...)

Return if the eigenvector v is smooth enough, according to the specified criteria.

eigvec_spectrum_filter(Feig::FilteredEigen, m::Integer, ind::Integer; kw...)

Return if the eigenvector ind-th eigenvector of Feig[m] is smooth enough, according to the specified criteria. By default, the filter parameters are taken to be identical to those used to compute the solutions in Feig, but each of these parameters may be overridden.

spatial_filter(v::AbstractVector{<:Number}, m::Integer;
    operators,
    filtercache = RossbyWaveSpectrum.allocate_filter_caches(m; operators),
    θ_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:θ_cutoff],
    equator_power_cutoff_frac = RossbyWaveSpectrum.DefaultFilterParams[:equator_power_cutoff_frac],
    pole_cutoff_angle = RossbyWaveSpectrum.DefaultFilterParams[:pole_cutoff_angle],
    pole_power_cutoff_frac = RossbyWaveSpectrum.DefaultFilterParams[:pole_power_cutoff_frac],
    filterfieldpowercutoff = RossbyWaveSpectrum.DefaultFilterParams[:filterfieldpowercutoff],
    radial_topbotpower_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:radial_topbotpower_cutoff],
    kw...)

Return if the eigenvector v for the azimuthal order m is spatially localized according to the specified parameters.

spatial_filter(Feig::FilteredEigen, m::Integer, ind::Integer; kw...)

Return if the ind-th eigenvector of Feig[m] is spatially localized according to the specified parameters. By default, the filter parameters are taken to be identical to those used to compute the solutions in Feig, but each of these parameters may be overridden.

eigenvalue_filter(ω_over_Ω0, m;
    eig_imag_unstable_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:eig_imag_unstable_cutoff],
    eig_imag_to_real_ratio_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:eig_imag_to_real_ratio_cutoff],
    eig_imag_stable_cutoff = RossbyWaveSpectrum.DefaultFilterParams[:eig_imag_stable_cutoff])

Return if the eigenvalue ω_over_Ω0 represents a damped eigenfunction.

Keyword arguments

• eig_imag_unstable_cutoff: lower cutoff below which modes are considered growing. This is slightly less than zero to account for numerical errors.
• eig_imag_to_real_ratio_cutoff: cutoff ratio of linewidth to central frequency to filter out strongly damped modes.
• eig_imag_stable_cutoff: absolute cutoff on linewidth to filter out strongly damped modes.

boundary_condition_filter(v::StructVector{<:Complex},
    BC::AbstractMatrix{<:Real},
    BCVcache::StructVector{<:Complex} = RossbyWaveSpectrum.allocate_BCcache(size(BC,1)),
    atol = 1e-5)

Return if the eigenfunction with spectral coefficients v satisfies the boundary conditions within the absoute tolerance atol. The boundary-condition matrix BC may be obtained from constraintmatrix.

constant_differential_rotation_terms!(A, m::Integer;
    operators,
    V_symmetric::Bool,
    ΔΩ_frac = 0.01)

Compute the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m, assuming that the Sun is rotating like a solid body, albeit at a rotation speed $Ω$ that exceeds the tracking rate $Ω_0$ by the specified ratio ΔΩ_frac (i.e. Ω = Ω0 * (1 + ΔΩ_frac)). This function serves as a validation test for the differentially rotating model.

This function overwrites the input matrix A that is obtained as the output to uniform_rotation_matrix, and adds the extra terms corresponding to the difference in rotation rates to it.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• ΔΩ_frac: Ratio by which the background rotation rate $Ω$ exceeds the reference rate $Ω_0$. The relation between the two is given by Ω = Ω_0 * (1 + ΔΩ_frac).

DefaultFilterParams

Default list of parameters that are used in the filtering process. The full list is:

const DefaultFilterParams = Dict(
    # boundary condition filter: Filters.BC
    :bc_atol => 1e-5, # absolute tolerance to which boundary conditions must be satisfied
    # eigval filter: Filters.EIGVAL
    :eig_imag_unstable_cutoff => -1e-6, # absolute lower cutoff below which modes are considered to be exponentially growing
    :eig_imag_to_real_ratio_cutoff => 3, # relative upper cutoff on linewidths
    :eig_imag_stable_cutoff => 0.2, # absolute upper cutoff on linewidths
    # eigensystem satisfy filter: Filters.EIGEN
    :eigen_rtol => 0.01, # tolerance above which modes are rejected
    # smooth eigenvector filter: Filters.EIGVEC
    :Δl_cutoff => 7, # spherical harmonic degree
    :n_cutoff => 10, # chebyshev order
    :eigvec_spectrum_power_cutoff => 0.5, # fraction of total power required below cutoffs
    :eigvec_spectrum_low_n_power_fraction_cutoff => 1, # maximum ratio of total power above n_cutoff to total power below n_cutoff
    # spatial localization filter: Filters.SPATIAL
    :θ_cutoff => deg2rad(45), # colatitude below which a mode is considered equatorial
    :equator_power_cutoff_frac => 0.4, # minimum fraction of power that lie below θ_cutoff for equatorial modes
    :pole_cutoff_angle => deg2rad(25), # colatitude above which a mode is considered polar
    :pole_power_cutoff_frac => 0.05, # maximum fraction of power that may be near poles for equatorial modes
    :radial_topbotpower_cutoff => 0.7, # upper cutoff for power in the top and bottom 10% of the radial domain
    # radial nodes filter # Filters.NODES
    :nnodesmax => 5, # maximum number of radial nodes in either the real or imaginary component of the eigenfunction
    :nodessmallpowercutoff => 0.05, # minimum fractional power between roots beyond which a zero-crossing is considered real
    # exclude a field from filters if relative power is below a cutoff
    :filterfieldpowercutoff => 1e-2,
)

Each parameter may be passed as a keyword argument to filter_eigenvalues.

filterfn(λ::Number,
    v::AbstractVector{<:Number},
    m::Integer,
    (A,B);
    operators,
    V_symmetric,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators),
    filtercache = RossbyWaveSpectrum.allocate_filter_caches(m; operators, constraints),
    filterflags = RossbyWaveSpectrum.DefaultFilter,
    filterparams...)

Return whether one eigenvalue and eigenvector pair (λ, v) for the azimuthal order m satisfy the chosen filters. The matrices (A,B) represent the pencil.

Keyword arguments

• operators: obtained as the output of radial_operators.
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• constraints: boundary condition constraints on the spectral coefficients, optional.
• filtercache: pre-allocated workspace used in the filtering process, optional.
• filterflags: the flags that specify which filters are used, optional. See Filters for a list of possible flags.
• filterparams: additional filter parameters to override the default ones from RossbyWaveSpectrum.DefaultFilterParams. The parameters specified in filterparams will override the default values, and the ones left unspecified will be set to the default values.

solar_differential_rotation_terms!(A, m::Integer;
    operators,
    V_symmetric::Bool,
    rotation_profile::Symbol,
    ΔΩ_frac = 0.01, # only used to test the constant case
    ΔΩ_scale = 1, # scale the diff rot profile
    kw...)

Compute the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ for a specific m, assuming that the Sun is rotating with a solar-like rotation profile. The keyword argument rotation_profile specifies the specific model that is chosen. Optionally, additional keyword arguments kw may be passed to alter the specifics of the model. These are detailed below.

This function overwrites the input matrix A that is obtained as the output to uniform_rotation_matrix, and adds the extra terms corresponding to differential rotation to it.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• rotation_profile: label to select the rotation profile used to compute the spectrum. See solar_differential_rotation_profile_derivatives_grid for the list of possible options.
• ΔΩ_frac: ratio by which the background rotation rate $Ω$ exceeds the reference rate $Ω_0$, optional. The relation between the two is given by Ω = Ω_0 * (1 + ΔΩ_frac). This is only used if rotation_profile == :constant, and is chosen to be 0.01 by default.
• ΔΩ_scale: factor by which the differential rotation profile may be scaled. This may be used to tune the differential rotation to trace the results back to the solid-body rotation case.

solar_differential_rotation_profile_derivatives_grid(;
    operators,
    rotation_profile,
    ΔΩ_frac = 0.01,
    kw...)

Return the smoothed spatial profile of rotation, corresponding to the model specified by rotation_profile. This function reads in rotation model files, and interpolates the profile on to Chebyshev nodes in both the radial and angular coordinates through smoothing splines. This spatial profile is subsequently transformed to spectral coefficients in solar_differential_rotation_profile_derivatives_Fun.

Typically, this is used to load the solar rotation profile, but other profiles may be specified for testing the code.

Keyword arguments

• operators: obtained as the output of radial_operators
• rotation_profile: the flag that chooses the model of the rotation profile. Possible options are:
  • :solar_latrad: Smoothed solar rotation profile, limited to the radial domain
  • :solar_latrad_squished: Same as latrad, except the solar surface is projected to the outer boundary of the radial domain.
  • :solar_radial_equator: smoothed radial profile of the solar equatorial rotation rate, but limited to the radial domain. This corresponds to the radial profile of the latrad model at the equator.
  • :solar_radial_equator_squished: same profile as radial_equator, except that the solar surface is projected onto the outer boundary of the radial domain. This corresponds to the radial profile of the latrad_squished model at the equator.
  • :solar_constant: background medium rotates like a solid body, at a rate Ω that differs from the tracking rate Ω0
• ΔΩ_frac: In case rotation_profile == :solar_constant, the factor by which the tracking rate Ω0 is increased to obtain the rotation rate of the background medium.

solar_differential_rotation_profile_derivatives_Fun(;
    operators,
    rotation_profile,
    ΔΩ_smoothing_param = 5e-2,
    kw...)

Compute the Chebyshev-Ultraspherical decomposition of the rotation profile specified by rotation_profile, which typically is a smoothed version of the solar rotation profile. The keyword arguments kw are passed on to solar_differential_rotation_profile_derivatives_grid.

Returns a collection (; ΔΩ, dr_ΔΩ, d2r_ΔΩ, dz_ΔΩ), where the first term is the spectral approximation to the profile of differential rotation, and the subsequent terms are radial and z-derivatives.

Keyword arguments

• operators: obtained as the output of radial_operators
• rotation_profile: the flag that chooses the model of the rotation profile. See solar_differential_rotation_profile_derivatives_grid for possible options.
• ΔΩ_smoothing_param: the extent of smoothing applied to the rotation profile before performing the spectral transform. Typically, this should be small to match the original profile closely, with the tradeoff being in the number of spectral coefficients that are necessary to represent the model.
• kw: Optional keyword arguments that are passed on to solar_differential_rotation_profile_derivatives_grid.

solar_differential_rotation_profile_derivatives_Fun(Feig::FilteredEigen; kw...)

Return the profile of background rotation that was used to compute the spectrum contained in Feig. Additional keyword arguments kw may be supplied to override the ones in Feig.

solar_differential_rotation_vorticity_Fun(; operators, ΔΩprofile_deriv)

Return the spectral approximation to the vorticity associated with differential rotation $\left(\boldsymbol{ω}_\Omega = ∇ × (Ω \hat{z})\right)$, given the rotation profile $Ω$ and its derivatives.

Keyword arguments

• operators: obtained as the output of radial_operators
• ΔΩprofile_deriv: Collection containing the rotation profile and its derivatives, obtained from solar_differential_rotation_profile_derivatives_Fun.

solar_differential_rotation_vorticity_Fun(Feig::FilteredEigen; kw...)

Return the profile of vorticity associated with background rotation that was used to compute the spectrum contained in Feig. Additional keyword arguments in kw may be used to override the ones in Feig.

constraintmatrix(operators)

Return a collection containing the spectral representation of the boundary condition operator BC, and a basis ZC whose columns lie in its null-space. The operators may be obtained using radial_operators.

The constraint satisfied by the eigenfunctions may be expressed as $C\mathbf{v}=0$. Given a matrix $Z$ that satisfies $CZ=0$, we may express the eigenfunction as $\mathbf{v}=Z\mathbf{w}$ for an arbitrary $\mathbf{w}$. This function uses the variable names BC for $C$ and ZC for $Z$.

The matrices may be destructured as (; BC, ZC) = RossbyWaveSpectrum.constraintmatrix(operators).

superadiabaticity(r::Real;
    r_out = Rsun,
    δcz = 3e-6,
    δtop = 3e-5,
    δrad = -1e-3,
    dtrans = 0.05Rsun,
    dtop = 0.05Rsun,
    r_sub = 0.8Rsun,
    r_tran = 0.725Rsun)

Compute the superadiabaticity profile $δ$ parameterized by the keyword arguments, following Rempel (2005), ApJ 622:1320 –1332, section 2.3.

Keyword arguments

• δcz: representative value of superadiabaticity $δ$ in the bulk of the convection zone
• δtop: representative value of superadiabaticity $δ$ at the surface
• δrad: representative value of superadiabaticity $δ$ in the radiative interior
• dtrans: scale over which the profile in the radiative zone transitions to that in the convection zone
• dtop: scale over which the profile in the convection zone transitions to that at the surface
• r_tran: radius below which the radiative profile transitions to strongly subadiabatic (in the overshoot region)
• r_sub: radius below which the radiative profile transitions to weakly subadiabatic within the convection zone

eigenfunction_realspace(v, m; operators, V_symmetric::Bool,
    n_lowpass_cutoff = nothing,
    Δl_lowpass_cutoff = nothing)

Compute the real-space eigenfunctions from the spectral coefficients v for the azimuthal order m.

Keyword arguments

• operators: obtained as the output of radial_operators
• V_symmetric: set to true or false depending on whether the stream function V is latitudinally symmetric or antisymmtric about the equator.
• n_lowpass_cutoff: maximum chebyshev order n to retain in the inverse transform, optional. If no value is provided, all the coefficients are used (default).
• Δl_lowpass_cutoff: maximum harmonic order ℓ to retain in the inverse transform, optional. If no value is provided, all the coefficients are used (default).

eigenfunction_realspace(Feig::FilteredEigen, m, ind; kw...)

Return the real-space eigenfunction from the ind-th eigenvector for the azimuthal order m, taken from the full list of solutions Feig.

compute_constrained_matrix(A, constraints)

Compute the constrained matrix $Z^T A Z$ given A. Here the columns of Z = constraints.ZC represents the spectral coefficients of a basis that satifies the boundary conditions.

compute_constrained_matrix!(out, constraints, A)

Compute the constrained matrix $Z^T A Z$ given A, and overwrite out with the result. Here the columns of Z = constraints.ZC represents the spectral coefficients of a basis that satifies the boundary conditions.

constrained_eigensystem((A, B);
    operators,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators),
    cache = RossbyWaveSpectrum.constrained_matmul_cache(constraints),
    temp_projectback = RossbyWaveSpectrum.allocate_projectback_temp_matrices(size(constraints.ZC)),
    kw...)

Solve the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$ subject to the boundary-condition constraint $C\mathbf{v}=0$, and return (λ, v, (A, B)). The input arguments are the matrix pencil (A, B).

Internally, this computes a matrix Z = constraints.ZC that represents the spectral coefficients of a basis that satifies the boundary conditions (i.e. this satisfies $CZ=0$). The eigenfunction $\mathbf{v}$ may therefore be expressed as $\mathbf{v} = Z\mathbf{w}$ for an arbitrary $\mathbf{w}$. Substituting this, and projecting the eigenvalue problem on to the basis $Z$, we obtain $(Z^T A Z)\mathbf{w} = (\omega/\Omega_0)(Z^T B Z)\mathbf{w}$, which is solved numerically. Finally, the eigenfunctions are transformed to the Chebyshev basis using $\mathbf{v} = Z\mathbf{w}$.

Keyword arguments

• operators: obtained as the output of radial_operators
• constraints: projection matrices to enforce the boundary conditions, optional.
• cache: temporary arrays used in evaluating the constrained matrix pencil (Z^T A Z, Z^T B Z), optional. See constrained_matmul_cache
• temp_projectback: temporary arrays used in computing v from w, optional.

constrained_matmul_cache(constraints)

Allocate temporary matrices that are overwritten with $Z^T A Z$ and $Z^T B Z$, given the eigenvalue pencil $(A, B)$, where the columns of Z = constraints.ZC represent the spectral coefficients of a basis that satifies the boundary conditions.

allocate_operator_matrix(operators)

Allocate the operator matrix $A$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$. The argument operators is obtained as the output of radial_operators.

allocate_mass_matrix(operators)

Allocate the mass matrix $B$ that features in the eigenvalue problem $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$. The argument operators is obtained as the output of radial_operators.

allocate_operator_mass_matrices(operators)

Allocate the eigenvalue pencil matrices $(A,B)$. The argument operators is obtained as the output of radial_operators.

allocate_filter_caches(m; operators,
    constraints = RossbyWaveSpectrum.constraintmatrix(operators))

Allocate temporary matrices used in the filtering process, given the azimuthal order m.

Keyword arguments

• operators: obtained as the output of radial_operators
• constraints: boundary condition constraints on the spectral coefficients, optional

allocate_projectback_temp_matrices(sz)

Allocate temporary matrices used to project the eigenfunctions for $(Z^T A Z)\mathbf{w} = (\omega/\Omega_0)(Z^T B Z)\mathbf{w}$ to those for $A\mathbf{v}=(\omega/\Omega_0)B\mathbf{v}$.

operator_matrices(Feig::FilteredEigen, m; kw...)

Return the pencil (A, B) for the azimuthal order m, using the parameters used to compute the eigenvalues in Feig.

rossbyeigenfilename(; operators, V_symmetric, rotation_profile, modeltag = "")

Return the file name to which the results should be written. The output is a combination of the keyword arguments specifed.

FilteredEigenSingleOrder

Struct containing the filtered set of eigenvalues and eigenfunctions for a single m. A FilteredEigenSingleOrder is equivalent to a Vector of serial => (eigenvalue, eigenfunction) pairs, with additional metadata stored along with this. Indexing into a FilteredEigenSingleOrder with a serial number should return a (eigenvalue, eigenfunction) pair.

colatitude_grid(m::Integer, operators)

Return the colatitude $\left(\theta\right)$ grid that is used to compute eigenfunctions in real space for the specific azimuthal order m.

allocate_MVcache(nrows)

Allocate temporary matrices that are used to compute matrix products in-place. These are used in filtering the solutions (See Filters.EIGEN).

allocate_BCcache(nrows)

Allocate temporary matrices that are used to compute matrix products in-place. These are used in filtering the solutions (See Filters.BC).