Add an SDE solver using the "LogWright" function

[kandersolar]

Is your feature request related to a problem? Please describe.
The lambertw method has to take some pains to appropriately handle numerical overflow in an intermediate calculation.

Describe the solution you'd like
An implementation of [1,2], which reformulates the SDE so that it uses a different form of the Lambert W function and does not require calculating such large intermediate values.

Additionally, according to [3] and also my own brief testing, it is also computationally faster than the lambertw method.

Additional context
[1] "On Calculating the Current-Voltage Characteristic of Multi-Diode Models for Organic Solar Cells": https://arxiv.org/abs/1601.02679
[2] "A Robust Approximation to a Lambert-Type Function": https://arxiv.org/abs/1504.01964
[3] "Solar Cells, Lambert W and the LogWright Functions": https://arxiv.org/abs/2307.08099

(I notice also that @markcampanelli is mentioned in the acknowledgements of [2])

Here is a rough implementation of the method, which I observe to be 2x faster than _lambertw_v_from_i while staying within $\pm$ 1e-10 of its output.

def _g(x, n=2):
    y = np.where(x <= -np.e, x, 1)
    y = np.where((-np.e < x) & (x < np.e), -np.e + (1 + np.e) * (x + np.e) / (2 * np.e), y)
    np.log(x, out=y, where=x >= np.e)
    
    for _ in range(n):
        ey = np.exp(y)
        ey_plus_one = 1 + ey
        y_ey_x = y + ey - x
        y = y - 2 * (y_ey_x) * (ey_plus_one) / ((2 * (ey_plus_one)**2 - (y_ey_x)*ey))
    
    return y


def _logwright_v_from_i(current, photocurrent, saturation_current,
                        resistance_series, resistance_shunt, nNsVth):
    output_is_scalar = all(map(np.isscalar,
                               (current, photocurrent, saturation_current,
                                resistance_series, resistance_shunt, nNsVth)))

    I, IL, I0, Rs, Rsh, a = \
        np.broadcast_arrays(current, photocurrent, saturation_current,
                            resistance_series, resistance_shunt, nNsVth)

    x = np.log(I0 * Rsh / a) + Rsh * (-I + IL + I0) / a
    V = -I * Rs - a * np.log(I0 * Rsh / a) + a * _g(x)

    if output_is_scalar:
        return V.item()
    else:
        return V

[markcampanelli]

@kandersolar Oh I fondly recall my brief interaction with Ken Roberts. I'm glad he and others have been able to develop his idea, and I didn't know it has apparently been extended to 2-diode models.

It seems like the LogWright function could be a good addition to scipy's special functions, perhaps alongside https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.wrightomega.html#scipy.special.wrightomega