Let's examine the code and describe what it does.
To compute values of a modular form, we typically work with its $q$-expansion, which looks like
$$ f(z) = \sum_{n \geq 0} a(n) q^n. \qquad \qquad (q = e^{2 \pi i z}). \tag{1} $$
This is a Fourier expansion, but for old reasons that I don't know it is very common to write $q$. The coefficients are the Fourier coefficients [1]. In the code snippet, we get the $q$-expansion of a modular form in the line
f = ModularForms(group=1, weight=12).newforms()[0].q_expansion(100).truncate()
In more detail: ModularForms(group, weight) returns the modular form space for a given congruence subgroup and given weight. Here the congruence subgroup is $\Gamma(1) = \mathrm{SL}(2, \mathbb{Z})$, and the weight is $12$. The method .newforms() returns the newforms in the cuspidal subspace. These are forms that don't come from simpler modular forms on simpler groups. Here this is a bit silly, as this is the "simplest group" in a sense and so all cusp forms are newforms. Further, there is exactly one Eisenstein series and exactly one cuspform.
The next portion, [0].q_expansion(100) chooses the unique cuspidal newform in this space and computes its first 100 Fourier coefficients as a power series in $q$. This gives computes that
$$ f(z) = q - 24q^2 + \cdots - 60754911516q^{99} + O(q^{100}). $$
The remaining part, .truncate(), is a method to make sage drop the $O(q^{100})$ part and consider the computed polynomial approximation to the power series. It truncates the remaining terms at the indicated precision.
This is one way to generate the approximating $q$-polynomial
$$ f(z) \approx q - 24q^2 + \cdots - 60754911516q^{99}, \tag{2} $$
which we'll use to compute the values of the modular form $f$.[2]
Htoq = lambda x: exp(2*CDF.pi()*CDF.0*x)
DtoH = lambda x: (-CDF.0*x + 1)/(x - CDF.0)
Dtoq = lambda x: Htoq(DtoH(CDF(x)))
Now that we have a (polynomial approximation to the) $q$-expansion, we can compute (an approximation to) $f(q)$ at any point $q$. But to plot the value of $f(z)$ at a point in the upper half-plane, we need to compute $f(z)$. These letters represent different choices of coordinates. The three lines of code above are coordinate transformations.[3]
Htoq means "given a point $z$ in the the upper halfplane $\mathcal{H}$, return the corresponding point $q$ for the $q$-expansion". As in $(1)$, we know how to write $q$ in terms of $z$: we use $q = e^{2 \pi i z}$. This is exactly what Htoq does.
In more detail: the code uses 53 bit precision complex fields for speed. In sage, this is CDF (standing for "complex double field"). The two constants CDF.pi() and CDF.0 stand for $\pi$ and $i$, represented in double precision. Thus Htoq computes precisely $q = e^{2 \pi i z}$.
You can make plots on portions of the upper halfplane now, using something like
complex_plot(
lambda z: f(Htoq(z)), # function to plot
(-1, 1), # x range
(0.01, 2) # y range
)
The way to think of this is as follows: for each $x$ value in $(-1, 1)$, and for each $y$ value in $(0.01, 2)$, write $z = x + iy$ and compute $f($Htoq(z)$)$. Then for each $z$ value, Htoq(z) gives the corresponding $q$ value, and this is the expression we know how to compute $f$ via $(2)$.
The second coordinate transformation, DtoH, represents changing from a model of the disk to the upper halfplane. This is the Cayley transformation
$$ \frac{-iz + 1}{z - i}, $$
mapping the unit disk to $\mathcal{H}$. It sends the point $i$ (the top of the unit disk) to $\infty$, $0$ (the center of the disk) to $i$, and $-i$ (the bottom of the unit disk) to $0$. Noting that this is a Mobius transform and thus sends lines+circles to lines+circles, this shows that this transformation preserves the apparent imaginary orientation along the line $\mathrm{Re } z = 0$. That's how I chose this transform.
The last coordinate transform $Dtoq$ takes a coordinate on the disk and outputs its $q$ expansion. It does this via composition: first map the coordinate to $\mathcal{H}$ (via the Cayley transformation), and then write that $z$ as $q$ (via $q = e^{2 \pi i z}$).
Finally, we look at the last line. I write this on several different lines and discuss each.
P = ccomplex_plot(
lambda x: +Infinity if abs(x) >= 0.99 else f(Dtoq(x)),
(-1,1),
(-1,1),
plot_points=500, aspect_ratio = 1, figsize=[5,5]
)
The first bit, ccomplex_plot, uses beta code before I added my complex routines to sage. In Sage9.6+, one can just use complex_plot.
The next three lines indicate that the plot should range over all $z$ in the square $[-1, 1] \times [-1, 1]$. The function calculation is a bit of an oddity. It turns out that sage's complex_plot colors NaN values or inf values as white. So one way to plot only the disk and not do anything with masking or something more clever is to explicitly define your function to be infinite outside the disk — in plots it will be white, like the default background. Inside the disk, for each $z$ in the disk, $f($Dtoq(z)$)$ converts $z$ from disk coordinates to $q$-coordinates as above, and plots it.
The remaining bits are sage plotting options. They set the size of the plot and the aspect ratio. The argument plot_points is important. One cannot compute every value of $f$ on the disk, as there are infinitely many. Instead, one computes values on a mesh of points and then interpolates between them for colors, and plot_points controls how fine the mesh is. More plot points means more points are actually computed.
Additional notes and references.
[1]: These are the Fourier coefficients at the cusp at $\infty$. There are possibly other cusps, with possibly different $q$-expansions. This can be confusing. For more about this, I suggest reading the first chapter of Diamond and Shurman — Iwaniec's Topics in Classical Modular Forms for a more rigorous, less friendly description.
[2]: There are other ways. This particular modular form is particularly famous, and you can get the exact same representation with delta_qexp(100).truncate(). These are also already computed and stored as Modular form 1.12.a.a in the LMFDB, and it's possible to retrieve these expansions directly from there.
[3]: This is dense code. Using lambdas like this to give one-line function definitions is a code smell. Explicitly giving names to anonymous functions is a code smell. But I think it expresses the mathematics very cleanly.
