from astropy.coordinates.matrix_utilities import rotation_matrix
from scipy.spatial.transform import Rotation
R_arr = rotation_matrix(30*u.deg, 'z')
R_scipy = Rotation.from_euler('z', 30, degrees=True)

fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharex=True, sharey=True)

grid = np.linspace(-5, 5, 100)

bar1 = gp.LongMuraliBarPotential(m=1e10, a=3.5, b=0.5, c=0.5,
                                 units=galactic)
bar2 = gp.LongMuraliBarPotential(m=1e10, a=3.5, b=0.5, c=0.5,
                                 units=galactic, R=R_arr)
bar3 = gp.LongMuraliBarPotential(m=1e10, a=3.5, b=0.5, c=0.5,
                                 units=galactic, R=R_scipy)

bar1.plot_contours(grid=(grid, grid, 0.), ax=axes[0])
bar2.plot_contours(grid=(grid, grid, 0.), ax=axes[1])
bar3.plot_contours(grid=(grid, grid, 0.), ax=axes[2])

axes[0].set_xlabel("$x$ [kpc]") # doctest: +SKIP
axes[0].set_ylabel("$y$ [kpc]") # doctest: +SKIP
axes[1].set_xlabel("$x$ [kpc]") # doctest: +SKIP
axes[2].set_xlabel("$x$ [kpc]") # doctest: +SKIP

fig.tight_layout()