Update readme

.gitignore

@@ -2,3 +2,8 @@
 *.pyc
 # emacs backup files
 *~
+# python egg info
+*.egg-info/*
+.cache/*
+.eggs/*
+.pytest_cache/*

README.md

-# cosmic_dustbox
-Toolbox for cosmic dust.
+# A toolbox for cosmic dust.
+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
+Galaxy. Countless observations have since then confirmed this and the
+"interstellar extinction" remains one of the observational pillars of the study
+of cosmic dust. It was later discovered that there is dust also in
+interplanetary space and that it is important for planet formation, there is
+dust in other galaxies, even far away galaxies, and maybe there also is dust in
+intergalactic space. It was confirmed that dust grains do indeed come in
+different sizes and, moreover, in different shapes and different chemical
+compositions. In short, while many things have been understood, it also became
+apparent that cosmic dust is very complicated. The goal of this software
+package is to make work with models of cosmic dust on the computer easier while
+still providing control over all the involved intricacies. An ambitious goal
+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
+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
+having to store information on every single grain. If we investigated the
+properties of single grains, at first the diversity would probably be
+overwhelming since hardly any two grains would be exactly alike. But after
+looking at enough grains, we would start discovering similarities. Grains that
+have similar structure, grains that look like larger or smaller versions of
+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)
+
+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)
+
+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
+simply a grain species combined with a size distribution and we call it a grain
+population. If we think back to the abstract representation of our collection
+of grains, it is not very different from one of the groups. It simply has more
+than one grain species, each with its associated size distribution. So to keep
+it simple, we will also call this a grain population.
+
+So, to recap, what we need to describe a collection of grains are:
+
+- grain species
+- size distributions
+- grain populations
+
+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.
+
+# Literature
+"Physics of the Interstellar and Intergalactic Medium" by B. Draine provides
+a good introduction (and more than that) to cosmic dust.

cosmic_dustbox/models/MRN/mrn.py

+from cosmic_dustbox import sdist
+import astropy.units as _u
+
+sdists = {
+    'gra': sdist.PowerLaw(
+        50*_u.angstrom, 0.25*_u.micron, 10**-25.13*_u.cm**2.5, -3.5),
+    'sil': sdist.PowerLaw(
+        50*_u.angstrom, 0.25*_u.micron, 10**-25.11*_u.cm**2.5, -3.5)
+}

cosmic_dustbox/sdist.py

@@ -41,12 +41,11 @@ class SizeDist(object):
     Examples
     --------
     >>> import astropy.units as u
-
     >>> def f(s): \
             return 1.0/s.unit
     >>> a = SizeDist(3.5*u.angstrom, 1*u.micron, f)
     >>> a(_np.logspace(-11, -5, 10)*u.m)
-    <Quantity [ 0., 0., 0., 1., 1., 1., 1., 1., 0., 0.] 1 / m>
+    <Quantity [0., 0., 0., 1., 1., 1., 1., 1., 0., 0.] 1 / m>
     """
     def __init__(self, sizeMin, sizeMax, func):
         """
@@ -78,6 +77,96 @@ class SizeDist(object):
         r[ind] = self.func(sizes[ind])
         return r
 
+    def __add__(self, other):
+        """
+        Commutative addition of size distributions.
+
+        Overloading of ``+`` operator. Can add an instance of ``SizeDist`` (or
+        subclass thereof), any callable or a scalar to ``self.func``. No check
+        on units is performed. Returns an instance of ``self.__class__``, the
+        ``func`` attribute of which is defined to return the corresponding sum.
+
+        If ``other`` is an instance of ``SizeDist`` (or subclass), take maximum
+        of the two ``sizeMin`` attributes and minimum of the two ``sizeMax``
+        attributes.
+
+        Parameters
+        ----------
+        other : SizeDist, callable or scalar-valued Quantity [1/L]
+            Is added to ``self.func``.
+
+        Returns
+        -------
+         : ``self.__class__``
+            Instance of own class with corresponding ``func`` attribute.
+
+        Examples
+        --------
+        >>> import astropy.units as u
+
+        >>> def f(s): \
+                return 1.0/s.unit
+        >>> a = SizeDist(1*u.angstrom, 1*u.micron, f)
+        >>> b = SizeDist(10*u.angstrom, 10*u.micron, f)
+        >>> c = a + b
+        >>> c(_np.logspace(-11, -4, 10)*u.m)
+        <Quantity [0., 0., 0., 2., 2., 2., 2., 0., 0., 0.] 1 / m>
+        """
+        if issubclass(other.__class__, SizeDist):
+            # find new size limits
+            sizeMin = max(self.sizeMin, other.sizeMin)
+            sizeMax = min(self.sizeMax, other.sizeMax)
+
+            # new differential number density is sum
+            def func(sizes):
+                return self.func(sizes) + other.func(sizes)
+
+            return self.__class__(sizeMin, sizeMax, func)
+        elif callable(other):
+            def func(sizes):
+                return self.func(sizes) + other(sizes)
+            return self.__class__(self.sizeMin, self.sizeMax, func)
+        else:
+            def func(sizes):
+                return other + self.func(sizes)
+            return self.__class__(self.sizeMin, self.sizeMax, func)
+
+    # make addition commutative
+    __radd__ = __add__
+
+    def __mul__(self, other):
+        if callable(other):
+            def func(sizes):
+                return self.function(sizes) * other(sizes)
+            return self.__class__(self.sizeMin, self.sizeMax, func)
+        else:
+            def func(sizes):
+                return other * self.function(sizes)
+            return self.__class__(self.sizeMin, self.sizeMax, func)
+
+    # make multiplication commutative
+    __rmul__ = __mul__
+
+
+class PowerLaw(SizeDist):
+    """
+    Parameters
+    ----------
+    sizeMin : scalar Quantity [L]
+        Low end size cutoff of the distribution.
+    sizeMax : scalar Quantity [L]
+        High end size cutoff of the distribution.
+    power : float
+        Log-slope of the size distribution.
+    C : scalar Quantity [L**(-1-power)]
+        Normalization of the size distribution.
+    """
+    def __init__(self, sizeMin, sizeMax, power, C):
+        def f(a):
+            return C*a**power
+        super(self).__init__(sizeMin, sizeMax, f)
+        return
+
 
 ###############################################################################
 if __name__ == "__main__":

cosmic_dustbox/tests/test_sdist.py

+from unittest import TestCase
+from .. import sdist
+import numpy as np
+import astropy.units as u
+
+
+class TestSdist(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        cls._sizes = np.logspace(-11, -5)*u.m
+        cls._unitLessSizes = cls._sizes.value
+        return
+
+    @classmethod
+    def tearDownClass(cls):
+        return
+
+    def test_init(self):
+        """
+        Make sure we can create an instance.
+        """
+        sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1/x.unit)
+
+    def test_call(self):
+        low = 4.0*u.angstrom
+        high = 1*u.micron
+        a = sdist.SizeDist(low, high, lambda x: 1/x.unit)
+        r = a(self.__class__._sizes)
+        self.assertTrue(
+            np.all(r[self.__class__._sizes > high] == 0.0))
+        self.assertTrue(
+            np.all(r[self.__class__._sizes < low] == 0.0))
+        self.assertTrue(
+            np.all(r[np.where(r > 0)] == 1.0/u.m))
+        return
+
+    def test_unitless_sizes(self):
+        a = sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1/x.unit)
+        with self.assertRaises(AttributeError) as context:
+            a(self.__class__._unitLessSizes)
+        return
+
+    def test_unitless_func(self):
+        a = sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1)
+        with self.assertRaises(u.core.UnitConversionError) as context:
+            a(self.__class__._sizes)
+        return
+
+    def test_scalar_sizes(self):
+        """
+        This is something to be fixed in the future.
+        """
+        a = sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1/x.unit)
+        with self.assertRaises(TypeError) as context:
+            a(0.5*u.micron)
+        return
+
+    def test_add_scalar_quantity_int(self):
+        scalar = 1/u.micron
+        a = sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1/x.unit)
+        b = a + scalar
+        r = a(self.__class__._sizes)
+        r[np.where(r != 0)] = r[np.where(r != 0)] + scalar
+        self.assertTrue(
+            np.all(b(self.__class__._sizes) == r))
+        return
+
+    def test_add_scalar_quantity_float(self):
+        scalar = 2.5/u.micron
+        a = sdist.SizeDist(3.5*u.angstrom, 1*u.micron, lambda x: 1/x.unit)
+        b = a + scalar
+        r = a(self.__class__._sizes)
+        r[np.where(r != 0)] = r[np.where(r != 0)] + scalar
+        self.assertTrue(
+            np.all(b(self.__class__._sizes) == r))
+        return

docs/index.md

+# 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.

figs/grain-pop.tex

+\documentclass[crop, tikz]{standalone}
+
+\usetikzlibrary{math}
+\begin{document}
+\begin{tikzpicture}
+  \draw[thick] (0, 0) rectangle (10, 10);
+
+  % red circles
+  % import numpy as np
+  % for x, y, r in zip(10*(np.random.random(5)), 10*(np.random.random(5)), 0.5*(np.random.random(5))):
+  %     print(str(x)+'/'+str(y)+'/'+str(r)+',')
+
+  \foreach \xcord / \ycord / \rcord in {
+    7.679061860689958/3.682979874355471/0.4565838404947177,
+    9.79222926473049/5.283254108481699/0.014712651759808348,
+    9.228163854796406/0.465467574267856/0.12417262955376379,
+    8.466903599867722/2.6335052370495404/0.15078170372073274,
+    3.8702372331492265/5.903399027966185/0.23838775192054612,
+    4.986562338321869/1.7117275570073676/0.044592439438311426,
+    7.448863277700564/8.226461665494632/0.06244327990940718,
+    6.306622283446704/0.5491178028038102/0.3074134849629631,
+    6.831041908490814/8.333630001447492/0.05653784106297699
+  }{
+    \draw[red, thick] (\xcord, \ycord) circle (\rcord);
+  }
+
+  % blue circles
+  \foreach \x / \y / \r in {
+    1.2293473513913522/4.117458330953575/0.8282351326975819,
+    7.687835189158383/6.525741647740903/0.4959151367782896,
+    5.07067407779401/4.413937901291682/0.32343724210759717,
+    6.252736357650197/8.899430585033288/0.11194328613448434
+  }{
+    \draw[blue, thick] (\x, \y) circle (\r);
+  }
+
+  % green ellipses
+  \foreach \x / \y / \a / \b / \ang in {
+    5.140384055187862/5.996278145733781/0.44458636535075624/0.13337590960522686/42.8959206218665,
+    3.234418862988054/4.945278056645622/0.36019401697947695/0.10805820509384308/132.77692008629268,
+    0.6563379428642213/6.504964785235348/0.19991550807517933/0.059974652422553794/124.22338658231189,
+    9.197961905303377/2.205439358917957/0.27526797473586584/0.08258039242075975/100.8813649933555,
+    0.3062812177004437/1.4677046604730348/0.42017352508362876/0.12605205752508863/68.71741736618897
+  }{
+    \draw[green, thick, rotate around={\ang:(\x,\y)}] (\x, \y) ellipse (\a cm and  \b cm);
+  }
+
+  % blue triangles
+  % import numpy as np
+  % for x, y, a, ang in zip(10*(np.random.random(5)), 10*(np.random.random(5)), 0.5*(np.random.random(5)), 180*(np.random.random(5))):
+  %    print(str(x)+'/'+str(y)+'/'+str(a)+'/'+str(0.3*a)+'/'+str(ang)+',')
+
+  \foreach \x / \y / \a / \b / \ang in {
+    3.1825266671917642/8.739378860463935/0.3076868629010771/0.09230605887032313/55.44224580874666,
+    5.9953999833826535/1.1586268023457924/0.31485139062877743/0.09445541718863322/71.98786362988174,
+    9.810321336977507/9.65425711309965/0.09875441625998671/0.02962632487799601/76.7939584508893,
+    9.05699195598717/5.39917055694529/0.419809206167695/0.1259427618503085/22.929429244839067,
+    4.835379117540814/2.649414224653106/0.13385616813162254/0.04015685043948676/13.633175956536425,
+    6.858701607071196/0.6665273841406083/0.19102094661889218/0.05730628398566765/142.1722416667769,
+    0.25352036670998834/9.508868483290078/0.10883725448703785/0.03265117634611135/87.79022903558271
+  }{
+    \draw[blue, thick, rotate around={\ang:(\x,\y)}] (\x, \y) \foreach \s in {120, 240}{
+      -- ++ (\s:\a)
+    } -- cycle ++ (90:\a);
+  }
+\end{tikzpicture}
+\end{document}

figs/grain-species.tex

+\documentclass[crop, tikz]{standalone}
+
+\usetikzlibrary{math}
+\usetikzlibrary{calc}
+\usetikzlibrary{positioning}
+\begin{document}
+\begin{tikzpicture}
+  \draw[thick] (0, 0) rectangle (5, 5);
+  \draw[thick] (5, 0) rectangle (10, 5);
+  \draw[thick] (0, 5) rectangle (5, 10);
+  \draw[thick] (5, 5) rectangle (10, 10);
+
+  % red circles
+  % import numpy as np
+  % for x, y, r in zip(10*(np.random.random(5)), 10*(np.random.random(5)), 0.5*(np.random.random(5))):
+  %     print(str(x)+'/'+str(y)+'/'+str(r)+',')
+
+  \foreach \xcord / \ycord / \rcord in {
+    {1*5/4} / {3*5/4} / 0.4565838404947177,
+    {2*5/4} / {3*5/4} / 0.3074134849629631,
+    {3*5/4} / {3*5/4} / 0.23838775192054612,
+    {1*5/4} / {2*5/4} / 0.15078170372073274,
+    {2*5/4} / {2*5/4} / 0.12417262955376379,
+    {3*5/4} / {2*5/4} / 0.06244327990940718,
+    {1*5/4} / {1*5/4} / 0.044592439438311426,
+    {2*5/4} / {1*5/4} / 0.05653784106297699,
+    {3*5/4} / {1*5/4} / 0.014712651759808348
+  }{
+    \draw[red, thick] (\xcord, \ycord) circle (\rcord);
+  }
+
+  % blue circles
+  \foreach \x / \y / \r in {
+    {5+1*5/3} / {2*5/3} / 0.8282351326975819,
+    {5+2*5/3} / {2*5/3} / 0.4959151367782896,
+    {5+1*5/3} / {1*5/3} / 0.32343724210759717,
+    {5+2*5/3} / {1*5/3} / 0.11194328613448434
+  }{
+    \draw[blue, thick] (\x,\y) circle (\r);
+  }
+
+  % green ellipses
+  \foreach \x / \y / \a / \b in {
+    {1*5/4} / {5+2*5/3} / 0.44458636535075624 / 0.13337590960522686  ,
+    {2*5/4} / {5+2*5/3} / 0.42017352508362876 / 0.12605205752508863  ,
+    {3*5/4} / {5+2*5/3} / 0.36019401697947695 / 0.10805820509384308  ,
+    {1*5/3} / {5+1*5/3} / 0.27526797473586584 / 0.08258039242075975  ,
+    {2*5/3} / {5+1*5/3} / 0.19991550807517933 / 0.059974652422553794
+  }{
+    \draw[green, thick] (\x, \y) ellipse (\b cm and  \a cm);
+  }
+
+  % blue triangles
+  % import numpy as np
+  % for x, y, a, ang in zip(10*(np.random.random(5)), 10*(np.random.random(5)), 0.5*(np.random.random(5)), 180*(np.random.random(5))):
+  %    print(str(x)+'/'+str(y)+'/'+str(a)+'/'+str(0.3*a)+'/'+str(ang)+',')
+
+  % This computation of the x coords is ugly as all hell. Need to right
+  % shift the triangles by \a/2 since they're anchored at the bottom right
+  % corner. Couldn't get it to work in the foreach loop.
+  \foreach \x / \y / \a in {
+    {5+1*5/4+0.419809206167695/2}   / {5+3*5/4} / 0.419809206167695  ,
+    {5+2*5/4+0.31485139062877743/2} / {5+3*5/4} / 0.31485139062877743,
+    {5+3*5/4+0.3076868629010771/2}  / {5+3*5/4} / 0.3076868629010771 ,
+    {5+1*5/4+0.19102094661889218/2} / {5+2*5/4} / 0.19102094661889218,
+    {5+2*5/4+0.13385616813162254/2} / {5+2*5/4} / 0.13385616813162254,
+    {5+3*5/4+0.10883725448703785/2} / {5+2*5/4} / 0.10883725448703785,
+    {5+2*5/4+0.09875441625998671/2} / {5+1*5/4} / 0.09875441625998671
+  }{
+    \draw[blue, thick] (\x,\y) \foreach \s in {120, 240}{
+      -- ++ (\s:\a)
+    } -- cycle ++ (90:\a);
+  }
+\end{tikzpicture}
+\end{document}

mkdocs.yml

+site_name: cosmic_dustbox documentation
+site_author: Martin Glatzle
+theme: readthedocs
+site_dir: docs/html_docs
+pages:
+- Introduction: 'index.md'
+- Quick start: 'quickstart.md'

setup.py

+from setuptools import setup
+
+with open('README.md') as fh:
+    long_desc = fh.read()
+
+setup(
+    name='cosmic_dustbox',
+    version='0.1.dev',
+    author_email='findessp@yandex.ru',
+    description='Toolbox for cosmic dust.',
+    long_description=long_desc,
+    long_description_content_type='text/markdown',
+    url='https://github.com/findesgh/cosmic_dustbox',
+    license='GPL-3.0',
+    packages=['cosmic_dustbox'],
+    install_requires=[
+        'numpy',
+        'astropy'
+    ],
+    tests_require=[
+        'pytest'
+    ],
+    zip_safe=False,
+)