diff --git a/nbodykit/cosmology/power/__init__.py b/nbodykit/cosmology/power/__init__.py
index 1d7ecf017..7eacc160c 100644
--- a/nbodykit/cosmology/power/__init__.py
+++ b/nbodykit/cosmology/power/__init__.py
@@ -1,4 +1,5 @@
+from .galaxy_opt import FNLGalaxyPower
 from .linear import LinearPower, EHPower, NoWiggleEHPower
 from .zeldovich import ZeldovichPower
 from .halofit import HalofitPower
-from . import transfers
+from . import transfers
\ No newline at end of file
diff --git a/nbodykit/cosmology/power/galaxy.py b/nbodykit/cosmology/power/galaxy.py
new file mode 100644
index 000000000..bbc3861b3
--- /dev/null
+++ b/nbodykit/cosmology/power/galaxy.py
@@ -0,0 +1,129 @@
+import numpy
+from . import transfers
+from ..cosmology import Cosmology
+from .linear import LinearPower
+
+class FNLGalaxyPower(object):
+    """
+    An object to compute the galaxy/redshift power spectrum and related quantities,
+    using a transfer function from the CLASS code or the analytic
+    Eisenstein & Hu approximation.
+
+    Parameters
+    ----------
+    cosmo : :class:`Cosmology`, astropy.cosmology.FLRW
+        the cosmology instance; astropy cosmology objects are automatically
+        converted to the proper type
+    redshift : float
+        the redshift of the power spectrum
+    transfer : str, optional
+        string specifying the transfer function to use; one of
+        'CLASS', 'EisensteinHu', 'NoWiggleEisensteinHu'
+    b0 : float
+        the linear bias of the galaxy in a gaussian field
+    fnl : float
+        primordial non-gaussian parameter
+    p : float
+        correction term due to galaxy selection beyond a Poisson sampling of the halos of a given mass.
+        p takes a between 1 and 1.6. 
+        p=1 for objects populating a fair sample of all the halos
+        p=1.6 for objects that populate only recently merged halos
+    Omega_m : float
+        the matter density parameter
+    H0 : float
+        the Hubble constant (in units of km/s Mpc/h)
+    c : float
+        speed of light (in units of km/s)
+        
+
+    Attributes
+    ----------
+    cosmo : class:`Cosmology`
+        the object giving the cosmological parameters
+    sigma8 : float
+        the z=0 amplitude of matter fluctuations
+    redshift : float
+        the redshift to compute the power at
+    transfer : str
+        the type of transfer function used
+    """
+    
+    def __init__(self, cosmo, redshift ,b0 ,fNL ,p ,c=3e5 ,transfer='CLASS'):
+        from astropy.cosmology import FLRW
+
+        # convert astropy
+        if isinstance(cosmo, FLRW):
+            from nbodykit.cosmology import Cosmology
+            cosmo = Cosmology.from_astropy(cosmo)
+
+        # store a copy of the cosmology
+        self.cosmo = cosmo.clone()
+        
+        #get the linear bias,p,fNL,omega_m,H0,c
+        self.b=b0
+        self.p=p
+        self.fnl=fNL
+        self.omega_m=cosmo.Omega0_m
+        self.H0=cosmo.H0
+        self.c=c
+
+        self.transfer = transfer
+
+        # set redshift
+        self.redshift = redshift
+        
+    
+    def total_bias(self, k):
+        """
+        Returns the total galaxy bias in a non-gaussian field at the redshift
+        specified by :attr:`redshift`.
+        
+        Parameters
+        ---------
+        k : float, array_like
+            the wavenumber in units of :math:`h Mpc^{-1}`
+
+        Returns
+        -------
+        total_bias : float, array_like
+            the total bias
+        
+        """
+        Plin=LinearPower(self.cosmo, self.redshift, transfer=self.transfer)
+        Pk=Plin(k)
+        crit_density=1.686
+        Dz=1/(1+self.redshift)         
+        del_b= 3*self.fnl*(self.b-self.p)* crit_density * self.omega_m/(k**2*Pk**0.5*Dz) * (self.H0/self.c)**2
+                
+        total_bias=self.b + del_b
+        
+        return total_bias
+
+    def __call__(self, k):
+        """
+        Return the galaxy power spectrum in units of
+        :math:`h^{-3} \mathrm{Mpc}^3` at the redshift specified by
+        :attr:`redshift`.
+
+        The transfer function used to evaluate the power spectrum is
+        specified by the ``transfer`` attribute.
+
+        Parameters
+        ---------
+        k : float, array_like
+            the wavenumber in units of :math:`h Mpc^{-1}`
+
+        Returns
+        -------
+        Pk : float, array_like
+            the galaxy power spectrum evaluated at ``k`` in units of
+            :math:`h^{-3} \mathrm{Mpc}^3`
+        """
+
+        Plin=LinearPower(self.cosmo, self.redshift, transfer=self.transfer)
+        Pk=Plin(k)
+        total_bias=self.total_bias(k)      
+        
+        Pgal = Pk * total_bias**2
+        
+        return Pgal
diff --git a/nbodykit/cosmology/tests/test_power.py b/nbodykit/cosmology/tests/test_power.py
index 5df44fd1a..20d05fd52 100644
--- a/nbodykit/cosmology/tests/test_power.py
+++ b/nbodykit/cosmology/tests/test_power.py
@@ -1,7 +1,7 @@
-from nbodykit.cosmology import Cosmology, LinearPower, HalofitPower, ZeldovichPower
+from nbodykit.cosmology import Cosmology, LinearPower, HalofitPower, ZeldovichPower ,FNLGalaxyPower
 from nbodykit.cosmology import EHPower, NoWiggleEHPower
 import numpy
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_equal
 import pytest
 
 def test_bad_transfer():
