<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](https://drive.google.com/file/d/1awDvFpp3YnNi6Xp_TAGqxNg6ADznR3xW/view?usp=sharing)(I have the Luminosity L[Ha] in erg/s).

I would like to plot like this figure from this [paper](https://iopscience.iop.org/article/10.3847/1538-4357/aab620/pdf).

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).

[![enter image description here](https://i.sstatic.net/5osi7.png)](https://i.sstatic.net/5osi7.png)

[![enter image description here](https://i.sstatic.net/gg6fd.png)](https://i.sstatic.net/gg6fd.png)

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

