Commit 8d030fe2 authored by Martin Glatzle's avatar Martin Glatzle

Merge branch 'dev' into 'master'

Preparation for Dusting the Universe

See merge request !2
parents d61c15e8 2fbc72e7
......@@ -2,6 +2,7 @@
*.pyc
# emacs backup files
*~
\#*\#
# python egg info
*.egg-info/*
.cache/*
......
The contents of the files in the "data" directory have been obtained from the
scientific literature and are the intelectual property of their respective
authors. If you are the author of the contents of one of those files and feel
you have not been properly credited or have any other concern, please contact
mglatzle@mpa-garching.mpg.de.
All other files in this repository are distributed under GPLv3 as detailed
below.
GNU GENERAL PUBLIC LICENSE
Version 3, 29 June 2007
......
# A toolbox for cosmic dust.
# `cosmic_dustbox`
A package that aims at bringing models of cosmic dust to Python to make working
with them easier.
*Note:* very early development at the moment so it does not yet do too many
useful things and the API _will_ change.
## Installation
- Clone the repo: `git clone
git@gitlab.mpcdf.mpg.de:mglatzle/cosmic_dustbox.git`
- `cd cosmic_dustbox`
- `python setup.py install`
## Examples
Note: `cosmic_dustbox` uses `astropy.units` throughout.
Evaluating the Weingartner & Draine 2001 size distributions for carbonaceous
and silicate grains for R\_V=3.1 and b\_C=6.0e-5 in case 'A' at a range of
grain sizes between 1 nm and 1 micron:
```
>>> import astropy.units as u
>>> import cosmic_dustbox as cb
>>> import numpy as np
>>> sd_car, sd_sil = cb.sdist.WD01(3.1, 6.0, 'A')
>>> sizes = np.logspace(-9, -6, num=10)*u.m
>>> sd_car(sizes)
<Quantity [7.39939133e+00, 9.02469859e-02, 1.94996300e-02, 1.07158763e-03,
1.18483502e-04, 1.30850238e-05, 1.24648409e-06, 9.43835731e-08,
2.41782877e-09, 2.43258206e-15] 1 / m>
>>> sd_sil(sizes)
<Quantity [7.86318728e+000, 6.70676535e-001, 5.73425242e-002,
4.92807748e-003, 4.28146066e-004, 3.80362296e-005,
3.52990312e-006, 3.09011411e-007, 7.19725776e-020,
9.29346597e-263] 1 / m>
```
## Philosophy
In 1930, after confirming that distant stars in the Milky Way are dimmed more
strongly than one would expect from geometry, Trumpler conjectured the
existence of "dust particles of various sizes" in the interstellar space of our
......@@ -18,8 +55,8 @@ for sure, but we'll try our best.
To introduce the objects central to how we describe cosmic dust, assume we are
interested in the dust in a certain environment. Assume then that we cut some
volume out of said environment. The volume should, of course, be representative
the environment. So if we are interested in the dust in the Solar System as a
whole, a cubic meter of space somewhere between Earth and Mars will be too
of the environment. So if we are interested in the dust in the Solar System as
a whole, a cubic meter of space somewhere between Earth and Mars will be too
small and a box the size of the Milky Way will be too large. Suppose we
collected all of the dust grains contained in this volume. We would like to be
able to describe this collection with a reasonable degree of accuracy without
......@@ -32,14 +69,14 @@ other grains, grains that have similar chemical composition and so on. So if we
are fine with some level of abstraction and generalization, our grain
collection could be represented like this:
![](figs/grain-pop.png)
![](docs/figs/grain-pop.png)
with color-coded chemical composition. How do we go about categorizing these
grains so that we can think in these categories without worrying about
individual grains? Two obvious properties to group by are the chemical
composition and the shape:
![](figs/grain-species.png)
![](docs/figs/grain-species.png)
The "blueprint grain", which we can rescale to obtain any of the grains in one
of the groups, is referred to as a grain species. So one of our groups is
......@@ -58,8 +95,10 @@ So, to recap, what we need to describe a collection of grains are:
This approach can be refined as needed by simply splitting grain species. In
the extreme, we would have one species per grain and each species would have a
delta-peak size distribution and we would be back to the overwhelming diversity
we started out from.
we started out from. With this package, we aim to translate the approach to
software. The size distribution class in `sdist.py` is already in a somewhat
advanced state so feel free to check it out.
# Literature
"Physics of the Interstellar and Intergalactic Medium" by B. Draine provides
a good introduction (and more than that) to cosmic dust.
## Literature
"Physics of the Interstellar and Intergalactic Medium" by B. Draine provides a
good introduction (and much more than that) to cosmic dust.
......@@ -22,11 +22,6 @@ class Crefin(object):
array shall correspond to grain size and the second axis to
wavelength. Note that this axis ordering is reversed wrt to the output
of `np.meshgrid`.
Attributes
----------
_f : callable
Directly taken from parameters.
"""
def __init__(self, f):
......@@ -146,9 +141,10 @@ class SGPAHCrefin(Crefin):
@classmethod
def parseCrefinFile(cls, path):
with open(path) as f:
data = _np.loadtxt(f, skiprows=5)
data = _np.loadtxt(f, skiprows=6)
f.seek(0)
next(f)
# skip tow lines
[next(f) for i in range(2)]
size = float(f.readline().split('=')[0])
return size, data
......
Downloaded from B. Draine's website at https://www.astro.princeton.edu/~draine/dust/dust.diel.html
ICOMP=18: Draine 2003 graphite E parallel to c
0.0100 = grain radius (micron)
20.00 = grain temperature (K)
......
Downloaded from B. Draine's website at https://www.astro.princeton.edu/~draine/dust/dust.diel.html
ICOMP=18: Draine 2003 graphite E parallel to c
0.1000 = grain radius (micron)
20.00 = grain temperature (K)
......
Downloaded from B. Draine's website at https://www.astro.princeton.edu/~draine/dust/dust.diel.html
ICOMP=19: Draine 2003 graphite E perp. c
0.0100 = grain radius (micron)
20.00 = grain temperature (K)
......
Downloaded from B. Draine's website at https://www.astro.princeton.edu/~draine/dust/dust.diel.html
ICOMP=19: Draine 2003 graphite E perp. c
0.1000 = grain radius (micron)
20.00 = grain temperature (K)
......
Downloaded from B. Draine's website at https://www.astro.princeton.edu/~draine/dust/dust.diel.html
ICOMP=11: Draine 2003 astrosilicate
0.1000 = grain radius (micron)
20.00 = grain temperature (K)
......
from cosmic_dustbox import sdist
from cosmic_dustbox import Gspcs as gs
import astropy.units as _u
import os as _os
import numpy as _np
import sys
import pandas
Path = _os.path.dirname(_os.path.dirname(_os.path.realpath(__file__))) + \
'\\lib\\paramsWD01a.csv'
sys.path.insert(0, Path)
class WD01(object):
with open(Path)as f:
R_V = [] ; bc=[] ; alpha_g = [] ; beta_g = [] ; a_tg = [] ; a_cg = []
C_g=[] ; alpha_s = []; beta_s = []; a_ts = [] ; a_cs = [] ; C_s =[]
df = pandas.read_csv(f)
R_V.append(df['R_V'])
bc.append(df['b_C'])
alpha_g.append(df['\\alpha_g'])
beta_g.append(df['\\beta_g'])
a_tg.append(df['a_{t,g} [m]'])
a_cg.append(df['a_{c,g} [m]'])
C_g.append(df['C_g'])
alpha_s.append(df['\\alpha_s'])
beta_s.append(df['\\beta_s'])
a_ts.append(df['a_{t,s} [m]'])
a_cs.append(df['a_{c,s} [m]'])
C_s.append(df['C_s'])
rho_g = 2.24 *_u.g/_u.cm**3
b = _np.array([0.75 * bc[1],0.25 * bc[2]])
#
def __init__(self, Indx):
def load_file(Indx, st):
sdists = {
'gra': sdist.Log_Normal(3.5*_u.Angstrom, _np.inf, rho_g,0.4,bc,b,a0) + \
sdist.WD01_dst(3.5*_u.Angstrom, _np.inf , a_tg[Indx],
beta_g[Indx],a_cg[Indx], alpha_g[Indx],C_g[Indx]) , \
'sil': sdist.WD01(
sdist.WD01_dst(3.5*_u.Angstrom, _np.inf,a_ts[Indx], beta_s[Indx], \
a_cs[Indx],alpha_s[Indx],C_s[Indx])
}
return sdist[st]
self.sdist_gr = load_file(Indx,'gra')
self.sdist_sil = load_file(Indx,'sil')
return
#
\ No newline at end of file
......@@ -171,9 +171,9 @@ class PowerLaw(SizeDist):
return
class DoubleLogNormal(SizeDist):
class LogNormal(SizeDist):
"""
Sum of two log normal distributions as in Eq. (2) of [1]_.
Log normal distribution as in Eq. (2) of [1]_.
References
----------
......@@ -184,28 +184,23 @@ class DoubleLogNormal(SizeDist):
# mass of carbon atom
m_C = 12.0107*_u.u
def B_val(i):
nominator = (3 * _np.exp(-4.5 * sgma**2) * bc[i] * m_C)
denominator = (
(2*_np.pi**2)**1.5 * rho * a0[i]**3 * sgma *
(1 +
_erf((3 * sgma/_np.sqrt(2)) +
(_np.log((a0[i]/3.5/_u.angstrom).decompose().value) /
(sgma * _np.sqrt(2)))
)
)
nominator = (3 * _np.exp(-4.5 * sgma**2) * bc * m_C)
denominator = (
(2*_np.pi)**1.5 * rho * a0**3 * sgma *
(1 +
_erf((3 * sgma/_np.sqrt(2)) +
(_np.log(a0/3.5/_u.angstrom) /
(sgma * _np.sqrt(2)))
)
)
# FIXME: what do we do if the denominator is zero
if denominator != 0:
B = (nominator/denominator).decompose().value
return B
_B = [B_val(i) for i in range(2)]
)
# FIXME: what do we do if the denominator is zero
if denominator != 0:
B = (nominator/denominator).decompose()
def f(a):
return sum([_B[i]/a * _np.exp(-0.5*(
_np.log((a/a0[i]).decompose().value)/sgma)**2)
for i in range(2)])
return B/a * _np.exp(-0.5*(
_np.log(a/a0)/sgma)**2)
super().__init__(sizeMin, sizeMax, f)
return
......@@ -224,19 +219,19 @@ class WD01ExpCutoff(SizeDist):
if beta >= 0:
def F(a):
return 1 + (beta * a) / a_t
return 1.0 + (beta * a) / a_t
else:
def F(a):
return 1 / (1 - (beta * a) / a_t)
return 1.0 / (1 - (beta * a) / a_t)
def exp_func(a):
r = _np.ones_like(a.value)
ind = _np.where(a > a_t)
r[ind] = _np.exp(-(((a[ind] - a_t)/a_c).decompose().value)**3)
r[ind] = _np.exp(-(((a[ind] - a_t)/a_c))**3)
return r
def f(a):
return C*(a / a_t)**alpha/a \
return C/a*(a / a_t)**alpha \
* F(a) * exp_func(a)
super().__init__(sizeMin, sizeMax, f)
......@@ -273,13 +268,21 @@ def WD01(RV, bc, case):
sizeMin = 3.5*_u.angstrom
sizeMax = 10*_u.micron
rho = 2.24*_u.g/_u.cm**3
s_car = DoubleLogNormal(
s_car = LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
0.75*params[2],
3.5*_u.angstrom
) + \
LogNormal(
sizeMin,
sizeMax,
rho,
0.4,
[0.75*params[2], 0.25*params[2]],
[3.5*_u.angstrom, 30*_u.angstrom]
0.25*params[2],
30*_u.angstrom
)
l_car = WD01ExpCutoff(
sizeMin,
......@@ -301,6 +304,7 @@ def WD01(RV, bc, case):
)
return s_car + l_car, sil
###############################################################################
if __name__ == "__main__":
import doctest
......
......@@ -81,6 +81,53 @@ class TestSdist(TestCase):
return
class TestWD01ExpCutoff(TestCase):
def test_simple(self):
sd = sdist.WD01ExpCutoff(
3.5*u.angstrom, 1*u.micron, 0, 0, 1*u.micron, 1*u.micron, 1)
sizes = [
1*u.micron,
1*u.nm,
1*u.angstrom,
2*u.micron,
]
res = [
1./(1*u.micron),
1./(1*u.nm),
0,
0,
]
for s, r in zip(sizes, res):
self.assertEqual(
sd(np.array([s.value])*s.unit)[0], r
)
return
def test_simple_2(self):
sd = sdist.WD01ExpCutoff(
3.5*u.angstrom, 1*u.micron, 0, 0, 1*u.nm, 1*u.micron, 1)
sizes = [
1*u.micron,
2*u.nm,
0.5*u.nm,
1*u.angstrom,
2*u.micron,
]
res = [
1./(1*u.micron) * np.exp(-(1-1e-3)**3),
1./(2*u.nm),
2/u.nm * np.exp(1e-9/8),
u.Quantity(0),
u.Quantity(0),
]
for s, r in zip(sizes, res):
self.assertAlmostEqual(
sd(np.array([s.value])*s.unit)[0].value, r.value
)
return
class TestWD01(TestCase):
def test_smoke_test(self):
......
# Welcome to MkDocs
For full documentation visit [mkdocs.org](https://mkdocs.org).
## Commands
* `mkdocs new [dir-name]` - Create a new project.
* `mkdocs serve` - Start the live-reloading docs server.
* `mkdocs build` - Build the documentation site.
* `mkdocs help` - Print this help message.
## Project layout
mkdocs.yml # The configuration file.
docs/
index.md # The documentation homepage.
... # Other markdown pages, images and other files.
with open('README.md') as fh:
long_desc = fh.read()
name = 'cosmic_dustbox'
version = '0.1.dev0'
authors = 'Martin Glatzle, Ayman Noureldin'
author_email = 'findessp@yandex.ru'
description = 'Toolbox for cosmic dust.'
url = 'https://github.com/findesgh/'+name
license = 'GPL-3.0'
site_name: cosmic_dustbox documentation
site_author: Martin Glatzle
theme: readthedocs
site_dir: docs/html_docs
pages:
- Introduction: 'index.md'
- Quick start: 'quickstart.md'
from setuptools import setup
import metadata as md
with open('README.md') as fh:
long_desc = fh.read()
name = 'cosmic_dustbox'
setup(
name=name,
version='0.1.dev0',
author_email='findessp@yandex.ru',
description='Toolbox for cosmic dust.',
long_description=long_desc,
name=md.name,
version=md.version,
author=md.authors,
author_email=md.author_email,
description=md.description,
long_description=md.long_desc,
long_description_content_type='text/markdown',
url='https://github.com/findesgh/'+name,
license='GPL-3.0',
packages=[name],
url=md.url,
license=md.license,
packages=[md.name],
install_requires=[
'numpy',
'astropy',
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment