1D RM-Synthesis¶

[1]:

from __future__ import annotations

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.visualization import quantity_support
from rm_lite.tools_1d import rmsynth

plt.rcParams["figure.dpi"] = 150

_ = quantity_support()
rng = np.random.default_rng(42)

First generate some synthetic data

[2]:

from rm_lite.utils.simulate import faraday_simple_spectrum, faraday_slab_spectrum
from rm_lite.utils.synthesis import freq_to_lambda2

Here we’ll simulate RACS-all frequency coverage

[3]:

bw_low = 288
bw_mid = 144
bw_high = 288
low = np.linspace(943.5 - bw_low / 2, 943.5 + bw_low / 2, 36) * u.MHz
mid = np.linspace(1367.5 - bw_mid / 2, 1367.5 + bw_mid / 2, 9) * u.MHz
high = np.linspace(1655.5 - bw_high / 2, 1655.5 + bw_high / 2, 9) * u.MHz
freqs = np.concatenate([low, mid, high])
freq_hz = freqs.to(u.Hz).value

Now we make a Faraday simple spectrum with a single RM component. We will use the following parameters:

[4]:

delta_rm_radm2 = 30
rm_radm2 = 100
frac_pol = 0.5
psi0_deg = 10
complex_data_noiseless = faraday_simple_spectrum(
    freq_to_lambda2(freq_hz),
    frac_pol=frac_pol,
    psi0_deg=psi0_deg,
    rm_radm2=rm_radm2,
)

[5]:

fig, ax = plt.subplots()
ax.plot(
    freq_hz, np.real(complex_data_noiseless), ".", label="Stokes Q", color="tab:red"
)
ax.plot(
    freq_hz, np.imag(complex_data_noiseless), ".", label="Stokes U", color="tab:blue"
)
ax.legend()
ax.set(
    xlabel=rf"$\nu$ / {u.Hz:latex_inline}",
    ylabel="Flux density",
    title="Stokes Q and U",
)

[5]:

[Text(0.5, 0, '$\\nu$ / $\\mathrm{Hz}$'),
 Text(0, 0.5, 'Flux density'),
 Text(0.5, 1.0, 'Stokes Q and U')]

Now we can run RM-synthesis by calling rmsynth.run_rmsynth

[6]:

help(rmsynth.run_rmsynth)

Help on function run_rmsynth in module rm_lite.tools_1d.rmsynth:

