How to use the starry._core.utils.autocompile function in starry

    @autocompile
    def RAxisAngle(self, axis=[0, 1, 0], theta=0):
        """Wigner axis-angle rotation matrix."""

        def compute(axis=[0, 1, 0], theta=0):
            axis = tt.as_tensor_variable(axis)
            axis /= axis.norm(2)
            cost = tt.cos(theta)
            sint = tt.sin(theta)

            return tt.reshape(
                tt.as_tensor_variable(
                    [
                        cost + axis[0] * axis[0] * (1 - cost),
                        axis[0] * axis[1] * (1 - cost) - axis[2] * sint,
                        axis[0] * axis[2] * (1 - cost) + axis[1] * sint,
                        axis[1] * axis[0] * (1 - cost) + axis[2] * sint,

    @autocompile
    def X_point_source(
        self, theta, xs, ys, zs, xo, yo, zo, ro, inc, obl, u, f, alpha, sigr
    ):
        """Compute the light curve design matrix for a point source."""
        # Determine shapes
        rows = theta.shape[0]
        cols = (self.ydeg + 1) ** 2
        X = tt.zeros((rows, cols))

        # Compute the occultation mask
        bo = tt.sqrt(xo ** 2 + yo ** 2)
        b_rot = tt.ge(bo, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
        b_occ = tt.invert(b_rot)
        i_rot = tt.arange(bo.size)[b_rot]
        i_occ = tt.arange(bo.size)[b_occ]

    @autocompile
    def compute_rect_grid(self, res):
        """Compute the polynomial basis on a rectangular lat/lon grid."""
        # See NOTE on tt.mgrid bug in `compute_ortho_grid`
        dx = np.pi / (res - 0.01)
        lat, lon = tt.mgrid[
            -np.pi / 2 : np.pi / 2 : dx, -3 * np.pi / 2 : np.pi / 2 : 2 * dx
        ]
        x = tt.reshape(tt.cos(lat) * tt.cos(lon), [1, -1])
        y = tt.reshape(tt.cos(lat) * tt.sin(lon), [1, -1])
        z = tt.reshape(tt.sin(lat), [1, -1])
        R = self.RAxisAngle(tt.as_tensor_variable([1.0, 0.0, 0.0]), -np.pi / 2)
        return (
            tt.concatenate(
                (tt.reshape(lat, [1, -1]), tt.reshape(lon, [1, -1]))
            ),
            tt.dot(R, tt.concatenate((x, y, z))),

    @autocompile
    def X(
        self,
        t,
        pri_r,
        pri_m,
        pri_prot,
        pri_t0,
        pri_theta0,
        pri_amp,
        pri_inc,
        pri_obl,
        pri_u,
        pri_f,
        pri_alpha,
        sec_r,
        sec_m,

    @autocompile
    def rv(self, theta, xo, yo, zo, ro, inc, obl, y, u, veq, alpha):
        """Compute the observed radial velocity anomaly."""
        # Compute the velocity-weighted intensity
        f = self.compute_rv_filter(inc, obl, veq, alpha)
        Iv = self.flux(theta, xo, yo, zo, ro, inc, obl, y, u, f, alpha)

        # Compute the inverse of the intensity
        f0 = tt.zeros_like(f)
        f0 = tt.set_subtensor(f0[0], np.pi)
        I = self.flux(theta, xo, yo, zo, ro, inc, obl, y, u, f0, alpha)
        invI = tt.ones((1,)) / I
        invI = tt.where(tt.isinf(invI), 0.0, invI)

        # The RV signal is just the product
        return Iv * invI

    @autocompile
    def compute_illumination_point_source(
        self, xyz, xs, ys, zs, sigr, on94_exact
    ):
        """Compute the illumination profile for a point source."""
        # Get cos(theta_i)
        x = tt.shape_padright(xyz[0])
        y = tt.shape_padright(xyz[1])
        z = tt.shape_padright(xyz[2])
        r2 = xs ** 2 + ys ** 2 + zs ** 2
        b = -zs / tt.sqrt(r2)  # semi-minor axis of terminator
        invsr = 1.0 / tt.sqrt(xs ** 2 + ys ** 2)
        cosw = ys * invsr
        sinw = -xs * invsr
        xrot = x * cosw + y * sinw
        yrot = -x * sinw + y * cosw
        bc = tt.sqrt(1.0 - b ** 2)

    @autocompile
    def spotYlm(self, amp, sigma, lat, lon):
        return self._spotYlm(amp, sigma, lat, lon)

    @autocompile
    def dotR(self, matrix, ux, uy, uz, theta):
        return self._dotR(matrix, ux, uy, uz, theta)

    @autocompile
    def sT(self, b, r):
        return self._sT(b, r)

    @autocompile
    def tensordotD(self, matrix, wta):
        return self._tensordotD(matrix, wta)