Replicating Complex Plot (branch cuts, poles, shading)

Question

I am trying to recreate the following 2 plots:

where the function being plotted is $|e^{2ipR}S_0^{\text{well}}(p)|$ where $|...|$ is the absolute value, $R=1$, and $S_0^{\text{well}}(p)$ is given as \begin{equation} S_0^{\text{well}}(p) = e^{-2ipR}\frac{1+ipR\cdot\text{tanc}(R\sqrt{p^2-U})}{1+ipR\cdot\text{tanc}(R\sqrt{p^2-U})} \end{equation} where $\text{tanc}(x)\equiv \tan(x)/x$ The plots are for $U = 10$ and $U = -50$ (both at $R=1$).

Here is my attempt in Mathematica:

tanc[x_] := Tan[x]/x; 
S0well[p_, R_, U_] := Exp[-2*I*p*R]*((1 + I*p*R*tanc[R*Sqrt[p^2 - U]])/(1 - 
I*p*R*tanc[R*Sqrt[p^2 - U]])); 
ComplexPlot[Abs[Exp[2*I*p*1]*S0well[p, 1, -50]], {p, -10 - 10*I, 10 + 10*I}]

Which reproduces the following image:

Obviously, there are a few glaring issues. The first, the color. Second, why do I have branch cuts? I want to plot the branch cuts next but by converting the $S(p)$ function to one of $S(E)$ via $E = p^2/2m$ (by setting $2m = 1$). I am also missing the "blue" colored poles corresponding to zeros, but I think I have the "red dots" which are the poles. The images and equations are originally taken from section 2.1.2 (page 36) of The Analytic S-Matrix.

Any help and/or suggestions is appreciated.

Answer 1

I don't know how to remove the "brunch cut", but I've almost recreated the picture. Try DensityPlot instead:

DensityPlot[
 Log@Abs[Exp[2*I*p*1]*S0well[p, 1, -50]] /. p -> x + I y // 
  Evaluate, {x, -10, 10}, {y, -10, 10}, PlotPoints -> 50, 
 ColorFunction -> ColorData["DarkRainbow"]]

Update:

According to the comment from @user64494, now the following codes works:

DensityPlot[
 Log@Abs[Exp[2*I*p*1]*S0well[p, 1, -50]] /. p -> x + I y // 
  Evaluate, {x, -10, 10}, {y, -10, 10}, PlotPoints -> 100, 
 ColorFunction -> ColorData["DarkRainbow"], Exclusions -> None]

Or, we can compile the code to avoid the brunch cut handling and get faster speed:

fun = Compile[{x, y}, 
  Log@Abs[Exp[2*I*p*1]*S0well[p, 1, -50]] /. p -> x + I y // Evaluate];
DensityPlot[fun[x, y], {x, -10, 10}, {y, -10, 10}, PlotPoints -> 50, 
 ColorFunction -> ColorData["DarkRainbow"]]