Note: In this question I am concerned only with real-valued variables and functions.
DLMF, §9.8 Airy Functions, Modulus and Phase, formula $9.8.4$ defines the phase of Airy functions: $$\theta(x)=\arctan\frac{\operatorname{Ai}x}{\operatorname{Bi}x}.$$ The corresponding Mathematica definition is
θ[x_] := ArcTan[AiryAi[x]/AiryBi[x]]
Apparently, in this form the phase has a countable infinite number of jump discontinuities for negative values of $x$.
I need to construct a Mathematica expression that represents a continuous monotonic version of the phase $\theta(x)$ (let's name it $\vartheta$), where all pieces between discontinuities are properly shifted and stitched, as shown on the following graph:

One approach would be to get the derivative of $\theta(x)$: $$\theta'(x)=-\frac1\pi\frac1{\operatorname{Ai}^2 x+\operatorname{Bi}^2 x},$$ so that all information about jumps is lost, and then consider the definite integral: $$\vartheta(x)=\frac1\pi\int_x^\infty\frac{dz}{\operatorname{Ai}^2 z+\operatorname{Bi}^2 z}.$$ Or, in Mathematica notation:
1/π Integrate[1/(AiryAi[z]^2 + AiryBi[z]^2), {z, x, ∞}]
Unfortunately, Mathematica leaves this integral unevaluated, even if an explicit value for the boundary $x$ is provided.
Is it possible to write a closed-form Mathematica expression, possibly containing special functions, that represents the continuous monotonic phase $\vartheta(x)$?