run_rmsynth(freq_arr_hz: 'NDArray[np.float64]', complex_pol_arr: 'NDArray[np.complex128]', complex_pol_error: 'NDArray[np.complex128]', stokes_i_arr: 'NDArray[np.float64] | None' = None, stokes_i_error_arr: 'NDArray[np.float64] | None' = None, stokes_i_model_arr: 'NDArray[np.float64] | None' = None, stokes_i_model_error: 'NDArray[np.float64] | None' = None, phi_max_radm2: 'float | None' = None, d_phi_radm2: 'float | None' = None, n_samples: 'float | None' = 10.0, weight_type: 'WeightType' = 'variance', robust: 'float | None' = None, do_fit_rmsf: 'bool' = False, do_fit_rmsf_real: 'bool' = False, fit_function: "Literal['log', 'linear']" = 'log', fit_order: 'int' = 2, ignore_stokes_i: 'bool' = False, moment_threshold_snr: 'float' = 5.0) -> 'RMSynth1DResults'
    Run RM-synthesis on 1D data

    Args:
        freq_arr_hz (NDArray[np.float64]): Frequencies in Hz
        complex_pol_arr (NDArray[np.complex128]): Complex polarisation values (Q + iU)
        complex_pol_error (NDArray[np.float64]): Complex polarisation errors (dQ + idU)
        stokes_i_arr (NDArray[np.float64] | None, optional): Total itensity values. Defaults to None.
        stokes_i_error_arr (NDArray[np.float64] | None, optional): Total intensity errors. Defaults to None.
        stokes_i_model_arr (NDArray[np.float64] | None, optional): Total intensity model array. Defaults to None.
        stokes_i_model_error (NDArray[np.float64] | None, optional): Total intensity model error. Defaults to None.
        phi_max_radm2 (float | None, optional): Maximum Faraday depth. Defaults to None.
        d_phi_radm2 (float | None, optional): Spacing in Faraday depth. Defaults to None.
        n_samples (float | None, optional): Number of samples across the RMSF. Defaults to 10.0.
        weight_type (WeightType, optional): Weighting: 'variance' (1/sigma^2), 'uniform' (equal per channel), 'uniform_lsq' (equal per lambda^2 interval, narrows the RMSF), 'briggs' (robust). Defaults to "variance".
        robust (float | None, optional): Briggs robust parameter, required for weight_type='briggs'. Defaults to None.
        do_fit_rmsf (bool, optional): Fit the RMSF main lobe. Defaults to False.
        do_fit_rmsf_real (bool, optional): Fit only the real part of the RMSF. Defaults to False.
        fit_function ("log" | "linear", optional): RMSF fit function. Defaults to "log".
        fit_order (int, optional): Polynomial fit order. Defaults to 2. Negative values will iterate until the fit is good.
        moment_threshold_snr (float, optional): SNR cut (times the theoretical FDF noise) applied to FDF amplitudes before computing the Faraday moments. Defaults to 5.0.

    Returns:
        RMSynth1DResults:
            fdf_parameters (pl.DataFrame): FDF parameters
            fdf_arrs (pl.DataFrame): RMSynth arrays
            rmsf_arrs (pl.DataFrame): RMSF arrays

[7]:

fdf_parameters, fdf_arrs, rmsf_arrs, stokes_i_arrs = rmsynth.run_rmsynth(
    freq_arr_hz=freq_hz,
    complex_pol_arr=complex_data_noiseless,
    complex_pol_error=np.zeros_like(complex_data_noiseless),
    do_fit_rmsf=True,
    n_samples=100,
)

WARNING rmsynth.run_rmsynth: Stokes I array/errors or model not provided. No fractional polarization will be calculated.
INFO synthesis.rmsynth_nufft: Running RM-synthesis using the NUFFTs over 2001 Faraday depth channels.
INFO synthesis.rmsynth_nufft: NUFFT complete in 0.0043 seconds.
INFO synthesis.get_rmsf_nufft: Fitting main lobe in each RMSF spectrum.
INFO rmsynth._run_rmsynth: RM-synthesis completed in 14.33ms.

The output values are Polars dataframes that can be inspected easily

[8]:

fdf_parameters

[8]:

shape: (1, 36)

fdf_error_mad	peak_pi_fit	peak_pi_error	peak_pi_fit_debias	peak_pi_fit_snr	peak_pi_fit_index	peak_rm_fit	peak_rm_fit_error	peak_q_fit	peak_u_fit	peak_pa_fit_deg	peak_pa_fit_deg_error	peak_pa0_fit_deg	peak_pa0_fit_deg_error	fit_function	lam_sq_0_m2	ref_freq_hz	fwhm_rmsf_radm2	phi_max_scale_radm2	fdf_error_noise	fdf_q_noise	fdf_u_noise	min_freq_hz	max_freq_hz	n_channels	median_d_freq_hz	frac_pol	frac_pol_error	sigma_add	sigma_add_minus	sigma_add_plus	mom0	mom0_debias	mom1_radm2	mom2_radm2	moment_threshold_snr
f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	str	f64	f64	f64	f64	f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64
NaN	0.500279	0.0	0.500279	inf	1296	99.999904	0.0	-0.202569	-0.457124	123.050018	0.0	10.000452	0.0	"log"	0.082563	1.0433e9	34.522155	113.191071	0.0	0.0	0.0	7.995e8	1.7995e9	54	8.2286e6	NaN	NaN	NaN	NaN	NaN	1.346747	1.346747	69.386287	138.885199	5.0

[9]:

