<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def geoid(self, potref, a=None, f=None, r=None, omega=None, order=2, lmax=None, lmax_calc=None, grid='DH2'):
 """ Create a global map of the height of the geoid and return an SHGeoid class instance. Usage ----- geoid = x.geoid(potref, [a, f, r, omega, order, lmax, lmax_calc, grid]) Returns ------- geoid : SHGeoid class instance. Parameters potref : float The value of the potential on the chosen geoid, in m2 / s2. a : optional, float, default = self.r0 The semi-major axis of the flattened ellipsoid on which the field is computed. f : optional, float, default = 0 The flattening of the reference ellipsoid: f=(a-b)/a. r : optional, float, default = self.r0 The radius of the reference sphere that the Taylor expansion of the potential is calculated on. order : optional, integer, default = 2 The order of the Taylor series expansion of the potential about the reference radius r. This can be either 1, 2, or 3. omega : optional, float, default = self.omega The angular rotation rate of the planet. lmax : optional, integer, default = self.lmax The maximum spherical harmonic degree that determines the number of samples of the output grid, n=2lmax+2, and the latitudinal sampling interval, 90/(lmax+1). lmax_calc : optional, integer, default = lmax The maximum spherical harmonic degree used in evaluating the functions. This must be less than or equal to lmax. grid : str, optional, default = 'DH2' 'DH' or 'DH1' for an equisampled lat/lon grid with nlat=nlon, or 'DH2' for an equidistant lat/lon grid with nlon=2*nlat. Description This method will create a global map of the geoid height, accurate to either first, second, or third order, using the method described in Wieczorek (2007; equation 19-20). The algorithm expands the potential in a Taylor series on a spherical interface of radius r, and computes the height above this interface to the potential potref exactly from the linear, quadratic, or cubic equation at each grid point. If the optional parameters a and f are specified, the geoid height will be referenced to a flattened ellipsoid with semi-major axis a and flattening f. The pseudo-rotational potential is explicitly accounted for by using the angular rotation rate omega of the planet in the SHGravCoeffs class instance. If omega is explicitly specified for this method, it will override the value present in the class instance. Reference Wieczorek, M. A. Gravity and topography of the terrestrial planets, Treatise on Geophysics, 10, 165-206, 2007. """

        if a is None:
            a = self.r0
        if f is None:
            f = 0.
        if r is None:
            r = self.r0
        if lmax is None:
            lmax = self.lmax
        if lmax_calc is None:
            lmax_calc = lmax

        if grid.upper() in ('DH', 'DH1'):
            sampling = 1
        elif grid.upper() == 'DH2':
            sampling = 2
        else:
            raise ValueError(
                    "grid must be 'DH', 'DH1', or 'DH2'. "
                    "Input value was {:s}".format(repr(grid)))

        if self.errors is not None:
            coeffs, errors = self.to_array(normalization='4pi', csphase=1)
        else:
            coeffs = self.to_array(normalization='4pi', csphase=1)

        if omega is None:
            omega = self.omega

        geoid = _MakeGeoidGridDH(coeffs, self.r0, self.gm, potref, lmax=lmax,
                                 omega=omega, r=r, order=order,
                                 lmax_calc=lmax_calc, a=a, f=f,
                                 sampling=sampling)

        return _SHGeoid(geoid, self.gm, potref, a, f, omega, r, order,
                        lmax, lmax_calc)