Spectral tansforms

class magic.spectralTransforms.SpectralTransforms(l_max=32, minc=1, n_theta_max=64, m_max=None, verbose=True)

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, 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)

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)

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)

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)

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)

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)

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)

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)

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=None, vmax=None, vmin=None, cbar=True, title=True, normRad=False, deminc=True, bounds=True, lines=False, linewidths=0.5, pcolor=False, rasterized=False)

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.)

• title (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 None, it tries to find it by itself.

• 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)

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

• linewidths (float) – the thickness of the solid lines, whenever plotted

• pcolor (bool) – when set to True, use pcolormesh instead of contourf

• rasterized (bool) – when set to True, the rasterization for vector graphics is turned on

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

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=None, vmax=None, vmin=None, cbar=True, title=True, fig=None, ax=None, bounds=True, lines=False, pcolor=False, linewidths=0.5, rasterized=False)

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.)

• title (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 None, it tries to guess by itself.

• 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

• linewidths (float) – the thickness of the solid lines, whenever plotted

• pcolor (bool) – when set to True, use pcolormesh instead of contourf

• rasterized (bool) – when set to True, the rasterization for vector graphics is turned on

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=None, cbar=True, title=True, lines=False, fig=None, ax=None, linewidths=0.5, pcolor=False, rasterized=False, gridLineStyle=':', gridColor='k', gridLineWidth=0.7)

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.)

• title (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

• linewidths (float) – the thickness of the solid lines, whenever plotted

• 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 None, it tries to find it by itself.

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

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

• pcolor (bool) – when set to True, use pcolormesh instead of contourf

• rasterized (bool) – when set to True, the rasterization for vector graphics is turned on

• gridColor (str) – this is used to set the color of the grid

• gridLineStyle (str) – this allows to set the line style of the grid (‘:’, ‘-‘, ‘–’)

• gridLineWidth (float) – this is used to tune the thickness of the lines used in the grid

Various useful functions

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

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)

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, tstop=None)

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

• tstart – the stopping 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)

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.cleanBIS(dir='.')

This function renames all files with a trailing ‘_BIS’ in a given directory.

Parameters

dir (str) – the working directory

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

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)

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'>)

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)

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)

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()

This function calculates the total CPU time of one run directory

Returns

the total RUN time

Return type

float

magic.libmagic.horizontal_mean(field, colat, std=False)

This function computes the horizontal mean (and the standard deviation) of an input array of size (nphi,ntheta) or (nphi,ntheta,nr)

Parameters

• field (numpy.ndarray) – the input array

• colat (numpy.ndarray) – an array that contains the colatitudes

• std (bool) – a boolean if one also wants to compute the standard deviation

Returns

the average value or radial profile

Return type

float

magic.libmagic.intcheb(f, nr, z1, z2)

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)

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)

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)

This function returns all prime factors of a number

Returns

all prime factors

Return type

list

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

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)

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)

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)

s derivative of an input array

>>> gr = MagiGraph()
>>> dvpds = sderavg(gr.vphi, gr.radius, 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)

This routine computes the second time derivative of an input array

>>> ts = MagicTs(field='e_kin')
>>> d2Ekdt2 = secondtimeder(ts, ts.ekin_pol)

Parameters

• time (numpy.ndarray) – an array that contains time

• y (numpy.ndarray) – an array which contains the field to be differentiated two times

Returns

an array that contains the second time derivative

Return type

numpy.ndarray

magic.libmagic.selectField(obj, field, labTex=True, ic=False)

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)

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)

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)

This routine computes the time derivative of an input array

>>> ts = MagicTs(field='e_kin')
>>> dEkdt = timeder(ts, ts.ekin_pol)

Parameters

• time (numpy.ndarray) – an array that contains time

• y (numpy.ndarray) – an array which contains the field to be differentiated

Returns

an array that contains the time derivative

Return type

numpy.ndarray

magic.libmagic.writeVpEq(par, tstart)

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)

z derivative of an input array

>>> gr = MagiGraph()
>>> dvrdz = zderavg(gr.vr, gr.radius, 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