def METHOD_NAME(self, n):
    """Generate semi-major axis (AU), eccentricity, geometric albedo, and
    planetary radius (earthRad)
    Semi-major axis and planetary radius are jointly distributed.
    Eccentricity is a Rayleigh distribution. Albedo is a constant value
    based on planetary radius.
    Args:
        n (integer):
            Number of samples to generate
    Returns:
        tuple:
        a (astropy Quantity array):
            Semi-major axis in units of AU
        e (float ndarray):
            Eccentricity
        p (float ndarray):
            Geometric albedo
        Rp (astropy Quantity array):
            Planetary radius in units of earthRad
    """
    n = self.gen_input_check(n)
    # generate semi-major axis and planetary radius samples
    a, Rp = self.gen_sma_radius(n)
    # check for constrainOrbits == True for eccentricity samples
    # constants
    C1 = np.exp(-self.erange[0] ** 2 / (2.0 * self.esigma**2))
    ar = self.arange.to("AU").value
    if self.constrainOrbits:
        # restrict semi-major axis limits
        arcon = np.array(
            [ar[0] / (1.0 - self.erange[0]), ar[1] / (1.0 + self.erange[0])]
        )
        # clip sma values to sma range
        sma = np.clip(a.to("AU").value, arcon[0], arcon[1])
        # upper limit for eccentricity given sma
        elim = np.zeros(len(sma))
        amean = np.mean(ar)
        elim[sma <= amean] = 1.0 - ar[0] / sma[sma <= amean]
        elim[sma > amean] = ar[1] / sma[sma > amean] - 1.0
        elim[elim > self.erange[1]] = self.erange[1]
        elim[elim < self.erange[0]] = self.erange[0]
        # additional constant
        C2 = C1 - np.exp(-(elim**2) / (2.0 * self.esigma**2))
        a = sma * u.AU
    else:
        C2 = self.enorm
    e = self.esigma * np.sqrt(-2.0 * np.log(C1 - C2 * np.random.uniform(size=n)))
    # generate albedo from planetary radius
    p = self.get_p_from_Rp(Rp)
    return a, e, p, Rp