Spectral tansforms

class magic.spectralTransforms.SpectralTransforms(l_max=32, minc=1, lm_max=561, n_theta_max=64, m_max=None, verbose=True)[source]

This python class is used to compute Legendre and Fourier transforms from spectral to physical space. It works in two steps: one first needs to initialize the transform

>>> sh = SpectralTransforms( l_max=256, lm_max=33153, n_theta_max=384)
>>> print(Tlm[:, 10].shape) # lm_max (Temperature at ir=10)
>>> T = sh.spec_spat(Tlm) # T[n_phi_max, n_theta_max]
spat_spec(*args)[source]

This subroutine computes a transfrom from spatial representation (n_phi,n_theta) to spectral representation (lm_max). It returns one complex 1-D array (dimension(n_phi_max))

>>> gr = MagicGraph()
>>> sh = SpectralTransforms(gr.l_max, gr.minc, gr.lm_max, gr.n_theta_max)
>>> vr = gr.vr[:,:,30] # Radius ir=30
>>> vrlm = sh.spat_spec(vr) # vrlm is a complex array (lm_max)
>>> # Caculation of the poloidal potential from vr:
>>> wlm = np.zeros_like(vrlm)
>>> wlm[1:] = vrlm[1:]/(sh.ell[1:]*(sh.ell[1:]+1))*gr.radius[30]**2
>>> # Spheroidal/Toroidal transform
>>> vtlm, vplm = spec_spat(gr.vtheta, gr.vphi)
Parameters:

input (numpy.ndarray) – input array in the physical space (n_phi,n_theta)

Returns:

output array in the spectral space (lm_max)

Return type:

numpy.ndarray

spec_spat(*args, **kwargs)[source]

This subroutine computes a transfrom from spectral to spatial for all latitudes. It returns either one or two 2-D arrays (dimension(n_phi_max,n_theta_max)) depending if only the poloidal or both the poloidal and the toroidal potentials are given as input quantities.

>>> print(wlmr.shape) # lm_max
>>> vr = spec_spat_equat(wlmr)
>>> print(vr.shape) # n_phi, n_theta
>>> vt, vp = spec_spat_equat(dwdrlmr, zlmr)
spec_spat_dphi(polo)[source]

This routine computes the phi-derivative and the transform from spectral to spatrial spaces. It returns a 2-D array of dimension (n_phi,n_theta)

>>> p = MagicPotential('V')
>>> vrlm = p.pol*p.ell*(p.ell+1)/p.radius[ir]**2/p.rho0[ir] # vr at r=ir
>>> dvrdp = p.sh.spec_spat_dphi(vrlm) # phi-derivative of vr
Parameters:

polo (numpy.ndarray) – the input array(lm_max) in spectral space

Returns:

the phi derivative in the physical space (n_phi, n_theta)

Return type:

numpy.ndarray

spec_spat_dtheta(polo, l_axi=False)[source]

This routine computes the theta-derivative and the transform from spectral to spatrial spaces. It returns a 2-D array of dimension (n_phi,n_theta)

>>> p = MagicPotential('V')
>>> vrlm = p.pol*p.ell*(p.ell+1)/p.radius[ir]**2/p.rho0[ir] # vr at r=ir
>>> dvrdt = p.sh.spec_spat_dtheta(vrlm) # theta-derivative of vr
Parameters:
  • polo (numpy.ndarray) – the input array(lm_max) in spectral space

  • l_axi (bool) – switch to True, if only the axisymmetric field is needed

Returns:

the theta derivative in the physical space (n_phi, n_theta)

Return type:

numpy.ndarray

spec_spat_equat(*args)[source]

This subroutine computes a transfrom from spectral to spatial at the equator. It returns either one or two 1-D arrays (dimension(n_phi_max)) depending if only the poloidal or both the poloidal and the toroidal potentials are given as input quantities.

>>> print(wlmr.shape) # lm_max
>>> vr = spec_spat_equat(wlmr)
>>> print(vr.shape) # n_phi
>>> vt, vp = spec_spat_equat(dwdrlmr, zlmr)

Plotting functions

magic.plotlib.cut(dat, vmax=None, vmin=None)[source]

This functions truncates the values of an input array that are beyond vmax or below vmin and replace them by vmax and vmin, respectively.

>>> # Keep only values between -1e3 and 1e3
>>> datNew = cut(dat, vmin=-1e3, vmax=1e3)
Parameters:
  • dat (numpy.ndarray) – an input array

  • vmax (float) – maximum upper bound

  • vmin (float) – minimum lower bound

Returns:

an array where the values >=vmax have been replaced by vmax and the values <=vmin have been replaced by vmin

Return type:

numpy.ndarray

magic.plotlib.default_cmap(field)[source]

This function selects a default colormap for an input field. This allows to have different colormaps depending on the quantity one wants to plot. This can make use of cmocean colormaps, whenever installed.

Parameters:

field (str) – the name of input field

Returns cm:

the name of the matplotlib colormap

Rtype cm:

str

magic.plotlib.diverging_cmap(field)[source]

This function determines whether the data are sequential or diverging (i.e. centered around zero). In the latter, colormaps will be by default centered.

Parameters:

field (str) – the name of input field

Returns diverging:

a boolean to say whether the colormap will be centered or not

Rtype cm:

bool

magic.plotlib.equatContour(data, radius, minc=1, label=None, levels=65, cm='seismic', normed=True, vmax=None, vmin=None, cbar=True, tit=True, normRad=False, deminc=True, bounds=True)[source]

Plot the equatorial cut of a given field

Parameters:
  • data (numpy.ndarray) – the input data (an array of size (nphi,nr))

  • radius (numpy.ndarray) – the input radius

  • minc (int) – azimuthal symmetry

  • label (str) – the name of the input physical quantity you want to display

  • normRad (bool) – when set to True, the contour levels are normalised radius by radius (default is False)

  • levels (int) – the number of levels in the contour

  • cm (str) – name of the colormap (‘jet’, ‘seismic’, ‘RdYlBu_r’, etc.)

  • tit (bool) – display the title of the figure when set to True

  • cbar (bool) – display the colorbar when set to True

  • vmax (float) – maximum value of the contour levels

  • vmin (float) – minimum value of the contour levels

  • normed (bool) – when set to True, the colormap is centered around zero. Default is True, except for entropy/temperature plots.

  • deminc (bool) – a logical to indicate if one wants do get rid of the possible azimuthal symmetry

  • bounds (bool) – a boolean to determine if one wants to plot the limits of the domain (True by default)

magic.plotlib.hammer2cart(ttheta, pphi, colat=False)[source]

This function is used to define the Hammer projection used when plotting surface contours in magic.Surf

>>> # Load Graphic file
>>> gr = MagicGraph()
>>> # Meshgrid
>>> pphi, ttheta = mgrid[-np.pi:np.pi:gr.nphi*1j, np.pi/2.:-np.pi/2.:gr.ntheta*1j]
>>> x,y = hammer2cart(ttheta, pphi)
>>> # Contour plots
>>> contourf(x, y, gr.vphi)
Parameters:
  • ttheta (numpy.ndarray) – meshgrid [nphi, ntheta] for the latitudinal direction

  • pphi – meshgrid [nphi, ntheta] for the azimuthal direction

  • colat (numpy.ndarray) – colatitudes (when not specified a regular grid is assumed)

Returns:

a tuple that contains two [nphi, ntheta] arrays: the x, y meshgrid used in contour plots

Return type:

(numpy.ndarray, numpy.ndarray)

magic.plotlib.merContour(data, radius, label=None, levels=65, cm='seismic', normed=True, vmax=None, vmin=None, cbar=True, tit=True, fig=None, ax=None, bounds=True, lines=False)[source]

Plot a meridional cut of a given field

Parameters:
  • data (numpy.ndarray) – the input data (an array of size (ntheta,nr))

  • radius (numpy.ndarray) – the input radius

  • label (str) – the name of the input physical quantity you want to display

  • levels (int) – the number of levels in the contour

  • cm (str) – name of the colormap (‘jet’, ‘seismic’, ‘RdYlBu_r’, etc.)

  • tit (bool) – display the title of the figure when set to True

  • cbar (bool) – display the colorbar when set to True

  • vmax (float) – maximum value of the contour levels

  • vmin (float) – minimum value of the contour levels

  • normed (bool) – when set to True, the colormap is centered around zero. Default is True, except for entropy/temperature plots.

  • bounds (bool) – a boolean to determine if one wants to plot the limits of the domain (True by default)

  • fig (matplotlib.figure.Figure) – a pre-existing figure (if needed)

  • ax (matplotlib.axes._subplots.AxesSubplot) – a pre-existing axis

  • lines (bool) – when set to True, over-plot solid lines to highlight the limits between two adjacent contour levels

magic.plotlib.radialContour(data, rad=0.85, label=None, proj='hammer', lon_0=0.0, vmax=None, vmin=None, lat_0=30.0, levels=65, cm='seismic', normed=True, cbar=True, tit=True, lines=False, fig=None, ax=None)[source]

Plot the radial cut of a given field

Parameters:
  • data (numpy.ndarray) – the input data (an array of size (nphi,ntheta))

  • rad (float) – the value of the selected radius

  • label (str) – the name of the input physical quantity you want to display

  • proj (str) – the type of projection. Default is Hammer, in case you want to use ‘ortho’ or ‘moll’, then Basemap is required.

  • levels (int) – the number of levels in the contour

  • cm (str) – name of the colormap (‘jet’, ‘seismic’, ‘RdYlBu_r’, etc.)

  • tit (bool) – display the title of the figure when set to True

  • cbar (bool) – display the colorbar when set to True

  • lines (bool) – when set to True, over-plot solid lines to highlight the limits between two adjacent contour levels

  • vmax (float) – maximum value of the contour levels

  • vmin (float) – minimum value of the contour levels

  • normed (bool) – when set to True, the colormap is centered around zero. Default is True, except for entropy/temperature plots.

  • fig (matplotlib.figure.Figure) – a pre-existing figure (if needed)

  • ax (matplotlib.axes._subplots.AxesSubplot) – a pre-existing axis

Various useful functions

magic.libmagic.ReadBinaryTimeseries(infile, ncols, datatype='f8', endianness='>')[source]

This function reads binary timeseries. It is then faster than the fast_read function.

Parameters:
  • infile (string) – the file to read

  • ncols (int) – number of columns of the file

  • datatype (string) – ‘f8’ = 64-bit floating-point number ‘f4’ = 32-bit floating-point number

  • endianness (string) – ‘>’ = big-endian ; ‘<’ = small-endian

Returns:

an array[nlines, ncols] that contains the data of the binary file

Return type:

numpy.ndarray

magic.libmagic.anelprof(radius, strat, polind, g0=0.0, g1=0.0, g2=1.0)[source]

This functions calculates the reference temperature and density profiles of an anelastic model.

>>> rad = chebgrid(65, 1.5, 2.5)
>>> temp, rho, beta = anelprof(rad, strat=5., polind=2.)
Parameters:
  • radius (numpy.ndarray) – the radial gridpoints

  • polind (float) – the polytropic index

  • strat (float) – the number of the density scale heights between the inner and the outer boundary

  • g0 (float) – gravity profile: g=g0

  • g1 (float) – gravity profile: g=g1*r/r_o

  • g2 (float) – gravity profile: g=g2*(r_o/r)**2

Returns:

a tuple that contains the temperature profile, the density profile and the log-derivative of the density profile versus radius

Return type:

(numpy.ndarray, numpy.ndarray, numpy.ndarray)

magic.libmagic.avgField(time, field, tstart=None, std=False, fix_missing_series=False)[source]

This subroutine computes the time-average (and the std) of a time series

>>> ts = MagicTs(field='misc', iplot=False, all=True)
>>> nuavg = avgField(ts.time, ts.topnuss, 0.35)
>>> print(nuavg)
Parameters:
  • time (numpy.ndarray) – time

  • field (numpy.ndarray) – the time series of a given field

  • tstart (float) – the starting time of the averaging

  • std (bool) – when set to True, the standard deviation is also calculated

  • fix_missing_series (bool) – when set to True, data equal to zero are ignored, this is done in case new columns have been added to the time series

Returns:

the time-averaged quantity

Return type:

float

magic.libmagic.chebgrid(nr, a, b)[source]

This function defines a Gauss-Lobatto grid from a to b.

>>> r_icb = 0.5 ; r_cmb = 1.5; n_r_max=65
>>> rr = chebgrid(n_r_max, r_icb, r_cmb)
Parameters:
  • nr (int) – number of radial grid points plus one (Nr+1)

  • a (float) – lower limit of the Gauss-Lobatto grid

  • b (float) – upper limit of the Gauss-Lobatto grid

Returns:

the Gauss-Lobatto grid

Return type:

numpy.ndarray

magic.libmagic.cylSder(radius, data, order=4)[source]

This function computes the s derivative of an input array defined on a regularly-spaced cylindrical grid.

>>> s = linspace(0., 1., 129 ; dat = cos(s)
>>> ddatds = cylSder(s, dat)
Parameters:
  • radius (numpy.ndarray) – cylindrical radius

  • data (numpy.ndarray) – input data

  • order (int) – order of the finite-difference scheme (possible values are 2 or 4)

Returns:

s derivative

Return type:

numpy.ndarray

magic.libmagic.cylZder(z, data)[source]

This function computes the z derivative of an input array defined on a regularly-spaced cylindrical grid.

>>> z = linspace(-1., 1., 129 ; dat = cos(z)
>>> ddatdz = cylZder(z, dat)
Parameters:
  • z (numpy.ndarray) – height of the cylinder

  • data (numpy.ndarray) – input data

Returns:

z derivative

Return type:

numpy.ndarray

magic.libmagic.fast_read(file, skiplines=0, binary=False, precision=<class 'numpy.float64'>)[source]

This function reads an input ascii table (can read both formatted or unformatted fortran)

>>> # Read 'e_kin.test', skip the first 10 lines
>>> data = fast_read('e_kin.test', skiplines=10)
Parameters:
  • file (str) – name of the input file

  • skiplines (int) – number of header lines to be skept during reading

  • binary (bool) – when set to True, try to read an unformatted binray Fortran file (default is False)

  • precision (str) – single (np.float32) or double precision (np.float64)

Returns:

an array[nlines, ncols] that contains the data of the ascii file

Return type:

numpy.ndarray

magic.libmagic.fd_grid(nr, a, b, fd_stretching=0.3, fd_ratio=0.1)[source]

This function defines a stretched grid between a and b

>>> r_icb = 0.5 ; r_cmb = 1.5; n_r_max=64
>>> rr = fd_grid(n_r_max, r_cmb, r_icb)
Parameters:
  • nr (int) – number of radial grid points

  • a (float) – upper boundary of the grid

  • b (float) – lower boundary of the grid

  • fd_stretching (float) – fraction of points in the bulk

  • fd_ratio (float) – ratio of minimum to maximum spacing

Returns:

the radial grid

Returns:

the radial grid

Return type:

numpy.ndarray

magic.libmagic.getCpuTime(file)[source]

This function calculates the CPU time from one given log file

Parameters:

file (file) – the log file you want to analyze

Returns:

the total CPU time

Return type:

float

magic.libmagic.getTotalRunTime()[source]

This function calculates the total CPU time of one run directory

Returns:

the total RUN time

Return type:

float

magic.libmagic.intcheb(f, nr, z1, z2)[source]

This function integrates an input function f defined on the Gauss-Lobatto grid.

>>> print(intcheb(f, 65, 0.5, 1.5))
Parameters:
  • f – an input array

  • nr (int) – number of radial grid points

  • z1 (float) – lower limit of the Gauss-Lobatto grid

  • z2 (float) – upper limit of the Gauss-Lobatto grid

Type:

numpy.ndarray

Returns:

the integrated quantity

Return type:

float

magic.libmagic.matder(nr, z1, z2)[source]

This function calculates the derivative in Chebyshev space.

>>> r_icb = 0.5 ; r_cmb = 1.5; n_r_max=65
>>> d1 = matder(n_r_max, r_icb, r_cmb)
>>> # Chebyshev grid and data
>>> rr = chebgrid(n_r_max, r_icb, r_cmb)
>>> f = sin(rr)
>>> # Radial derivative
>>> df = dot(d1, f)
Parameters:
  • nr (int) – number of radial grid points

  • z1 (float) – lower limit of the Gauss-Lobatto grid

  • z2 (float) – upper limit of the Gauss-Lobatto grid

Returns:

a matrix of dimension (nr,nr) to calculate the derivatives

Return type:

numpy.ndarray

magic.libmagic.phideravg(data, minc=1, order=4)[source]

phi-derivative of an input array

>>> gr = MagicGraph()
>>> dvphidp = phideravg(gr.vphi, minc=gr.minc)
Parameters:
  • data (numpy.ndarray) – input array

  • minc (int) – azimuthal symmetry

  • order (int) – order of the finite-difference scheme (possible values are 2 or 4)

Returns:

the phi-derivative of the input array

Return type:

numpy.ndarray

magic.libmagic.prime_factors(n)[source]

This function returns all prime factors of a number

Returns:

all prime factors

Return type:

list

magic.libmagic.progressbar(it, prefix='', size=60)[source]

Fancy progress-bar for loops

for i in progressbar(range(1000000)):
    x = i
Parameters:
  • prefix (str) – prefix string before progress bar

  • size (int) – width of the progress bar (in points of xterm width)

magic.libmagic.rderavg(data, rad, exclude=False)[source]

Radial derivative of an input array

>>> gr = MagiGraph()
>>> dvrdr = rderavg(gr.vr, gr.radius)
Parameters:
  • data (numpy.ndarray) – input array

  • rad (numpy.ndarray) – radial grid

  • exclude (bool) – when set to True, exclude the first and last radial grid points and replace them by a spline extrapolation (default is False)

Returns:

the radial derivative of the input array

Return type:

numpy.ndarray

magic.libmagic.scanDir(pattern, tfix=None)[source]

This function sorts the files which match a given input pattern from the oldest to the most recent one (in the current working directory)

>>> dat = scanDir('log.*')
>>> print(log)
Parameters:
  • pattern (str) – a classical regexp pattern

  • tfix (float) – in case you want to add only the files that are more recent than a certain date, use tfix (computer 1970 format!!)

Returns:

a list of files that match the input pattern

Return type:

list

magic.libmagic.sderavg(data, rad, colat=None, exclude=False)[source]

s derivative of an input array

>>> gr = MagiGraph()
>>> dvpds = sderavg(gr.vphi, eta=gr.radratio, colat=gr.colatitude)
Parameters:
  • data (numpy.ndarray) – input array

  • rad (numpy.ndarray) – radial grid

  • exclude (bool) – when set to True, exclude the first and last radial grid points and replace them by a spline extrapolation (default is False)

  • colat (numpy.ndarray) – colatitudes (when not specified a regular grid is assumed)

Returns:

the s derivative of the input array

Return type:

numpy.ndarray

magic.libmagic.secondtimeder(time, y)[source]

second time derivative of an input array

computed with central differences (numpy.gradient)

>>> ts = MagicTs(field='e_kin')
>>> dt_ekinpol = secondtimeder(ts,field='ekin_pol')
magic.libmagic.selectField(obj, field, labTex=True, ic=False)[source]

This function selects for you which field you want to display. It actually allows to avoid possible variables miss-spelling: i.e. ‘Bphi’=’bp’=’Bp’=’bphi’

Parameters:
  • obj (magic.MagicGraph) – a graphic output file

  • field (str) – the name of the field one wants to select

  • labTex (bool) – when set to True, format the labels using LaTeX fonts

Returns:

a tuple that contains the selected physical field and its label

Return type:

(numpy.ndarray, str)

magic.libmagic.symmetrize(data, ms, reversed=False)[source]

Symmetrise an array which is defined only with an azimuthal symmetry minc=ms

Parameters:
  • data (numpy.ndarray) – the input array

  • ms (int) – the azimuthal symmetry

  • reversed (bool) – set to True, in case the array is reversed (i.e. n_phi is the last column)

Returns:

an output array of dimension (data.shape[0]*ms+1)

Return type:

numpy.ndarray

magic.libmagic.thetaderavg(data, order=4)[source]

Theta-derivative of an input array (finite differences)

>>> gr = MagiGraph()
>>> dvtdt = thetaderavg(gr.vtheta)
Parameters:
  • data (numpy.ndarray) – input array

  • order (int) – order of the finite-difference scheme (possible values are 2 or 4)

Returns:

the theta-derivative of the input array

Return type:

numpy.ndarray

magic.libmagic.timeder(time, y)[source]

time derivative of an input array

computed with central differences (numpy.gradient)

>>> ts = MagicTs(field='e_kin')
>>> dt_ekinpol = timeder(ts,field='ekin_pol')
magic.libmagic.writeVpEq(par, tstart)[source]

This function computes the time-averaged surface zonal flow (and Rolc) and format the output

>>> # Reads all the par.* files from the current directory
>>> par = MagicTs(field='par', iplot=False, all=True)
>>> # Time-average
>>> st = writeVpEq(par, tstart=2.1)
>>> print(st)
Parameters:
  • par (magic.MagicTs) – a MagicTs object containing the par file

  • tstart (float) – the starting time of the averaging

Returns:

a formatted string

Return type:

str

magic.libmagic.zderavg(data, rad, colat=None, exclude=False)[source]

z derivative of an input array

>>> gr = MagiGraph()
>>> dvrdz = zderavg(gr.vr, eta=gr.radratio, colat=gr.colatitude)
Parameters:
  • data (numpy.ndarray) – input array

  • rad (numpy.ndarray) – radial grid

  • exclude (bool) – when set to True, exclude the first and last radial grid points and replace them by a spline extrapolation (default is False)

  • colat (numpy.ndarray) – colatitudes (when not specified a regular grid is assumed)

Returns:

the z derivative of the input array

Return type:

numpy.ndarray