Add functionality for a custom outer boundary condition on `br`

This PR adds the option to solve using an arbitrary map for the outer boundary condition on br, a functionality which is already present in Anthony Yeates's original code. Here, it is implemented in a similar manner. (For mathematical details, see the text following Equation 24 in the "Notes on the Numerical Methods" document.)

An example use would be as follows:

import sunpy.map

gong_fname = get_gong_map()
gong_map = sunpy.map.Map(gong_fname)
nrho = 16
rss = 2.5
pfss_in = pfsspy.Input(gong_map, nrho, rss, brout="closed", br1=gong_map.data)
pfss_out = pfsspy.pfss(pfss_in)

This should give a solution where the outer boundary condition is the same source field as the inner. Alternatively, allowing br1=0*gong_map.data should give an outer boundary condition where the source field Br is always zero. This latter example may be useful to some, e.g., Caplan et al 2021.

dstansby / pfsspy

Add functionality for a custom outer boundary condition on `br` #357