<p>I would like to plot the Schechter Luminosity function of galaxies. Here is my <a href="https://drive.google.com/file/d/1awDvFpp3YnNi6Xp_TAGqxNg6ADznR3xW/view?usp=sharing" rel="nofollow noreferrer">data</a>(I have the Luminosity L[Ha] in erg/s).</p>
<p>I would like to plot like this figure from this <a href="https://iopscience.iop.org/article/10.3847/1538-4357/aab620/pdf" rel="nofollow noreferrer">paper</a>.</p>
<p>In this paper, they used these equations for the plotting(there are four datasets in this case, but for my case, I want to plot only my data).</p>
<p><a href="https://i.sstatic.net/5osi7.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/5osi7.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/gg6fd.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gg6fd.png" alt="enter image description here" /></a></p>
<p>I used python to plot the function. Here is my code. But, I am not able to plot it. Kindly do help me with this. Also How to do the fitting with my data to get the required parameters.</p>
<pre><code>import numpy as np, math
import matplotlib.pyplot as plt
from astropy.table import Table
from matplotlib.colors import LogNorm
from scipy.interpolate import interp1d

data=Table.read('lum.fits')

def Schechter_Lum(L, L_star=None, Alpha=None, P_star=None):
'''
Schechter luminosity function in LUMINOSITY  
number per dex
'''
ratio = L/L_star
f = P_star* np.exp**(-ratio) *(ratio)**Alpha *ratio*np.log(10.) 
return(f)

#For the reference, I used this value from the paper
L_star=42.64
P_star=-3.32
Alpha= -1.75

schechter_lf=Schechter_Lum(data['ha_luminosity_log'],L_star,Alpha,P_star)

plt.plot(data['ha_luminosity_log'],schechter_lf)
plt.yscale('log')
plt.title('H<span class="math-container">$\alpha$</span> Luminosity Functions')
plt.xlabel(r'Log (L[H<span class="math-container">$\alpha$</span>]) erg/s', fontsize=14)
plt.ylabel(r'<span class="math-container">$\varphi$</span>(logL) per <span class="math-container">$Mpc$</span>^{-3}$ per logL',fontsize=14)
plt.tight_layout()
plt.show()
</code></pre>
<p>I am not able to get the figure properly and I am not able to perform fitting of my Luminosity function. Kindly do help with this. Thank You.</p>

CC BY-SA 4.0
2021-03-16T14:41:08.380Z
Without having tested, I think your L and L* should be in linear units, not log units. I.e. L* = 10**42.64, not 42.64.
pela
CC BY-SA 4.0
2021-03-17T05:02:34.820Z
@pela Thank You. I tried it but it is not getting a proper result like the plot. Kindly do help if you know the solution.
Edward
CC BY-SA 4.0
2021-03-17T10:56:52.743Z
Your syntax `np.exp**(-ratio)` is wrong, and should be `np.exp(-ratio)` (you're trying to lift a _function_ to some power, rather than the number $e$). If that's not it, then maybe add a little more info about why you're not able to plot it, i.e. what happens when you try, can you compile, etc.
pela
<galaxy><data-analysis><luminosity><python>
Schechter Luminosity Function fitting plot of galaxies
Edward
# Schechter Luminosity Function fitting plot of galaxies

(Asked by: username_0 on 2021-03-16)

I would like to plot the Schechter Luminosity function of galaxies. Here is my data(I have the Luminosity L[Ha] in erg/s).
I would like to plot like this figure from this paper.
In this paper, they used these equations for the plotting(there are four datasets in this case, but for my case, I want to plot only my data).


I used python to plot the function. Here is my code. But, I am not able to plot it. Kindly do help me with this. Also How to do the fitting with my data to get the required parameters.
import numpy as np, math
import matplotlib.pyplot as plt
from astropy.table import Table
from matplotlib.colors import LogNorm
from scipy.interpolate import interp1d

data=Table.read('lum.fits')

def Schechter_Lum(L, L_star=None, Alpha=None, P_star=None):
'''
Schechter luminosity function in LUMINOSITY  
number per dex
'''
ratio = L/L_star
f = P_star* np.exp**(-ratio) *(ratio)**Alpha *ratio*np.log(10.) 
return(f)

#For the reference, I used this value from the paper
L_star=42.64
P_star=-3.32
Alpha= -1.75

schechter_lf=Schechter_Lum(data['ha_luminosity_log'],L_star,Alpha,P_star)

plt.plot(data['ha_luminosity_log'],schechter_lf)
plt.yscale('log')
plt.title('H$\alpha$ Luminosity Functions')
plt.xlabel(r'Log (L[H$\alpha$]) erg/s', fontsize=14)
plt.ylabel(r'$\varphi$(logL) per $Mpc$^{-3}$ per logL',fontsize=14)
plt.tight_layout()
plt.show()

I am not able to get the figure properly and I am not able to perform fitting of my Luminosity function. Kindly do help with this. Thank You.

**Tags:** galaxy, data-analysis, luminosity, python

> **Comments:**
> **username_1**: Without having tested, I think your L and L* should be in linear units, not log units. I.e. L* = 10**42.64, not 42.64.  
> **username_0**: @username_1 Thank You. I tried it but it is not getting a proper result like the plot. Kindly do help if you know the solution.  
> **username_1**: Your syntax `np.exp**(-ratio)` is wrong, and should be `np.exp(-ratio)` (you're trying to lift a _function_ to some power, rather than the number $e$). If that's not it, then maybe add a little more info about why you're not able to plot it, i.e. what happens when you try, can you compile, etc.  

---

## 0 answers