@@ -160,3 +160,45 @@ def test_zeldovich():
     D2 = c.scale_independent_growth_factor(0.)
     D3 = c.scale_independent_growth_factor(0.55)
     assert_allclose(Pk2.max()/Pk3.max(), (D2/D3)**2, rtol=1e-2)
+    
+def test_galaxy():
+    
+    # initialize the power
+    c = Cosmology().match(sigma8=0.82)
+    P = FNLGalaxyPower(c, redshift=0, b0=1, fNL=10, p=1, c=3e5, transfer='CLASS')
+
+    # desired wavenumbers (in h/Mpc)
+    k = numpy.logspace(-3, 2, 500)
+
+    # initialize EH power
+    P1 = FNLGalaxyPower(c, redshift=0, b0=1, fNL=10, p=1, c=3e5, transfer="CLASS")
+    P2 = FNLGalaxyPower(c, redshift=0, b0=1, fNL=10, p=1, c=3e5, transfer='EisensteinHu')
+    P3 = FNLGalaxyPower(c, redshift=0, b0=1, fNL=10, p=1, c=3e5, transfer='NoWiggleEisensteinHu')
+
+    # check different transfers (very roughly)
+    Pk1 = P1(k)
+    Pk2 = P2(k)
+    Pk3 = P3(k)
+    assert_allclose(Pk1 / Pk1.max(), Pk2 / Pk2.max(), rtol=0.1)
+    assert_allclose(Pk1 / Pk1.max(), Pk3 / Pk3.max(), rtol=0.1)
+
+    # also try scalar
+    Pk = P(0.1)
+    
+    # check if bias and total bias are same for a gaussian field (i.e.,fnl=0) for any k
+    bias=1
+    P = FNLGalaxyPower(c, redshift=0, b0=bias, fNL=0, p=1, c=3e5, transfer='CLASS')
+    assert_equal(P.total_bias(0.1),bias)
+    assert_equal(P.total_bias(0.3),bias)
+    
+def test_gaussian():
+    
+    # initialize linear power
+    c = Cosmology().match(sigma8=0.82)
+    Plin = LinearPower(c, redshift=0, transfer='CLASS')
+    
+    # initialize galaxy power for a gaussian field (i.e.,fnl=0) with bias=1
+    Pgal = FNLGalaxyPower(c, redshift=0, b0=1, fNL=0, p=1, c=3e5, transfer='CLASS')
+    
+    # check if they are equal
+    assert_equal(Pgal(0.1),Plin(0.1))
\ No newline at end of file