Riemann surface of incomplete gamma function

Question

Michael Trott in his book The Mathematica GuideBook for Symbolics, p. 1004, has illustrated a nice way to visualize the Riemann surface of the incomplete gamma function $\Gamma(\alpha, z)$. To plot the Riemann surface of, for example, $\Gamma\left(\frac{3+i}{10}, z\right)$, we have:

With[{α = 0.3 + 0.1 I, ε = 10^-12},
     GraphicsRow[Graphics3D[{EdgeForm[Thickness[0.002]],
                             SurfaceColor[Hue[0.09], Hue[0.18], 2.3],
                             Table[Last /@ Partition[Cases[
                 ParametricPlot3D[{r Cos[φ], r Sin[φ],
                                   #[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                   (1 - Exp[2 k π I α]) Gamma[α]]},
                                  {r, 0, 2}, {φ, -π + ε, π - ε},
                                  PlotPoints -> {30, 40}], _Polygon, ∞], 2],
                                   {k, -2, 2}]}, BoxRatios -> {1, 1, 2.5},
                            PlotRange -> All] & /@ {Re, Im}]]

where we have used the identity:

$$\Gamma\left(\alpha, \exp(2 k \pi i)z\right) = \exp(2 k \pi i \alpha)\, \Gamma\left(\alpha, z\right) + (1- \exp(2 k \pi i \alpha)) \Gamma\left(\alpha\right)$$

However, Mathematica does not return anything.

Any help is appreciated!

Answer 1

I believe that code is obsolete and was written before the use of GraphicsComplex. It could be updated, but perhaps the easiest way to solve your problem is to simply remove the Cases operator (and other constructs) which incorrectly extracts the polygon data and just plot the 5 sheets of the function using the code below. Note that I explicitly extract the real part of the function with the Re operator. Can plot the imaginary with Im. You can then experiment with the code. For example, plot one sheet, let $k=0$ for example, or plot the sheets with different colors, etc.

α = 0.3 + 0.1 I;
ee = 10^(-12)
myTable = 
  Table[
    ParametricPlot3D[
      {r Cos[φ], 
       r Sin[φ], 
       Re[Exp[2 k Pi I α] Gamma[α, r Exp[I φ]] + (1 - Exp[2 k Pi I α]) Gamma[α]]}, 
      {r, 0, 2}, {φ, -Pi + ee, Pi - ee}, 
      PlotPoints -> {30, 40}, 
      BoxRatios -> {1, 1, 1}], {k, -2, 2}]; 
Show[myTable]

Answer 2

It seems to me that to fully reproduce what Michael Trott's code does, we have to write something like the following which shows both the real and the imaginary surfaces.

α = 0.3 + 0.1 I;
ee = 10^(-12);
indxs = Range[-2, 2];
colors = Most[Hue /@ Subdivide[Length[indxs]]];
plotF =
  Function[u,
    MapThread[
      ParametricPlot3D[
        {r Cos[φ],
         r Sin[φ],
         u[Exp[2 #1 Pi I α] Gamma[α, r Exp[I φ]] + (1 - Exp[2 #1 Pi I α]) Gamma[α]]}, 
        {r, 0, 2}, {φ, -Pi + ee, Pi - ee},
        PlotPoints -> {30, 40},
        BoxRatios -> {1, 1, 1},
        PlotStyle -> #2] &, 
      {indxs, colors}]];
Column[
  Legended[
    Show[plotF[#], PlotRange -> All, ImageSize -> 400],
    SwatchLegend[colors, indxs, LegendLabel -> Style["k", "TBI", 16]]] & /@ {Re, Im}]

Answer 3

For reference, this is what the code produces in version 5.2 (the latest version at the time the bulk of Trott's book was written):

Now, the reason for the Cases[(* stuff *), _Polygon, ∞] and the Last /@ Partition[(* stuff *)] business is that it was Trott's method to help reveal the structure of the Riemann surface, by extracting all polygons from the plots, and then deleting every other polygon in the plotted surfaces.

In the current version of Mathematica, a little more work is necessary to reproduce this figure. Instead of manually deleting polygons, we can instead fiddle with the Mesh and MeshShading settings, like so:

With[{α = 0.3 + 0.1 I, ε = 1*^-12},
     GraphicsRow[Table[
                 ParametricPlot3D[Table[{r Cos[φ], r Sin[φ],
                                         ff[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                            (1 - Exp[2 k π I α]) Gamma[α]]},
                                        {k, -2, 2}],
                                  {r, 0, 2}, {φ, -π + ε, π - ε},
                                  BoxRatios -> {1, 1, 2.5}, Lighting -> "Classic", 
                                  Mesh -> {31, 41},
                                  MeshShading -> {{None, Automatic}, {Automatic, None}}, 
                                  MeshStyle -> Thickness[0.002], PlotRange -> All,
                                  PlotStyle -> Directive[Hue[0.09],
                                                         Specularity[Hue[0.18], 2.3]]],
                            {ff, {Re, Im}}]]]

---

Nowadays, however, we have at our disposal a graphics directive that Trott didn't: Opacity[]. We can use it to have a good look at the Riemann surface's structure by making the surface translucent, and forgo all the polygon removal business:

With[{α = 0.3 + 0.1 I, ε = 1*^-12},
     GraphicsRow[Table[ParametricPlot3D[Table[{r Cos[φ], r Sin[φ],
                                               ff[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                                               (1 - Exp[2 k π I α]) Gamma[α]]},
                                              {k, -2, 2}],
                                        {r, 0, 2}, {φ, -π + ε, π - ε},
                                        BoxRatios -> {1, 1, 2.5}, Mesh -> None,
                                        PlotRange -> All,
                                        PlotStyle -> Directive[Hue[0.09], Opacity[0.6]]],
                       {ff, {Re, Im}}]]]

Answer 4

Here is code to plot two sheet separately. You can add to it:

k = 0; α = 0.3 + 0.1 I; ε = 10^(-12)
p1 = ParametricPlot3D[{r Cos[φ], r Sin[φ], 
                       Re[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                          (1 - Exp[2 k π I α]) Gamma[α]]},
                      {r, 0, 2}, {φ, -π + ε, π - ε}, PlotPoints -> {30, 40}, 
                      BoxRatios -> {1, 1, 1}, PlotStyle -> Red];
k = -1;
p2 = ParametricPlot3D[{r Cos[φ], r Sin[φ], 
                       Re[Exp[2 k π I α] Gamma[α, r Exp[I φ]] +
                          (1 - Exp[2 k π I α]) Gamma[α]]},
                      {r, 0, 2}, {φ, -π + ε, π - ε}, PlotPoints -> {30, 40}, 
                      BoxRatios -> {1, 1, 1}, PlotStyle -> Blue];
Show[{p1, p2}, PlotRange -> All]