(Image: Nick84, Wikimedia Commons.)

— Notes by Steven G. Johnson, December 2019. If you find this work useful, please cite: *Steven G. Johnson, "Accurate solar-power integration: Solar-weighted Gaussian quadrature," arXiv:1912.06870 (2019).*

For many problems involving solar energy, e.g. simulating the efficiency of solar cells, it is necessary to compute integrals (also called *quadrature* in numerical contexts) over wavelength λ of the form:

for some desired wavelength interval $(a,b)$, where $f(\lambda)$ is application-specific and $S(\lambda)$ is a solar irradiance spectrum: the **intensity of sunlight measured on Earth**. For example, to evaluate the efficiency of solar cells, typically $f(\lambda) = \lambda A(\lambda)$ where $A(\lambda)$ is the absorption spectrum of the photovoltaic device.

As is shown by the red curve above, S(λ) is an extremely complicated function tabulated experimentally by the American Society for Testing and Materials (ASTM), displaying many sharp dips and oscillations due to atmosopheric absorption by different molecules.

The challenge is that the application-specific f(λ) function is often **extremely expensive** to evaluate, e.g. for solar cells it involves solving Maxwell's equations (solving a partial differential equation) for each λ.

Even if f(λ) is a smooth function, after multiplying by the solar weight S(λ) one obtains a very badly behaved integrand, requiring a **huge number of function evaluations** to obtain good accuracy.

We can solve this problem *very* efficiently because the **weight function** S(λ) is **fixed**, allowing us to **precompute** a quadrature scheme that takes the "nastiness" of S(λ) into account *analytically*, yielding a formula of the form:

where the $N$ points $x_i$ are known as *quadrature points* and the coefficients $w_i$ are the corresponding *quadrature weights*. The goal is to find these points and weights so that the integral is accurately approximated for smooth functions f(λ), where the effect of S(λ) is *precomputed* in $x_i$ and $w_i$.

A well-known method to accomplish this goal is Gaussian quadrature, in which clever theoretical methods are used to construct an $N$-point quadrature scheme that integrates *polynomial* functions $f(\lambda)$ *exactly* up to degree $2N-1$. Moreover, Gaussian quadrature converges extremely quickly with increasing $N$ for *smooth* functions f(λ). All of this is true *regardless* of the smoothness (or lack thereof) of S(λ), since S(λ) is taken into account in the construction of the quadrature scheme.

The construction of Gaussian quadrature schemes is reviewed extensively elsewhere. In brief, however, we first define an inner product $\langle p_1, p_2 \rangle = \int_a^b S(x) p_1(x) p_2(x) dx$ of real polynomials on $(a,b)$. Then we perform a Gram–Schmidt orthogonalization of polynomials $\{1,x,x^2,\ldots\}$ to find orthogonal polynomials $\{q_0,q_1,\ldots\}$ with *respect to this weighted inner product*. The Gaussian quadrature points $x_i$, amazingly, are the roots of $q_N(x)$. Moreover, as is reviewed in e.g. chapter 37 of Trefethen and Bau, these roots (and also the quadrature weights) can be obtained by solving a tridiagonal $N\times N$ eigenvalue problem whose entries are obtained by a Lanczos algorithm. Essentially, we just need to integrate $2N$ polynomials against S(λ) — since these integrals only need to be done *once*, and are of cheap functions (polynomials), we can use expensive brute-force methods to perform them.

All of this process, for an *arbitrary* weight function, is implemented for us in the Julia QuadGK.jl package. I should also mention one other wrinkle: for numerical stability when constructing Gaussian quadrature schemes computationally for arbitrary weight functions, QuadGK applies the Lanczos process to polynomials in the basis $\{T_0,T_1,T_2,\ldots\}$ of Chebyshev polynomials, after rescaling the interval $(a,b)$ to $(-1,1)$, rather than using "textbook" basis $\{1,x,x^2,\ldots\}$.

Below just need to supply QuadGK the solar-irradiance S(λ), and give it a bit of information on how to do the brute-force integrals on interpolated data by breaking it up into intervals. Then it will construct our solar-weighted Gaussian quadrature scheme to any desired order $N$ and for any wavelength interval!

The following Julia packages must be installed for this notebook to work: QuadGK.jl (numerical integration), Dierckx.jl (cubic splines), PyPlot.jl (plots), and ZipFile.jl (reading zip-compressed data).

Uncomment and execute the following line to install them. (In IJulia, cells starting with `]`

are commands sent to the Julia package manager.)

In [1]:

```
# ] add QuadGK Dierckx PyPlot ZipFile
```

To begin, let us download and import the ASTM solar irradiance data. We will use the ZipFile.jl package to unpack the `.zip`

file containing the comma-delimited data, and then the `DelimitedFiles`

standard-library package to read the data.

In [2]:

```
using ZipFile, DelimitedFiles
ASTM = let filename = download("https://www.nrel.gov/grid/solar-resource/assets/data/astmg173.zip")
r = first(ZipFile.Reader(filename).files)
data, header = readdlm(IOBuffer(read!(r, Array{UInt8}(undef, r.uncompressedsize))), ',', header=true, skipstart=1)
data
end
λ, S = ASTM[:,1], ASTM[:,4]
```

Out[2]:

The data consists of 2002 data points for wavelengths (λ) from 280nm to 4000nm:

In [3]:

```
@show extrema(λ)
@show length(λ);
```

Let's plot the data to see the "messy" irradiance spectrum:

In [4]:

```
using PyPlot
figure(figsize=(10,4))
subplot(1,2,1)
plot(λ, S)
xlabel("wavelength (nm)")
ylabel("solar irradiance (W/sm/nm)")
title("Terrestrial solar-irradiance spectrum")
subplot(1,2,2)
semilogy(λ, S)
xlabel("wavelength (nm)")
title("(log scale)")
ylim(1e-5,10)
```

Out[4]:

Since we are planning on integrating this data, it is nice to see how far out in the decaying "tail" we need to go. To do this, let's compute the fraction of the integral:

$$ \mathrm{fraction}(\lambda) = \frac{\int_{280}^\lambda S(\lambda') d\lambda'}{\int_{280}^{4000} S(\lambda') d\lambda'} $$as a function of λ, superimposed on a normalized spectrum.

(For this plot, we will estimate the integral using a simple Euler approximation.)

In [5]:

```
title("Irradiance data and integral fraction")
dλ = diff(λ)
semilogy(λ[1:end-1], 1 .- cumsum(S[1:end-1] .* dλ) / sum(S[1:end-1] .* dλ), "-")
semilogy(λ, S / maximum(S), "-")
xlabel("wavelength (nm)")
legend(["fraction of integral", "irradiance, normalized"])
ylim(1e-4, 1)
```

Out[5]:

From the above plot, we can see that a significant fraction (≈ 1%) of the integral is in the "tail" from 3000–4000nm, so we need the whole spectrum if we want to compute integrals to many decimal places (unless the application's function f(λ) is rapidly decaying too).

To compute integrals, we need more than just a set of data points — we need to **interpolate** this tabulated data so that we can evaluate S(λ) for any λ in (200nm,4000nm).

Since the underlying physical spectrum is a smooth function, unless the data are extremely noisy it makes sense to use a cubic spline) interpolation, which interpolates the data by a sequence of cubic polynomials so that the first and second derivatives are also continuous.

However, there is a wrinkle — for Gaussian quadrature, we need to make sure that $S(\lambda) \ge 0$ everywhere (since $S$ must produce a weighted inner product on polynomials). Our data points are all $> 0$, but it is conceivable that a spline interpolant will dip below zero at intermediate points. To ensure that this is not the case, we will fit our cubic spline to $\sqrt{S}$ at the data points, and then square the result.

We will do this using the Dierckx package in Julia, a wrapper around a popular Fortran library that supports irregularly-spaced data like the ASTM solar data.

In [6]:

```
using Dierckx
let Sinterp_sqrt = Spline1D(λ, sqrt.(S)) # spline for √S
global Sinterp
Sinterp(x) = Sinterp_sqrt(x)^2
end
```

Out[6]:

The "knots" of a spline are the points where the interpolation switches from one polynomial to the next, i.e. points at which some higher derivative is discontinuous. We will come back to this below, since knowing the points of discontinuity is crucial to computing our weight-function integrals quickly and accurately.

Let's plot the spline along with the underlying data in a smaller wavelength interval 400–500nm so that we can see the wiggles:

In [7]:

```
plot(λ, S, "b.")
xlabel("wavelength (nm)")
ylabel("solar irradiance (W/sm/nm)")
title("Cubic-spline interpolated spectrum")
xlim(400,500)
λinterp = range(400,500,length=1000)
plot(λinterp, Sinterp.(λinterp), "r-")
legend(["data", "spline"])
```

Out[7]:

The spline fit looks reasonable — notice that the rapid oscillations encompass several data points in a row, so they are presumably real features of the data and not simply noise.

The QuadGK package does almost all of the work for us in constructing a Gaussian quadrature rule for a given weight function, its `x, w = gauss(...)`

functions that returns arrays of quadrature points `x`

and weights `w`

.

All we need to do is to give it the weight function, in this case our spline-interpolated data. `gauss`

will then laboriously compute a sequence of integrals of this weight function against polynomials.

Although we don't care *too* much about the speed of this process, since we only need to do it once for a given quadrature rule, it turns out that there is something simple we can do to speed it up. Numerical integration is best for *smooth* functions, but our spline interpolant is *not* smooth: it has discontinuities (in some higher derivative) at the "knot" points (≈ data points). To integrate such a *piecewise-smooth* function accurately, we should simply perform the integral piecewise. The `gauss`

function allows us to pass a `quad`

argument with a custom integration routine, so below we pass a specialized `interpquad`

routine that performs integrals by breaking them up into pieces at the knots. (Remember: we are doing this work only during the *construction* of the quadrature rule for S(λ), not during the evaluation of the resulting integrals of f(λ).)

The result, below, is a `x, w = gaussquad_interpolant(N, X, W, a, b)`

function that returns the `N`

-point quadrature rule on the interval `(a, b)`

for a weight function interpolated from data points `(X[j],W[j])`

.

In [8]:

```
using QuadGK, Dierckx
function gaussquad_interpolant(N::Integer,
X::AbstractVector{<:Real}, W::AbstractVector{<:Real},
a::Real=minimum(X), b::Real=maximum(X);
rtol::Real=sqrt(eps(typeof(float(b-a)))),
interpdegree::Integer=3)
# some quick argument sanity checks:
all(w ≥ 0 for w in W) || throw(ArgumentError("weights must be nonnegative"))
length(X) == length(W) || throw(DimensionMismatch("points and weights must have same length"))
N > 1 || throw(ArgumentError("a positive order N is required"))
a < b || throw(ArgumentError("endpoints b > a are required"))
rtol > 0 || throw(ArgumentError("a positive rtol is required"))
interpdegree > 0 || throw(ArgumentError("a positive interpdegree is required"))
# construct a cubic-spline interpolation from the data.
# … we fit √W and then square it below to ensure non-negativity.
Winterp_sqrt = Spline1D(X, sqrt.(W), k=interpdegree)
# break integration interval at knots
xab = sort!(collect(Iterators.filter(x -> a < x < b, Dierckx.get_knots(Winterp_sqrt))))
push!(pushfirst!(xab, a), b) # add endpoints
# quadrature routine for W-weighted integrands, that breaks integrand up into each
# interpolation interval at the knots.
interpquad(f, _a, _b; kws...) = quadgk(f, xab...; kws...)
return gauss(x -> Winterp_sqrt(x)^2, N, a, b; quad=interpquad, rtol=rtol)
end
```

Out[8]:

Now that we have the short function above, let's apply it to compute a 15-point Gaussian quadrature rule for the S(λ) weight function. This will *exactly* integrate polynomial f(λ) up to degree 29, and should be good for any smooth functions as long as they aren't oscillating much faster than our 15 points can capture.

Here we calculate it, time how long the calculation takes with Julia's `@time`

macro, and show the quadrature points and weights:

In [9]:

```
@time x15, w15 = gaussquad_interpolant(15, λ, S)
println()
@show x15
println()
@show w15;
```

The calculation took about 22 seconds, which is not too bad since we only have to do this once. But actually the calculation is faster than that — most of these 22 seconds are spent by Julia in JIT compiling the computational routines.

Subsequently running the code, now that it has been compiled, is much quicker:

In [10]:

```
@time gaussquad_interpolant(15, λ, S);
@time gaussquad_interpolant(15, λ, S);
@time gaussquad_interpolant(15, λ, S);
```

We can also plot the points and weights. Unsurprisingly, they vaguely resemble the solar spectrum itself, with more weight where the solar irradiance is greater:

In [11]:

```
semilogy(x15, w15, "bo-")
xlabel(L"quadrature point $\lambda$ (nm)")
ylabel("quadrature weight")
title("15-point solar-weighted Gaussian quadrature")
```

Out[11]:

We'll also compute a **99-point** solar-weighted quadrature rule, suitable for integrating much more oscillatory functions.

In [12]:

```
@time x99, w99 = gaussquad_interpolant(99, λ, S)
println()
@show x99
println()
@show w99
semilogy(x99, w99, "go-")
xlabel(L"quadrature point $\lambda$ (nm)")
ylabel("quadrature weight")
title("99-point solar-weighted Gaussian quadrature");
```

Notice that it puts few quadrature points around wavelengths like 2500–3000nm where the solar irradiance is almost zero (due to strong absorption peaks). It still puts a lot of data points in the 3000–4000nm range where the irradiance is small but not tiny (≈ 1% of total), because integrating these ranges accurately is necessary if we want numeric integrals that are correct to many decimal places.

Computing the 99-point rule was more expensive, but still took less than 6 seconds. If we need a 999-point rule (to integrate *extremely* oscillatory f(λ) functions), we could do that too without too much trouble.

To illustrate these quadrature rules and see how they compare to simpler methods, let's integrate the test function: $$ f(\lambda) = \sin(2\pi \lambda / 500) $$ which oscillates with a period of 500nm (eight times in the ASTM data's 280–4000nm range).

(This is not too wiggly, but integrating it against S(λ) accurately is very challenging for simpler integration techniques.)

Let's plot it along with the solar spectrum:

In [13]:

```
f(x) = sin(2π*x/500) # test integrand
```

Out[13]:

In [14]:

```
plot(λ, f.(λ), "r-")
plot(λ, Sinterp.(λ), "k--")
plot(x15, f.(x15), "bo")
legend(["test integrand", "solar spectrum", "15-point quadrature"])
xlabel(L"\lambda")
```

Out[14]:

Above, the blue dots show the function f(x) sampled at the quadrature points $x_i$ of our 15-point solar-weighted Gaussian quadrature scheme. This is **not very many points**! Even without the S(λ) weight, extremely simple methods like the trapezoidal rule would not be able to integrate this f(x) accurately with so few points. Capturing the S(λ) accurately as well would be crazy, except that the effect of S(λ) is already *exactly* taken into account by the solar-weighted quadrature scheme.

To begin with, we'll just compute the "exact" integral (to nearly machine precision) by the brute-force Gauss–Kronrod quadrature function `quadgk`

provided by the QuadGK package.

We won't be *completely* naive, however: instead of passing `280, 4000`

as the integration limits, we'll pass `λ...`

, which "splats" all of the data points: this causes `quadgk`

to break the integral into 2001 sub-segments where our interpolant `Sinterp`

is smooth, giving an extremely accurate result at the cost of lots of function evaluations. We'll increment a global counter `cnt`

in our integrand routine to keep track of the number of times `quadgk`

evaluates our integrand.

The resulting integral `81.264730...`

is accurate to nearly 15 significant digits, at the cost of 30,000 function evaluations. That's fine for such a cheap `f(λ)`

, but not if evaluating `f(λ)`

required us to solve Maxwell's equations once for each λ!

In [15]:

```
cnt = 0
I,E = quadgk(x -> begin; global cnt += 1; Sinterp(x) * f(x); end, λ...)
println("got ∫S(λ)f(λ)dλ = $I ± $E after $cnt evaluations")
```

Next, we evaluate the same integral using our 15-point and 99-point solar-weighted Gaussian quadrature rules from above, and compare the results to the brute-force "exact" answer above.

Using **only 15 points**, we now obtain better than **1%** accuracy! Using only 99 points, our result is accurate to **13 digits**.

Without knowing about the magic of Gaussian quadrature rules with precomputed weights, it might seem unimaginable that one could compute this crazy S(λ)f(λ) integral so accurately with so few points.

In [16]:

```
I15 = sum(f.(x15) .* w15)
println("got ∫S(λ)f(λ)dλ = $I15 ± $(abs(I15-I)) after 15 evaluations")
```

In [17]:

```
I99 = sum(f.(x99) .* w99)
println("got ∫S(λ)f(λ)dλ = $I99 ± $(abs(I99-I)) after 99 evaluations")
```

If you are integrating against tabulated data like we have here, the first technique that might come to mind is the trapezoidal rule. If we have M data points for S(λ), we simply evaluate our integrand S(λ)f(λ) at **all M points** and interpolate with straight lines in between (adding up the area of trapezoids).

Here, that scheme is actually fairly accurate, to about 4 significant digits, but requires $M=2002$ function evaluations so it is horribly impractical for expensive integrands.

In [18]:

```
dλ = diff(λ)
fS = Sinterp.(λ) .* f.(λ)
Itrap = 0.5 * (sum(fS[1:end-1] .* dλ) + sum(fS[2:end] .* dλ))
println("got ∫S(λ)f(λ)dλ = $Itrap ± $(abs(I-Itrap)) after $(length(λ)) evaluations")
```

Another naive approach might be to simply use an adaptive quadrature method that evaluates the integrand at more and more points until a desired accuracy is achieved. In fact, the `quadgk`

routine provided by the QuadGK package is one such algorithm, based on a Gauss–Kronrod method.

However, to resolve the extremely rapid oscillations of S(λ), especially since the **interpolated** data is **not smooth** (the third derivative is discontinuous), adaptive quadrature again takes a huge number of function evaluations.

Here, we see that it takes 5295 evaluations for 1% accuracy, and 12135 evaluations for 0.003% accuracy. Again, this is not practical for expensive integrands.

In [19]:

```
cnt = 0
I′,E′ = quadgk(x -> begin; global cnt += 1; Sinterp(x) * f(x); end, minimum(λ), maximum(λ), rtol=1e-2)
println("got ∫S(λ)f(λ)dλ = $I′ ± $(abs(I-I′)) after $cnt evaluations")
```

In [20]:

```
cnt = 0
I′,E′ = quadgk(x -> begin; global cnt += 1; Sinterp(x) * f(x); end, minimum(λ), maximum(λ), rtol=1e-3)
println("got ∫S(λ)f(λ)dλ = $I′ ± $(abs(I-I′)) after $cnt evaluations")
```

As mentioned above, the accuracy of Gaussian quadrature typically improves *exponentially* ("spectral" convergence = faster than power-law) with the number $N$ of quadrature points if $f(\lambda)$ is a smooth function. Let's see that here by plotting the error vs. $N$, both for the $f(\lambda)$ above and for a second function $f_2(\lambda)$ that oscillates $10\times$ more rapidly.

In [21]:

```
f₂(x) = sin(2π*x/50) # oscillates 10× faster than f(λ) from above
# brute-force integral, for refernce:
cnt = 0
I₂,E₂ = quadgk(x -> begin; global cnt += 1; Sinterp(x) * f₂(x); end, λ...)
println("got ∫S(λ)f₂(λ)dλ = $I₂ ± $E₂ after $cnt evaluations")
```

In [22]:

```
Ns = [5:2:33; 40:10:140]
err = Float64[]
err₂ = Float64[]
rule = Dict() # save quadrature rules in case we want them later
print("Working on N: ")
@time begin
for N in Ns
print(N, ", ")
flush(stdout) # force immediate progress output
x, w = gaussquad_interpolant(N, λ, S)
push!(err, abs(sum(f.(x) .* w) - I) / abs(I))
push!(err₂, abs(sum(f₂.(x) .* w) - I₂) / abs(I₂))
rule[N] = (x,w)
end
println("\nDone.\n")
end
```

In [23]:

```
semilogy(Ns, err, "ro-")
semilogy(Ns, err₂, "b*-")
xlabel(L"quadrature order $N$")
ylabel(L"fractional error in $\int S(\lambda) f(\lambda) d\lambda$")
legend([L"f(\lambda) = \sin(2\pi\lambda/500)", L"f_2(\lambda) = \sin(2\pi\lambda/50)"])
title("convergence of solar-weighted Gaussian quadrature")
grid()
```

Notice that when the integrand $f(\lambda)$ is rapidly oscillating (or has other "sharp" features), you still need to sample it finely enough to capture the features of $f(\lambda)$. Only after you use a large enough $N$ [to roughly fit $f(\lambda)$ to a polynomial of degree $N$] does the spectral convergence "kick in." (Still, we get 9 digits with just 140 points.)

For this $f_2$, even a 99-point quadrature rule is not enough: sampling $f_2(\lambda)$ at just 99 points clearly leads to aliasing effects in which the function appears to be oscillating at a much lower frequency:

In [24]:

```
λ₂ = range(minimum(λ), maximum(λ), length=10^4) # need a finer grid to plot f₂
plot(λ₂, f₂.(λ₂), "r-")
plot(λ₂, Sinterp.(λ₂), "k--")
plot(x99, f.(x99), "bo--")
legend([L"test integrand $f_2(\lambda)$", "solar spectrum", "99-point quadrature"])
xlabel(L"\lambda")
title(L"aliasing of rapidly oscillating $f_2(\lambda)$")
```

Out[24]:

To evaluate the efficiency of solar cells, typically $f(\lambda) = \lambda A(\lambda)$ where $A(\lambda)$ is the absorption spectrum of the photovoltaic device. However, for silicon solar cells the absorption decays extremely (exponentially) rapidly for longer wavelengths, so there is not much point in integrating beyond 1100nm in many applications.

We might therefore want quadrature rules for **integrating from 280–1100nm**, and this is easy to do with our routines above. Here are 15-point and 99-point quadrature schemes, for example:

In [25]:

```
@time x15_280_1100, w15_280_1100 = gaussquad_interpolant(15, λ, S, 280, 1100)
println()
@show x15_280_1100
println()
@show w15_280_1100
semilogy(x15_280_1100, w15_280_1100, "go-")
xlabel(L"quadrature point $\lambda$ (nm)")
ylabel("quadrature weight")
title("15-point solar-weighted 280–1100nm quadrature");
```

In [26]:

```
@time x99_280_1100, w99_280_1100 = gaussquad_interpolant(99, λ, S, 280, 1100)
println()
@show x99_280_1100
println()
@show w99_280_1100
semilogy(x99_280_1100, w99_280_1100, "go-")
xlabel(L"quadrature point $\lambda$ (nm)")
ylabel("quadrature weight")
title("99-point solar-weighted 280–1100nm quadrature");
```