[AstroPy] Blackbody fit and Wien's displacement law

[Rudolf Baer]

In the following code I have done a black body + power law fit of the broad line free continuum of the NIR spectrum of BASS DR1 ID 641
(NGC 4748). The resulting parameters are shown in the  title line of the plot. The plot also shows the two components (blackbody and power law).  I have separately verified that the two components add up to the fit.
However, the lambda derived from Wien’s displacement law for the Temperature derived from black body fit ( 1’995.24 K) does not agree with the maximum of the black body plot.
Do I have a coding error or I am missing something on the science side?
Any comment will be highly appreciated.
R. Baer



from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np

from lmfit import Model
from lmfit import CompositeModel, Model

object='0641'
path_output='/Users/rudolf/Durham/output/'
path_object='baseline_BASS'+object+' 0407.csv'   #<----------change date as nessary
#path_object='baseline_BASS'+object+' 0615.csv'   #<----------change date as nessary
#path_object='baseline_BASS'+object+' 0625.csv'   #<----------change date as nessary

path=path_output+path_object
print ( 'line 19', object, path)


x_AA=  np.genfromtxt(path,delimiter=',',dtype="float64",autostrip=True,skip_header=1,skip_footer=1, usecols=(0))
x=x_AA/10000                                                            #<---------conversion AA to um
#x=x[725:5283]                                                          #<----------partial spectrum

y=  np.genfromtxt(path,delimiter=',',dtype="float64",autostrip=True,skip_header=1,skip_footer=1, usecols=(1))
#y=y[725:5283]                                                         #<----------partial spectrum

#=============================================================================================================
# fit bb + powerlaw

def bb(x, T, const):
    from scipy.constants import h,k,c
    x = 1e-6 * x # convert to metres from um
    return const*2*h*c**2 / (x**5 * (np.exp(h*c / (x*k*T)) - 1)) 

def powerlaw(x,A,p):
    return A*x**p


mod= Model(bb) + Model(powerlaw)
pars  = mod.make_params(T=1600,const=2*1e-21,A=2*np.amin(y),p=-1.0)                          #<-----automatice

result = mod.fit(y,pars,x=x)



print(f"{result.fit_report()[340:387]:<20}")


#============================================================================================================

#Parameters
T=    (result.params['T'].value)
const=(result.params['const'].value)
A=    (result.params['A'].value)
p=    (result.params['p'].value)

# λ Wien
Wien=(2897.77196/T)

#plot bb and PWL separately
#Blackbody

#from scipy.constants import h,k,c
x_bb = 1e-6 * x # convert to metres from um
y_bb=const*2*h*c**2 / (x_bb**5 * (np.exp(h*c / (x_bb*k*T )) - 1)) 

#Power law
y_pwl=A*x**p

plt.title('BB + Powerlaw fitting of continuum of  BASS '+object+str(f"{result.fit_report()[340:588]:<20}"), fontsize='8')
plt.xlabel(' log λ$_{rest}$  [μm]')
plt.ylabel((r' log $ \;\; {\lambda}  \;\; F_{\lambda} \;\;  (\rm erg\;\;s^{-1}\;\;cm^{-2}\;\;)$'))
plt.axis([np.log10(0.7),np.log10 (2.4), np.log10(np.amin(y*x*0.2)),np.log10(np.amax(y*x*1.2))])    #automatic

plt.plot(np.log10(x), np.log10(y_bb*x), 'b--', label='Blackbody component')
plt.plot(np.log10(x), np.log10(y_pwl*x), 'g--', label='Power law component')
plt.plot(np.log10(x),np.log10(y*x),'b-', label='Continumm')                                                  
plt.plot((np.log10(x)), (np.log10(result.best_fit*x)), 'r--', label='Best fit')
plt.vlines(np.log10(Wien),np.amin(np.log10(y*x*0.1)),np.amax(np.log10(y*x*1.2)),colors='red', linestyles='solid', label='Wien')  #Wien's displacement law
plt.legend()
pathfig2='/Users/rudolf/Durham/plots_bb/ '
date='_210725'
type='.png'
pathfig=pathfig2+object+date+type
print (pathfig)

#plt.savefig(pathfig)
plt.show()

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20210725/a90f242e/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot 2021-07-25 at 15.46.25.png
Type: image/png
Size: 84403 bytes
Desc: not available
URL: <https://mail.python.org/pipermail/astropy/attachments/20210725/a90f242e/attachment-0001.png>