diff --git a/readme.md b/readme.md
index 783cfdc..eeede75 100644
--- a/readme.md
+++ b/readme.md
@@ -57,7 +57,7 @@ flux = star.flux(phase)
 ```python
 plt.figure(figsize=(9, 3))
 plt.subplot(1, 5, (1, 2))
-star.plot()
+star.show()
 
 plt.subplot(1, 5, (3, 5))
 plt.plot(phase, flux, c="k")
diff --git a/spotter/star.py b/spotter/star.py
index dc70624..4548977 100644
--- a/spotter/star.py
+++ b/spotter/star.py
@@ -3,15 +3,28 @@
 import numpy as np
 
 
+def _wrap(*args):
+    n = len(args)
+    signature = ",".join(["()"] * n)
+    signature = f"{signature}->{signature}"
+    new_args = np.vectorize(lambda *args: args)(*args)
+    if new_args[0].shape == ():
+        return [np.array([n]) for n in new_args]
+    else:
+        return new_args
+
+
 class Star:
     def __init__(self, u=None, N=64):
         self.N = N
         self.u = u
         self.n = hp.nside2npix(self.N)
-        self._m = np.ones(self.n)
         self._phis, self._thetas = hp.pix2ang(self.N, np.arange(self.n))
         self._sin_phi = np.sin(self._phis)
 
+        self._spot_map = np.zeros(self.n)
+        self._faculae_map = np.zeros(self.n)
+
     def _z(self, phase=0):
         return self._sin_phi * np.cos(self._thetas - phase)
 
@@ -30,30 +43,57 @@ def _get_mask(self, phase=0):
         else:
             mask = (self._thetas > b) | (self._thetas < a)
 
-        return mask * self._ld(phase)
+        return mask
 
     def add_spot(self, theta, phi, radius, contrast):
-        @np.vectorize(signature="(),(),(),()->(),(),(),()")
-        def foo(a, b, c, d):
-            return a, b, c, d
-
-        for t, p, r, c in zip(*foo(theta, phi, radius, contrast)):
+        for t, p, r, c in zip(*_wrap(theta, phi, radius, contrast)):
             idxs = hp.query_disc(self.N, hp.ang2vec(t, p), r)
-            self._m[idxs] = 1 - c
+            self._spot_map[idxs] = c
+
+    def add_faculae(self, theta, phi, radius_in, radius_out, contrast):
+        for t, p, ri, ro, c in zip(*_wrap(theta, phi, radius_in, radius_out, contrast)):
+            inner_idxs = hp.query_disc(self.N, hp.ang2vec(t, p), ri)
+            outer_idxs = hp.query_disc(self.N, hp.ang2vec(t, p), ro)
+            idxs = np.setdiff1d(outer_idxs, inner_idxs)
+            self._faculae_map[idxs] = c
+
+    def add_spot_faculae(
+        self, theta, phi, radius_in, radius_out, contrast_spot, contrast_faculae
+    ):
+        for t, p, ri, ro, cs, cf in zip(
+            *_wrap(theta, phi, radius_in, radius_out, contrast_spot, contrast_faculae)
+        ):
+            inner_idxs = hp.query_disc(self.N, hp.ang2vec(t, p), ri)
+            outer_idxs = hp.query_disc(self.N, hp.ang2vec(t, p), ro)
+            facuale_idxs = np.setdiff1d(outer_idxs, inner_idxs)
+            self._faculae_map[facuale_idxs] = cf
+            self._spot_map[inner_idxs] = cs
 
     def flux(self, phase=0):
         def _flux(phase):
             mask = self._get_mask(phase)
-            return (self._m * mask).sum() / mask.sum()
+            limb_darkening = self._ld(phase)
+            # spot contribution
+            m = (1 - self._spot_map) * mask * limb_darkening
+            # facuale contribution will have a different limb darkening
+            # and another normalization will be needed
+            return m.sum() / (mask * limb_darkening).sum()
 
         return np.vectorize(_flux)(phase)
 
     def m(self, phase=0):
-        return self._m * self._get_mask(phase)
+        mask = self._get_mask(phase)
+        limb_darkening = self._ld(phase)
+        # spot contribution
+        m = 1 - self._spot_map * mask * limb_darkening
+        return m
 
-    def plot(self, phase=0, grid=False, return_img=False, **kwargs):
+    def show(self, phase=0, grid=False, return_img=False, **kwargs):
         kwargs.setdefault("cmap", "magma")
-        rotated_m = hp.Rotator(rot=[phase, 0], deg=False).rotate_map_pixel(self._m)
+        # only spot contribution for now
+        rotated_m = hp.Rotator(rot=[phase, 0], deg=False).rotate_map_pixel(
+            1 - self._spot_map
+        )
         projected_map = hp.orthview(
             rotated_m * self._ld(0), half_sky=True, return_projected_map=True
         )
@@ -63,3 +103,31 @@ def plot(self, phase=0, grid=False, return_img=False, **kwargs):
         else:
             plt.axis(False)
             plt.imshow(projected_map, **kwargs)
+
+    def covering_fraction(
+        self, phase: float = None, vmin: float = 0.01, vmax: float = 1.0
+    ):
+        """Return the covering fraction of active regions
+
+        Either computed for the whole star (`phase=None`) or for the stellar
+        disk given a phase
+
+        Parameters
+        ----------
+        phase : float, optional
+            stellar rotation phase, by default None
+        vmin : float, optional
+            minimum contrast value for spots, by default 0.01
+        vmax : float, optional
+            minimum contrast value for faculae, by default 1.0
+
+        Returns
+        -------
+        float
+            full star or disk covering fraction
+        """
+        if phase is None:
+            return np.sum(self._spot_map >= vmin) / self.n
+        else:
+            mask = self._get_mask(phase)
+            return np.sum(self._spot_map[mask] >= vmin) / mask.sum()