fdf_arrs

[9]:

shape: (2_001, 2)

phi_arr_radm2	fdf_dirty_complex_arr
f64	object
-336.725921	(-0.012941952987807103+0.005575905271668905j)
-336.389195	(-0.013584447466549732+0.005189884663585659j)
-336.052469	(-0.014213132463713843+0.004794461009568579j)
-335.715743	(-0.01482712771337132+0.004389828120960796j)
-335.379017	(-0.015425552941836488+0.003976168117102392j)
…	…
335.379017	(-0.0009008057626928257+0.008053597616970908j)
335.715743	(-0.0009576261370891269+0.007885479952290737j)
336.052469	(-0.0009822552939409935+0.007717971936861509j)
336.389195	(-0.0009751402168175126+0.007551392504192875j)
336.725921	(-0.0009367721104613681+0.007386049740302505j)

[10]:

rmsf_arrs

[10]:

shape: (4_003, 2)

phi2_arr_radm2	rmsf_complex_arr
f64	object
-673.788568	(-0.0329575974415113-0.0360831117177476j)
-673.451842	(-0.032467101832120024-0.03631435677037289j)
-673.115116	(-0.03194042093852627-0.036538602915399204j)
-672.77839	(-0.03137760921834162-0.03675582151570247j)
-672.441664	(-0.030778770230473887-0.03696597750309249j)
…	…
672.441664	(-0.030778770230891164+0.036965977502953294j)
672.77839	(-0.031377609218734684+0.036755821515558944j)
673.115116	(-0.031940420938894855+0.03653860291525102j)
673.451842	(-0.03246710183246473+0.03631435677022105j)
673.788568	(-0.03295760578020578+0.03608308730681604j)

Since we provided no Stokes \(I\) data, the stokes I model will just be unity with 0 error. The flag_arr array tells us which channels were not used in RM-synthesis or model fitting

[11]:

stokes_i_arrs

[11]:

shape: (54, 7)

freq_arr_hz	lambda_sq_arr_m2	stokes_i_model_arr	stokes_i_model_error	flag_arr	complex_pol_arr	complex_pol_error
f64	f64	f64	f64	bool	object	object
7.995e8	0.140606	null	null	true	(-0.49042944532403415-0.09735994638022451j)	0j
8.0773e8	0.137756	null	null	true	(-0.4654241884816873+0.18270283187778674j)	0j
8.1596e8	0.134992	null	null	true	(-0.3001379571186853+0.3998964949791661j)	0j
8.2419e8	0.13231	null	null	true	(-0.053617325136419355+0.497116870006657j)	0j
8.3241e8	0.129707	null	null	true	(0.20074123829142537+0.45793335240974226j)	0j
…	…	…	…	…	…	…
1.6555e9	0.032793	null	null	true	(0.4056251419940925+0.29234952399871j)	0j
1.6915e9	0.031412	null	null	true	(0.46997598588767076+0.17065336998990638j)	0j
1.7275e9	0.030117	null	null	true	(0.49801253407054746+0.04453668048509923j)	0j
1.7635e9	0.0289	null	null	true	(0.49406601726750204-0.07680345422849415j)	0j
1.7995e9	0.027755	null	null	true	(0.463743096038396-0.1869287053309979j)	0j

We can also easily visualise the data

[12]:

phi_arr_radm2 = fdf_arrs["phi_arr_radm2"].to_numpy()
fdf_dirty_arr = fdf_arrs["fdf_dirty_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
x1, x2, y1, y2 = 95, 105, 0.45, 0.55  # subregion of the original image
axins = ax.inset_axes(
    (0.9, 0.6, 0.4, 0.4), xlim=(x1, x2), ylim=(y1, y2), xticklabels=[], yticklabels=[]
)
for _ax in [ax, axins]:
    _ax.plot(
        phi_arr_radm2,
        fdf_dirty_arr.real,
        color="tab:red",
        label="Stokes Q",
    )
    _ax.plot(
        phi_arr_radm2,
        fdf_dirty_arr.imag,
        color="tab:blue",
        label="Stokes U",
    )
    _ax.plot(
        phi_arr_radm2,
        np.abs(fdf_dirty_arr),
        color="k",
        label="Polarized intensity",
    )

    _ax.errorbar(
        fdf_parameters["peak_rm_fit"],
        fdf_parameters["peak_pi_fit"],
        xerr=fdf_parameters["peak_rm_fit_error"],
        yerr=fdf_parameters["peak_pi_error"],
        fmt="o",
        lw=1,
        color="red",
        mfc="none",
        label="Fitted peak",
    )

ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="Flux density",
    title="Dirty FDF",
    # xlim=[50, 150],
)
ax.indicate_inset_zoom(axins, edgecolor="black")
ax.legend()

[12]:

<matplotlib.legend.Legend at 0x73b3ccd28fb0>

[13]:

phi2_arr_radm2 = rmsf_arrs["phi2_arr_radm2"].to_numpy()
rmsf_arr = rmsf_arrs["rmsf_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi2_arr_radm2,
    np.abs(rmsf_arr),
    color="k",
    label="Polarized intensity",
)
ax.legend()
ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="RMSF",
    title="RMSF",
)

[13]:

[Text(0.5, 0, '$\\phi$ / $\\mathrm{rad\\,m^{-2}}$'),
 Text(0, 0.5, 'RMSF'),
 Text(0.5, 1.0, 'RMSF')]

Now lets do a more complex example. We’ll add noise and a Stokes \(I\) spectrum

[14]:

from rm_lite.utils.fitting import power_law

[15]:

delta_rm_radm2 = 30
rm_radm2 = 100
frac_pol = 0.5
psi0_deg = 10
complex_data_noiseless = faraday_slab_spectrum(
    freq_to_lambda2(freq_hz),
    frac_pol=frac_pol,
    psi0_deg=psi0_deg,
    rm_radm2=rm_radm2,
    delta_rm_radm2=delta_rm_radm2,
)


stokes_i_flux = 1.0
spectral_index = -0.7
rms_noise = 0.1


stokes_i_model = power_law(order=1)
stokes_i_noiseless = stokes_i_model(
    freq_hz / (np.mean(freq_hz)), stokes_i_flux, spectral_index
)
stokes_i_noise = rng.normal(0, rms_noise, size=freq_hz.size)
stokes_i_noisy = stokes_i_noiseless + stokes_i_noise


stokes_q_noise = rng.normal(0, rms_noise, size=freq_hz.size)
stokes_u_noise = rng.normal(0, rms_noise, size=freq_hz.size)
complex_noise = stokes_q_noise + 1j * stokes_u_noise

complex_flux = complex_data_noiseless * stokes_i_noiseless
complex_data_noisy = complex_data_noiseless + complex_noise

Now we enable Stokes \(I\) model fitting through providing the data, and enabling fit_order. If fit_order<0 an iterative fit will be performed.

[16]:

fdf_parameters, fdf_arrs, rmsf_arrs, stokes_i_arrs = rmsynth.run_rmsynth(
    freq_arr_hz=freq_hz,
    complex_pol_arr=complex_data_noisy,
    complex_pol_error=np.ones_like(complex_data_noiseless)
    * (rms_noise + rms_noise * 1j),
    stokes_i_arr=stokes_i_noisy,
    stokes_i_error_arr=np.ones_like(stokes_i_noisy) * rms_noise,
    do_fit_rmsf=True,
    n_samples=100,
    fit_order=-3,
)

