notebook.community

Edit and run



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline



In [ ]:

    
from newdust import scatmodels
from newdust.graindist import composition

from scipy.integrate import trapz



In [ ]:

    
ENERGY = np.logspace(-1,2,100) # keV
AUM    = np.array(1.0) # um
AUM_CM = AUM * (1.e-6 *100.) # cm
CM     = composition.CmDrude()
CMS    = composition.CmSilicate()
#THETA  = np.logspace(-10.0, np.log10(np.pi), 1000) # radians
#THETA_ASEC = THETA * (360.*60.*60)/(2.*np.pi)

VLAM  = 4500. # angs

THETA_RAD  = np.logspace(-5., np.log10(np.pi), 1000)
THETA_ASEC = THETA_RAD * (360.0*60.*60.) / (2.0*np.pi)



In [ ]:

    
rgd = scatmodels.RGscat()
rgd.calculate(ENERGY, AUM, CM, unit='kev')



In [ ]:

    
plt.plot(ENERGY, rgd.qsca)
plt.loglog()



In [ ]:

    
rgd2 = scatmodels.RGscat()
rgd2.calculate(ENERGY, AUM, CM, unit='kev', theta=THETA_ASEC)



In [ ]:

    
np.shape(rgd2.qsca)



In [ ]:

    
np.shape(rgd2.diff)



In [ ]:

    
plt.plot(THETA_ASEC, rgd2.diff[0,0,:], 'b-', lw=2)
plt.plot(THETA_ASEC, rgd2.diff[-1,0,:], 'k--', lw=2)
plt.loglog()



In [ ]:

    
from scipy.integrate import trapz



In [ ]:

    
sigma_sca = rgd2.qsca[0,0] * np.pi * AUM_CM**2  # cm^2
test = trapz(rgd2.diff[0,0,:] * 2.0*np.pi*np.sin(THETA_RAD), THETA_RAD)



In [ ]:

    
print(test/sigma_sca)

Mie model



In [ ]:

    
mtest = scatmodels.Mie()
mtest.calculate(VLAM, AUM, CMS, unit='angs', theta=THETA_ASEC)



In [ ]:

    
np.shape(mtest.diff)



In [ ]:

    
plt.plot(THETA_ASEC, mtest.diff[0,0,:])
plt.loglog()
plt.xlim(1.e5, 1.e6)



In [ ]:

    
trapz(mtest.diff * 2.0*np.pi*np.sin(THETA_RAD), THETA_RAD)



In [ ]:

    
mtest.qsca * np.pi * (AUM_CM)**2

Test multi-dimensional input



In [ ]:

    
NE, NA = 2, 20
LAMVALS = np.linspace(1000.,5000.,NE)  # angs
AVALS   = np.linspace(0.1, 0.5, NA)    # um



In [ ]:

    
mtest2 = scatmodels.Mie()
mtest2.calculate(LAMVALS, AVALS, CMS, unit='angs', theta=THETA_ASEC)



In [ ]:

    
plt.plot(THETA_ASEC, mtest2.diff[0,0,:])
plt.plot(THETA_ASEC, mtest2.diff[-1,-1,:])
plt.semilogy()



In [ ]:

    
trapz(mtest2.diff[0,0,:] * 2.0*np.pi*np.sin(THETA_RAD), THETA_RAD)



In [ ]:

    
mtest2.qsca[0,0] * np.pi * (AVALS[0] * 1.e-4)**2



In [ ]:

    
trapz(mtest2.diff[-1,-1,:] * 2.0*np.pi*np.sin(THETA_RAD), THETA_RAD)



In [ ]:

    
mtest2.qsca[-1,-1] * np.pi * (AVALS[-1] * 1.e-4)**2



In [ ]:

    
import newdust.constants as c



In [ ]:

    
c.micron2cm



In [ ]: