How to use the floris.utilities.cosd function in floris

for j in range(len(yt)):
                for k in range(len(zt)):
                    x_grid[i, j, k] = xt[i]
                    y_grid[i, j, k] = yt[j]
                    z_grid[i, j, k] = zt[k]

                    xoffset = x_grid[i, j, k] - coord.x1
                    yoffset = y_grid[i, j, k] - coord.x2
                    x_grid[i, j, k] = (
                        xoffset * cosd(-1 * self.wind_map.turbine_wind_direction[i])
                        - yoffset * sind(-1 * self.wind_map.turbine_wind_direction[i])
                        + coord.x1
                    )

                    y_grid[i, j, k] = (
                        yoffset * cosd(-1 * self.wind_map.turbine_wind_direction[i])
                        + xoffset * sind(-1 * self.wind_map.turbine_wind_direction[i])
                        + coord.x2
                    )

        return x_grid, y_grid, z_grid

"""
        This property returns the axial induction factor of the yawed
        turbine calculated from the coefficient of thrust and the yaw
        angle.

        Returns:
            float: Axial induction factor of a turbine.

        Examples:
            To get the axial induction factor for a turbine:

            >>> aI = floris.farm.turbines[0].aI()
        """
        return (
            0.5
            / cosd(self.yaw_angle)
            * (1 - np.sqrt(1 - self.Ct * cosd(self.yaw_angle)))
        )

def initial_wake_expansion(turbine, U_local, veer, uR, u0):
        yaw = -1 * turbine.yaw_angle 
        sigma_z0 = turbine.rotor_diameter * 0.5 * np.sqrt( uR / (U_local + u0) )
        sigma_y0 = sigma_z0 * cosd(yaw) * cosd(veer)
        return sigma_y0, sigma_z0