INFO synthesis.create_fractional_spectra: Fitting Stokes I model to calculate fractional spectra.
INFO fitting.dynamic_fit: Iteratively fitting Stokes I model of type log with max order 3.
INFO fitting.static_fit: Fitting Stokes I model of type log with order 0.
INFO fitting.static_fit: Fit results: ['1.04 +/- 0.0136']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 1.
INFO fitting.static_fit: Fit results: ['1.07 +/- 0.0138', '-0.671 +/- 0.0608']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 2.
INFO fitting.static_fit: Fit results: ['1.08 +/- 0.0205', '-0.62 +/- 0.0835', '-0.593 +/- 0.672']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 3.
INFO fitting.static_fit: Fit results: ['1.08 +/- 0.0212', '-0.62 +/- 0.127', '-0.589 +/- 1.1', '-0.0254 +/- 6.9']
INFO fitting.dynamic_fit: Fit results for orders [0 1 2 3]:
INFO fitting.dynamic_fit: Best fit found with 2 parameters.
INFO synthesis.rmsynth_nufft: Running RM-synthesis using the NUFFTs over 2001 Faraday depth channels.
INFO synthesis.rmsynth_nufft: NUFFT complete in 0.00204 seconds.
INFO synthesis.get_rmsf_nufft: Fitting main lobe in each RMSF spectrum.
INFO rmsynth._run_rmsynth: RM-synthesis completed in 7.53ms.

[17]:

fdf_parameters

[17]:

shape: (1, 36)

fdf_error_mad	peak_pi_fit	peak_pi_error	peak_pi_fit_debias	peak_pi_fit_snr	peak_pi_fit_index	peak_rm_fit	peak_rm_fit_error	peak_q_fit	peak_u_fit	peak_pa_fit_deg	peak_pa_fit_deg_error	peak_pa0_fit_deg	peak_pa0_fit_deg_error	fit_function	lam_sq_0_m2	ref_freq_hz	fwhm_rmsf_radm2	phi_max_scale_radm2	fdf_error_noise	fdf_q_noise	fdf_u_noise	min_freq_hz	max_freq_hz	n_channels	median_d_freq_hz	frac_pol	frac_pol_error	sigma_add	sigma_add_minus	sigma_add_plus	mom0	mom0_debias	mom1_radm2	mom2_radm2	moment_threshold_snr
f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	str	f64	f64	f64	f64	f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64
NaN	0.204969	0.014589	0.203771	14.049493	1265	89.365093	1.228591	-0.198433	-0.050471	97.135175	2.039069	34.393588	11.803601	"log"	0.082563	1.0433e9	34.522155	113.191071	0.014589	0.014608	0.01457	7.995e8	1.7995e9	54	8.2286e6	0.191143	0.013685	2.054997	1.586183	2.535694	0.624902	0.615042	115.552737	58.845317	5.0

[18]:

fdf_arrs

[18]:

shape: (2_001, 2)

phi_arr_radm2	fdf_dirty_complex_arr
f64	object
-336.725921	(-0.014738173003141175+0.009372052470022232j)
-336.389195	(-0.015269780690890756+0.009370568081949573j)
-336.052469	(-0.015794276454245474+0.009361398402216321j)
-335.715743	(-0.016310968330576707+0.009344198836196797j)
-335.379017	(-0.016819164175972918+0.00931860431917085j)
…	…
335.379017	(-0.020917865122834278+0.03205000232658517j)
335.715743	(-0.021469957869800993+0.03132128925890664j)
336.052469	(-0.02199789950551498+0.030593761735830328j)
336.389195	(-0.022501544578229376+0.029868651890150415j)
336.725921	(-0.022980778790561224+0.029147173743289064j)

[19]:

stokes_i_arrs

[19]:

shape: (54, 7)

