# Spherical harmonic comparisons

Comparing analytical spherical harmonic solutions to PFSS output.

First, import required modules

```import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
from matplotlib.gridspec import GridSpec

import pfsspy
from pfsspy import analytic
```

Setup some useful functions for testing

```def theta_phi(nphi, ns):
# Return a theta, phi grid with a given numer of points
phi = np.linspace(0, 2 * np.pi, nphi)
s = np.linspace(-1, 1, ns)
s, phi = np.meshgrid(s, phi)
theta = np.arccos(s)

# Return the pfsspy solution for given input parameters
theta, phi = theta_phi(nphi, ns)

br_in = analytic.Br(l, m, rss)(1, theta, phi)

pfss_output = pfsspy.pfss(pfss_input)
return pfss_output.bc[0][:, :, -1].T

# Return the analytic solution for given input parameters
theta, phi = theta_phi(nphi, ns)
```

Compare the the pfsspy solution to the analytic solutions. Cuts are taken on the source surface at a constant phi value to do a 1D comparison.

```ls = [1, 2, 3]
fig = plt.figure(tight_layout=True)
gs = GridSpec(len(ls), len(ls) + 1)

for i, l in enumerate(ls):
axs = [fig.add_subplot(gs[i, j], sharey=ax0) for j in range(1, l + 1)]
axs = [ax0] + axs
for j, m in enumerate(list(range(l + 1))):
ax = axs[j]
nphi = 359
ns = 179
nrho = 20