beta = 0.5 * ((1 + np.sqrt(1 - Ct)) / np.sqrt(1 - Ct))

        # Calculate sigma_tilde (Eq 9, pp 5 of ref. [1] in docstring)
        sigma_tilde = (self.a_s * TI + self.b_s) * x_tilde + self.c_s * np.sqrt(beta)

        # Calculate n (Eq 13, pp 6 of ref. [1] in docstring)
        n = self.a_f * np.exp(self.b_f * x_tilde) + self.c_f

        # Calculate max vel def (Eq 5, pp 4 of ref. [1] in docstring)
        a1 = 2 ** (2 / n - 1)
        a2 = 2 ** (4 / n - 2)
        C = a1 - np.sqrt(
            a2
            - (
                (n * Ct)
                * cosd(yaw)
                / (
                    16.0
                    * gamma(2 / n)
                    * np.sign(sigma_tilde)
                    * (np.abs(sigma_tilde) ** (4 / n))
                )
            )
        )

        # Compute wake velocity (Eq 1, pp 3 of ref. [1] in docstring)
        velDef1 = U_local * C * np.exp((-1 * r_tilde ** n) / (2 * sigma_tilde ** 2))
        velDef1[x_locations < xR] = 0

        return (
            np.sqrt(velDef1 ** 2),
            np.zeros(np.shape(velDef1)),

Uinf = np.mean(flow_field.wind_map.input_speed) # TODO Is this right?
        
        dist = np.sqrt(yLocs**2 + zLocs**2)
        xLocs = np.abs(x_locations - turbine_coord.x1)
        idx = np.where((dist < D/2) & (xLocs < D/4) & (np.abs(yLocs) > 0.1))

        Gamma = V[idx] * ((2 * np.pi) * (yLocs[idx] ** 2 + zLocs[idx] ** 2)) / (
                yLocs[idx] * (1 - np.exp(-(yLocs[idx] ** 2 + zLocs[idx] ** 2) / ((eps) ** 2))))
        Gamma_wake_rotation = 1.0 * 2 * np.pi * D * (aI - aI ** 2) * turbine.average_velocity / TSR
        Gamma0 = np.mean(np.abs(Gamma))

        test_gamma = np.linspace(-30, 30, 61)
        minYaw = 10000
        for i in range(len(test_gamma)):
            tmp1 = 8 * Gamma0 / (np.pi * flow_field.air_density * D * turbine.average_velocity * Ct)
            tmp = np.abs((sind(test_gamma[i]) * cosd(test_gamma[i]) ** 2) - tmp1)
            if tmp < minYaw:
                minYaw = tmp
                idx = i
        try:
            return test_gamma[idx]
        except:
            print('ERROR',idx)
            return 0.0

def _rotated_grid(self, angle, center_of_rotation):
        """
        Rotate the discrete flow field grid.
        """
        xoffset = self.x - center_of_rotation.x1
        yoffset = self.y - center_of_rotation.x2
        rotated_x = (
            xoffset * cosd(angle) - yoffset * sind(angle) + center_of_rotation.x1
        )
        rotated_y = (
            xoffset * sind(angle) + yoffset * cosd(angle) + center_of_rotation.x2
        )
        return rotated_x, rotated_y, self.z

# determine the effective yaw angle
            yaw_effective = self.effective_yaw(x_locations, y_locations, z_locations, coord, turbine, flow_field)
            yaw = -turbine.yaw_angle  - yaw_effective # opposite sign convention in this model
            print('Effective yaw angle = ', yaw_effective, turbine.yaw_angle)
        else:
            yaw = -turbine.yaw_angle

        tilt = turbine.tilt_angle
        Ct = turbine.Ct

        # U_local = flow_field.wind_map.grid_wind_speed  #just a placeholder for now, should be initialized with the flow_field
        U_local = flow_field.u_initial

        # initial velocity deficits
        uR = U_local * Ct * cosd(tilt) * cosd(yaw) / (
            2. * (1 - np.sqrt(1 - (Ct * cosd(tilt) * cosd(yaw)))))
        u0 = U_local * np.sqrt(1 - Ct)

        # length of near wake
        x0 = D * (cosd(yaw) * (1 + np.sqrt(1 - Ct * cosd(yaw)))) / (
            np.sqrt(2) * (4 * alpha * TI + 2 * beta *
                          (1 - np.sqrt(1 - Ct)))) + coord.x1

        # wake expansion parameters
        ky = ka * TI + kb
        kz = ka * TI + kb

        C0 = 1 - u0 / wind_speed
        M0 = C0 * (2 - C0)
        E0 = C0**2 - 3 * np.exp(1. / 12.) * C0 + 3 * np.exp(1. / 3.)

        # initial Gaussian wake expansion

# turbine parameters
        D = turbine.rotor_diameter
        yaw = -turbine.yaw_angle  # opposite sign convention in this model
        tilt = turbine.tilt_angle
        Ct = turbine.Ct

        # U_local = flow_field.wind_map.grid_wind_speed  #just a placeholder for now, should be initialized with the flow_field
        U_local = flow_field.u_initial

        # initial velocity deficits
        uR = U_local * Ct * cosd(tilt) * cosd(yaw) / (
            2. * (1 - np.sqrt(1 - (Ct * cosd(tilt) * cosd(yaw)))))
        u0 = U_local * np.sqrt(1 - Ct)

        # length of near wake
        x0 = D * (cosd(yaw) * (1 + np.sqrt(1 - Ct * cosd(yaw)))) / (
            np.sqrt(2) * (4 * alpha * TI + 2 * beta *
                          (1 - np.sqrt(1 - Ct)))) + coord.x1

        # wake expansion parameters
        ky = ka * TI + kb
        kz = ka * TI + kb

        C0 = 1 - u0 / wind_speed
        M0 = C0 * (2 - C0)
        E0 = C0**2 - 3 * np.exp(1. / 12.) * C0 + 3 * np.exp(1. / 3.)

        # initial Gaussian wake expansion
        sigma_z0 = D * 0.5 * np.sqrt(uR / (U_local + u0))
        sigma_y0 = sigma_z0 * cosd(yaw) * cosd(veer)

        yR = y_locations - coord.x2