I'm trying to write a programme that draws an IV curve for a diode, using an expanded version of the Shockley diode equation that uses extra variables, but now the 'I' term appears on both sides of the equation and I can't work out how to rearrange it to do so. I understand that these things can be solved iteratively, but would prefer the rearrangement to make I the subject and only appear on one side.
$$I=I_0e^{\frac{q(V-IR_{se})}{nkT}}+\frac{V-IR_{se}}{R_{sh}} -I_L$$
Edit 26 Nov 2019
I don't have the old notebook, so
or, in usual notation and switching back to your capitalization and subscripting,
$$ I = \frac{k n T W\left(\frac{I_0 q R_{se} R_{sh} \mathrm{e}^{\frac{I_L q R_{se} R_{sh} + q R_{sh} V}{k n R_{se} T+k n R_{sh} T}}}{k n T (R_{se} + R_{sh})}\right)}{q R_{se}}+\frac{V-I_L R_{sh}}{R_{se}+R_{sh}} $$ which is slightly easier to read as
$$ I = \frac{V-I_L R_{sh}}{R_{se}+R_{sh}} + \frac{k n T}{q R_{se}} W\left(\frac{I_0 q R_{se} R_{sh} }{k n T (R_{se} + R_{sh})} \exp\left( \frac{I_L q R_{se} R_{sh} + q R_{sh} V}{k n R_{se} T+k n R_{sh} T}\right) \right) \text{.} $$
There seem to be typos in the version previously, probably related to trying to simplify using common subexpressions.
Let \begin{align*} a &= k q R_{sh} T \text{,} \\ b &= R_{se} + R_{sh} \end{align*} Then $$ I = \frac{-I_L R_{sh} + V}{b} + \frac{n}{k q R_{se} T} W\left( \frac{a \mathrm{e}^{a(I_L R_{se}+V)/bn}I_0 R_{se}}{bn} \right) \text{,} $$ where $W$ is the Lambert $W$ function.
Note that $W$ has several branches. You might notice that its graph does not satisfy the vertical line test. To get your entire graph, you may need to use different branches of $W$ in different places.
A little more about branches. The Lambert $W$ function has multiple branches (in the older literature, we would say it is a multivalued function). Before diving into this function, let's understand that a function having branches is not an unusual thing.
(Apologies in advance: I come from a math background, so I naturally use $\mathrm{i} = \sqrt{-1}$ instead of $\mathrm{j}$. Try not to think "current" when you see "$\mathrm{i}$" below. It should help that when I write the constant $\mathrm{i}$, it is upright, and when I write $i$ a variable frequently representing current, it is italicized.)
Consider $f(\theta) = \mathrm{e}^{\mathrm{i}\theta} = \cos(\theta) + \mathrm{i} \sin(\theta)$, where we have used Euler's formula to write the polar form of the complex number in rectangular form. As $\theta$ varies, this function runs around the unit circle in the complex plane. (Using rectangular coordinates for complex numbers, the real part of the number is the $x$-coordinate and the imaginary part is the $y$-coordinate. So to graph this, we would graph the points in the plane $(\cos \theta, \sin \theta)$, which trace out the unit circle.) Having said all that, notice that $$ \cdots = f(\theta - 2\pi) = f(\theta) = f(\theta + 2\pi) = f(\theta + 4\pi) = \cdots \text{,} $$ because sine and cosine are periodic with period $2\pi$. As exponentials, this says $$ \cdots = \mathrm{e}^{\mathrm{i}(\theta - 2\pi)} = \mathrm{e}^{\mathrm{i}\theta} = \mathrm{e}^{\mathrm{i}(\theta + 2\pi)} = \mathrm{e}^{\mathrm{i}(\theta + 4\pi)} = \cdots \text{.} $$ Let's take logarithms. The logarithm is the inverse function of the exponential, so this will be easy. (There is a small caveat I am delaying, so the first cou0ple of things I write won't be the best possible way to write them, because they are astoundingly wrong.) We see $$ \cdots = \mathrm{i}(\theta - 2\pi) = \mathrm{i}\theta = \mathrm{i}(\theta + 2\pi) = \mathrm{i}(\theta + 4\pi) = \cdots \text{.} $$ Cancel the common factor of $\mathrm{i}$ and subtract the common $\theta$ term and we discover $$ \cdots = - 2\pi = 0 = 2\pi = 4\pi = \cdots \text{.} $$ By now, we should have alarm bells going off, because we know that not all of these numbers are the same. (Certainly, $2\pi$ is not zero). What went wrong?
Recall what we noticed about $\mathrm{e}^{\mathrm{i}\theta}$: if we increase or decrease $\theta$ by a multiple of $2\pi$, the value of $\mathrm{e}^{\mathrm{i}(\theta + 2 \pi k)}$, where $k$ is any integer (positive, negative, or zero) is the same. So the logarithm cannot possibly start from a value given by evaluating a complex exponential and take us back to the angle that was used to generate it because many different angles give the same evaluation. The same thing happens when we try to construct inverse trigonometric functions. The trig functions take every value they take infinitely many times, so we conventionally restrict the range of angles that the inverse functions give back. In our setting, we might choose that the logarithm gives back the angle in one of the ranges $[-\pi, \pi)$, $(-\pi, \pi]$, $[0,2\pi)$, $(0, 2\pi]$, $[17, 17+2\pi)$, or any other half-open interval of length $2\pi$. (Actually, we could go crazy and for each angle in one circling of the origin, pick which multiple of $2\pi$ we will add to it to get the angle that the logarithm will produce. Typically, this is not very useful, but sometimes picking four multiples, one for each quadrant, is briefly useful.)
When we pick a range of angles for the logarithm, we are picking a branch of the logarithm. (Note that we are only focusing on the behaviour of the logarithm on the unit circle. In general, we would need to pick an angle interval to use for the whole complex plane, or pick a multiple of $2 \pi$ for each point in the plane. This can lead to some additional complexities, but these additional complexities are not critical for understanding branches of the Lambert $W$ function.) We can think of a set of spiral staircases running above and below the unit circle. As $\theta$ increases, we ascend, going anticlockwise along the stairs. As $\theta$ decreases, we descend, going clockwise along the stairs. Usually, we decide which range of altitudes for one time around stairs we are going to use for our angles. Every such choice is a choice of branch.
I used the logarithm above because the spiral staircase analogy is a really good one. The spiral staircase is called the Riemann surface for the logarithm (restricted to the unit circle). Other functions have different Riemann surfaces. For instance, the square root function has two "floors". (Warning, the upper end of the upper floor is glued to the lower end of the lower floor, so the image below necessarily has distortion -- you are looking at a three-dimensional shadow of something that does not intersect itself in four dimensions and naturally lives in four dimensions.)
It might help to notice that the real part axis is the bottom-right axis, the imaginary part axis is the bottom-left axis and the left axis is the real part of the square root of the complex point given by the other two coordinates. We see both square roots of $1$ on the right side of the graph, $\pm 1$. The two square roots of $-1$ are $\pm \mathrm{i}$, whose real part is zero, so this projection makes it look as if the two are the same, but they're not. Let's walk along these stairs. If we start at $(1,0,1)$, from the familiar fact that $\sqrt{1} = 1$ and walk anticlockwise along the surface (staying at radius $1$ in the input coordinates), we descend as we go around. Having circled the origin once, we are at the point $(1,0,-1)$ corresponding to the other square root of $1$. If we continue anticlockwise, we ascend back to $(1,0,1)$.
It should be familiar that the square root function cannot give you both square roots of a number. We conventionally have it give you only the nonnegative square root. That is, the convention is to use the upper half of this graph (all points above the horizontal plane $(x,y,0)$). Also, we decide that $\sqrt{-1} = \mathrm{i}$, so one copy of the (not really) self intersection is included and the other copy is not. By making this choice, we choose a branch of the square root function.
The short version is: when a function is not one-to-one (also called injective), its inverse function may have many choices for where to send a value back to. The various choices give different branches of the inverse function. To be a function (i.e., to have only one output for each input), we must pick a branch.
The Lambert $W$ function is an inverse function. The forward function is $f(z) = z \mathrm{e}^z$, giving $z = W(z \mathrm{e}^z)$. The forward function is not one-to-one. For instance, $$ f(18.9 + \mathrm{i}17.3) = f(15.290667\dots + \mathrm{i}946.398012\dots) \text{.} $$ (These tend to be transcendental numbers that are hard to write down except by mentioning branches of the $W$ function. The number on the right is $W_{151}(18.9 + \mathrm{i}17.3)$, where the branch number is the subscript.) If one computes $f$ at the two points above, using the digits shown, the values of $f$ agree in the seven most significant digits in both the real and imaginary parts. It is worth pointing out that $f$ magnifies its inputs, so small errors in the inputs produce larger errors in the outputs.
The Lambert $W$ function has infinitely many branches (so it is more like the logarithm than like the square root). Its Riemann surface is ... complicated. However, we do not have to understand the entire surface. It is conventional to make the $-1$ and $0$ branches be the only branches on which the $W$ function is real valued. For your purposes (computing real currents), you only care about these two branches. Putting these together, the graph of the part of the $W$ function having real output for real input is
We can immediately see that there are two branches. Since $f(z)$ for $z$ in $(-\infty, -1)$ gives the same set of values as $f(z)$ for $z$ in $(-1,0)$, we have two branches of $W$ to let us go back to either of those ranges of $z$. However, if $f(z) \geq 0$, we are on only one branch. The Mathematica numbering convention for these two branches is shown here:
So, using Mathematica, the $0^\text{th}$ branch of $W$ gives the currents you want. The NIST DLMF labels these as $Wp$ and $Wm$ and you would want $Wp$ for your current calculation. The SciPy function
lambertw()does not entirely specify its convention for all the branches of $W$, but it is likely that the default branch (withk=0in their description) is the one you would want.If your currents can go negative, you will want to calculate on both branches shown above to get a full picture of what is going on. This may mean that your system shows history dependent metastability. In particular for the quantity $q = \frac{a \mathrm{e}^{a(I_L R_{se}+V)/bn}I_0 R_{se}}{bn}$ near $-1/\mathrm{e}$, the system may randomly hop between the state where $I$ is given by $W_0$ and the state where $I$ is given by $W_{-1}$. If somehow $q = -1/\mathrm{e}$ is critically stable in your dynamics, when you lower $q$ to $-1/\mathrm{e}$, which branch of $W$ to use when raising it will be decided randomly at the instant of starting the increase and will control which branch of $W$ you use until you lower to $-1/\mathrm{e}$ and re-randomize.
I'm also a ham, so I have some idea what you are actually describing in your equations. Assuming I am interpreting your symbols as you intend, unless you are reverse biased, there is no way the argument to $W$ can drop below $0$. However, I can imagine that for certain very low resistances or very small $I_0$, $q$ can be very close to $0$ and in the presence of noise, $q$ might occasionally be negative, so you may see transient aggressive instability when $q$ is so small.