Compute SSS basis using non-optimized (but more readable) algorithms.
(exp, coils, mag_scale=100.0, method="standard")
| 1553 | |
| 1554 | |
| 1555 | def _sss_basis_basic(exp, coils, mag_scale=100.0, method="standard"): |
| 1556 | """Compute SSS basis using non-optimized (but more readable) algorithms.""" |
| 1557 | int_order, ext_order = exp["int_order"], exp["ext_order"] |
| 1558 | origin = exp["origin"] |
| 1559 | assert "extended_proj" not in exp # advanced option not supported |
| 1560 | # Compute vector between origin and coil, convert to spherical coords |
| 1561 | if method == "standard": |
| 1562 | # Get position, normal, weights, and number of integration pts. |
| 1563 | rmags, cosmags, ws, bins = _concatenate_coils(coils) |
| 1564 | rmags -= origin |
| 1565 | # Convert points to spherical coordinates |
| 1566 | rad, az, pol = _cart_to_sph(rmags).T |
| 1567 | cosmags *= ws[:, np.newaxis] |
| 1568 | del rmags, ws |
| 1569 | out_type = np.float64 |
| 1570 | else: # testing equivalence method |
| 1571 | rs, wcoils, ezs, bins = _concatenate_sph_coils(coils) |
| 1572 | rs -= origin |
| 1573 | rad, az, pol = _cart_to_sph(rs).T |
| 1574 | ezs *= wcoils[:, np.newaxis] |
| 1575 | del rs, wcoils |
| 1576 | out_type = np.complex128 |
| 1577 | del origin |
| 1578 | |
| 1579 | # Set up output matrices |
| 1580 | n_in, n_out = _get_n_moments([int_order, ext_order]) |
| 1581 | S_tot = np.empty((len(coils), n_in + n_out), out_type) |
| 1582 | S_in = S_tot[:, :n_in] |
| 1583 | S_out = S_tot[:, n_in:] |
| 1584 | coil_scale = np.ones((len(coils), 1)) |
| 1585 | coil_scale[_get_mag_mask(coils)] = mag_scale |
| 1586 | |
| 1587 | # Compute internal/external basis vectors (exclude degree 0; L/RHS Eq. 5) |
| 1588 | for degree in range(1, max(int_order, ext_order) + 1): |
| 1589 | # Only loop over positive orders, negative orders are handled |
| 1590 | # for efficiency within |
| 1591 | for order in range(degree + 1): |
| 1592 | S_in_out = list() |
| 1593 | grads_in_out = list() |
| 1594 | # Same spherical harmonic is used for both internal and external |
| 1595 | sph = sph_harm_y(degree, order, pol, az) |
| 1596 | sph_norm = _sph_harm_norm(order, degree) |
| 1597 | # Compute complex gradient for all integration points |
| 1598 | # in spherical coordinates (Eq. 6). The gradient for rad, az, pol |
| 1599 | # is obtained by taking the partial derivative of Eq. 4 w.r.t. each |
| 1600 | # coordinate. |
| 1601 | az_factor = 1j * order * sph / np.sin(np.maximum(pol, 1e-16)) |
| 1602 | pol_factor = ( |
| 1603 | -sph_norm |
| 1604 | * np.sin(pol) |
| 1605 | * np.exp(1j * order * az) |
| 1606 | * _alegendre_deriv(order, degree, np.cos(pol)) |
| 1607 | ) |
| 1608 | if degree <= int_order: |
| 1609 | S_in_out.append(S_in) |
| 1610 | in_norm = _mu_0 * rad ** -(degree + 2) |
| 1611 | g_rad = in_norm * (-(degree + 1.0) * sph) |
| 1612 | g_az = in_norm * az_factor |