freq_arr_hz	lambda_sq_arr_m2	stokes_i_model_arr	stokes_i_model_error	flag_arr	complex_pol_arr	complex_pol_error
f64	f64	f64	f64	bool	object	object
7.995e8	0.140606	1.275069	0.052855	true	(0.04338330037782169+0.08183877539714739j)	(0.07844773110714742+0.07850045294804807j)
8.0773e8	0.137756	1.266291	0.051246	true	(0.023954174896411554-0.16561639003877818j)	(0.07897675329112609+0.07925471095100153j)
8.1596e8	0.134992	1.257689	0.049497	true	(0.024580042563320165-0.02735655938948127j)	(0.0795167928357659+0.07951819725215176j)
8.2419e8	0.13231	1.249224	0.048165	true	(0.1276321078739654-0.0547943844919334j)	(0.08020080318272964+0.08007756462558462j)
8.3241e8	0.129707	1.240965	0.046527	true	(-0.09809548941561544-0.11249267088983378j)	(0.08066635817231832+0.08069277040078873j)
…	…	…	…	…	…	…
1.6555e9	0.032793	0.781665	0.04704	true	(0.5335278716959772+0.44225475096806127j)	(0.13189949646253907+0.1306710777956422j)
1.6915e9	0.031412	0.770364	0.048235	true	(0.40250870448586556+0.12534835081582515j)	(0.13223271970388215+0.13004589640679562j)
1.7275e9	0.030117	0.759361	0.048891	true	(0.5431221781032378-0.10505460325026317j)	(0.13625323808315237+0.13186320722625705j)
1.7635e9	0.0289	0.748901	0.049887	true	(0.45326910551159816-0.21905252847583304j)	(0.13690017582201144+0.13432394362881442j)
1.7995e9	0.027755	0.73907	0.05067	true	(0.511564507931407-0.3228225475746143j)	(0.13977689691222428+0.1371034334406259j)

[20]:

fig, ax = plt.subplots()
ax.plot(freq_hz, stokes_i_noiseless, label="Input model")
ax.plot(freq_hz, stokes_i_noisy, ".", label="Noisy data")
ax.plot(
    stokes_i_arrs["freq_arr_hz"],
    stokes_i_arrs["stokes_i_model_arr"],
    "k--",
    label="Fitted model",
)
ax.fill_between(
    stokes_i_arrs["freq_arr_hz"],
    stokes_i_arrs["stokes_i_model_arr"] - stokes_i_arrs["stokes_i_model_error"],
    stokes_i_arrs["stokes_i_model_arr"] + stokes_i_arrs["stokes_i_model_error"],
    alpha=0.3,
    color="k",
    label="Fitted model error",
)
ax.legend()
ax.set(
    xlabel=rf"$\nu$ / {u.Hz:latex_inline}",
    ylabel="Flux density",
    title="Stokes I",
)

[20]:

[Text(0.5, 0, '$\\nu$ / $\\mathrm{Hz}$'),
 Text(0, 0.5, 'Flux density'),
 Text(0.5, 1.0, 'Stokes I')]

[21]:

phi_arr_radm2 = fdf_arrs["phi_arr_radm2"].to_numpy()
fdf_dirty_arr = fdf_arrs["fdf_dirty_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()

ax.plot(
    phi_arr_radm2,
    fdf_dirty_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi_arr_radm2,
    fdf_dirty_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi_arr_radm2,
    np.abs(fdf_dirty_arr),
    color="k",
    label="Polarized intensity",
)

ax.errorbar(
    fdf_parameters["peak_rm_fit"],
    fdf_parameters["peak_pi_fit"],
    xerr=fdf_parameters["peak_rm_fit_error"],
    yerr=fdf_parameters["peak_pi_error"],
    fmt="o",
    lw=1,
    color="red",
    mfc="none",
    label="Fitted peak",
)

ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="Flux density",
    title="Dirty FDF",
)
ax.legend()

[21]:

<matplotlib.legend.Legend at 0x73b3cc1380e0>

[22]:

phi2_arr_radm2 = rmsf_arrs["phi2_arr_radm2"].to_numpy()
rmsf_arr = rmsf_arrs["rmsf_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi2_arr_radm2,
    np.abs(rmsf_arr),
    color="k",
    label="Polarized intensity",
)
ax.legend()
ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="RMSF",
    title="RMSF",
)

[22]:

[Text(0.5, 0, '$\\phi$ / $\\mathrm{rad\\,m^{-2}}$'),
 Text(0, 0.5, 'RMSF'),
 Text(0.5, 1.0, 'RMSF')]