keflavich / imf

Simple tools to work with the Initial Mass Function
MIT License
44 stars 17 forks source link

Numerical Integration problems in KoenConvolvedPowerLaw #12

Open teachristian opened 5 years ago

teachristian commented 5 years ago

image

Looks like there is a numerical problem with integration, these plots of pdf do not look very nice. Happened after I vectorized the Koen functions. The plots should have the same shape Generated these plots with:

mmin = .03
mmax = 120
gamma = .9
sigma = .4
m = np.logspace(-3,3)
def integrand(x,y):
    return (x**-(gamma+1)) * np.exp(-.5*((y-x)/sigma)**2)

coef = gamma/((sigma*np.sqrt(2*np.pi)) * ((mmin**-gamma) - (mmax**-gamma)))

def Int(y):
    I = quad(integrand, mmin,mmax,args=(y))[0]   
    return I

vector_I = np.vectorize(Int) 
pbb = coef* vector_I(m)
masses = np.logspace(-3,3)
noerror = imf.KoenTruePowerLaw(.03,120,.9)
pdf_noerror = noerror(masses)

smallerror = imf.KoenConvolvedPowerLaw(.03,120,.9,.05)
pdf_smallerror = smallerror(masses)
pl.loglog(masses,pdf_smallerror)
pl.loglog(masses,pdf_noerror)
pl.loglog(masses,pbb)