diff --git a/slsim/Sources/Supernovae/supernovae_pop.py b/slsim/Sources/Supernovae/supernovae_pop.py
index fa4572790..7c7f20746 100644
--- a/slsim/Sources/Supernovae/supernovae_pop.py
+++ b/slsim/Sources/Supernovae/supernovae_pop.py
@@ -80,6 +80,7 @@ def calculate_SNIa_rate(self, z, eta=0.04):
         """Calculates the rate of SN Ia. (Eq 15 - Oguri and Marshall 2010)
 
         :param z: redshift (z>=0)
+        :type z: array-like
         :param eta: canonical efficiency
 
         :return: SN Ia rate n(z) in [(h)yr^(-1)Mpc^(-3)]
@@ -99,3 +100,23 @@ def calculate_SNIa_rate(self, z, eta=0.04):
             SNIa_rate_list.append(SNIa_rate)
 
         return np.array(SNIa_rate_list)
+
+    def calculate_luminosity(self, M, z):
+        """Calculates the luminosity function of SN Ia. (Eq 18 - Oguri and Marshall 2010)
+
+        :param M: B-band absolute magnitude
+        :param z: redshift (z>=0)
+        :type z: array-like
+
+        :return: luminosity function dPhi_dM in [mag^(-1)Mpc^(-3)]
+        :return type: array-like
+        """
+        n = self.calculate_SNIa_rate(z)
+        M_star = -19.06
+        sigma = 0.56
+
+        dPhi_dM = (n * np.e ** (-((M - M_star) ** 2) / (2 * sigma**2))) / (
+            (1 + z) * np.sqrt(2 * np.pi) * sigma
+        )
+
+        return dPhi_dM
diff --git a/tests/test_Source/test_supernovae_pop.py b/tests/test_Source/test_supernovae_pop.py
index 791073007..a3afc3090 100644
--- a/tests/test_Source/test_supernovae_pop.py
+++ b/tests/test_Source/test_supernovae_pop.py
@@ -6,6 +6,7 @@
 from astropy.cosmology import FlatLambdaCDM
 import numpy.testing as npt
 import pytest
+import numpy as np
 
 
 def test_calculate_star_formation_rate():
@@ -60,13 +61,23 @@ def test_numerator_integrand(self):
 
     def test_calculate_SNIa_rate(self):
         # (Fig 2 - Oguri and Marshall 2010)
-        z_array = [0, 1, 2, 3]
-        rate_array = self.sne_rate.calculate_SNIa_rate(z_array)
+        z = np.array([0, 1, 2, 3])
+        rate = self.sne_rate.calculate_SNIa_rate(z)
 
-        npt.assert_almost_equal(rate_array[0], 0.000041006, decimal=3)
-        npt.assert_almost_equal(rate_array[1], 0.0001191, decimal=3)
-        npt.assert_almost_equal(rate_array[2], 0.0001349, decimal=3)
-        npt.assert_almost_equal(rate_array[3], 0.00008008, decimal=3)
+        npt.assert_almost_equal(rate[0], 0.000041006, decimal=3)
+        npt.assert_almost_equal(rate[1], 0.0001191, decimal=3)
+        npt.assert_almost_equal(rate[2], 0.0001349, decimal=3)
+        npt.assert_almost_equal(rate[3], 0.00008008, decimal=3)
+
+    def test_calculate_luminosity(self):
+        M = -28
+        z = np.array([1, 2, 4, 6])
+        luminosity = self.sne_rate.calculate_luminosity(M, z)
+
+        npt.assert_almost_equal(luminosity[0], 1.93040e-60, decimal=3)
+        npt.assert_almost_equal(luminosity[1], 1.45804e-60, decimal=3)
+        npt.assert_almost_equal(luminosity[2], 1.95788e-61, decimal=3)
+        npt.assert_almost_equal(luminosity[3], 1.89669e-62, decimal=3)
 
 
 if __name__ == "__main__":