# Open flux

Comparing total unsigned flux to analytic solutions. This is done as a function of the number of radial grid cells used in the PFSS calculation.

First, import required modules

```import functools

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.integrate
from helpers import brss_pfsspy
from tqdm import tqdm

from pfsspy import analytic

def open_flux_analytic(l, m, zss):
Br = analytic.Br(l, m, zss)
Br = functools.partial(Br, zss)

def absBr(theta, phi):
return np.abs(Br(theta * u.rad, phi * u.rad)) * np.sin(theta)

res = scipy.integrate.nquad(absBr, ranges=[[0, np.pi], [0, 2 * np.pi]])
return res[0]

def open_flux_numeric(l, m, zss, nrho):
nphi = 360
ns = 180
br = brss_pfsspy(nphi, ns, nrho, zss, l, m)
return np.sum(np.abs(br)) * (4 * np.pi) / nphi / ns
```

Set the source surface height, and the (l, m) values to investigate

```zss = 2
lms = [10]

fig, ax = plt.subplots()
# nrhos = np.arange(10, 51, 2)
nrhos = [10, 50]

df = pd.DataFrame(index=nrhos, columns=lms)

for lm in lms:
l = lm // 10
m = lm % 10
print(f"l={l}, m={m}")
flux_analytic = open_flux_analytic(l, m, zss)
flux_numeric = []
for nrho in tqdm(nrhos):
flux_numeric.append(open_flux_numeric(l, m, zss, nrho))
flux_numeric = np.array(flux_numeric)
flux_ratio = flux_numeric / flux_analytic
df[lm] = flux_ratio
```

Save a copy of the data

```df.to_csv('open_flux_results.csv')
df = pd.read_csv('open_flux_results.csv', index_col=0)

for lm in lms:
l = lm // 10
m = lm % 10
color = {1: 'tab:blue', 2: 'tab:orange', 3: 'tab:green'}[l]
marker = {0: 'o', 1: 10, 2: 11}[m]
ax.plot(nrhos, df[str(lm)],
marker=marker, color=color,
label=f'l={l}, m={m}')

ax.legend()
ax.set_ylim(1)
ax.set_xlim(8)
ax.yaxis.grid(linestyle='--')
ax.set_xlabel('\$n_{r}\$')
ax.set_ylabel(r'\$\Phi_{pfsspy} / \Phi_{analytic}\$')
plt.show()
```

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery