Re-write Galaxy modeling code

docs/astro/index.rst

@@ -21,6 +21,10 @@ The ``gammapy.astro`` namespace is empty ... use these import statements::
   source/index
   population/index
 
+Getting Started
+===============
+
+
 
 Example
 =======

docs/astro/population/index.rst

@@ -9,14 +9,75 @@ Astrophysical source population models (`gammapy.astro.population`)
 Introduction
 ============
 
-`gammapy.astro.population` implements some common astrophysical source and population models.
-
-TODO: describe
+The `gammapy.astro.population` module provides a simple framework for population synthesis of 
+gamma-ray sources. 
 
 Getting Started
 ===============
 
-TODO: describe
+TODO
+
+
+Radial surface density distributions
+====================================
+Here is a comparison plot of all available radial distribution functions of the surface density of pulsars
+and related objects used in literature:
+
+.. plot::
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    from gammapy.astro.population import radial_distributions
+    from gammapy.utils.distributions import normalize
+
+    max_radius = 20  # kpc
+    r = np.linspace(0, max_radius, 100)
+    colors = ['b', 'k', 'k', 'b', 'g', 'g']
+
+    for color, key in zip(colors, radial_distributions.keys()):
+        model = radial_distributions[key]()
+        if model.evolved:
+            linestyle = '-'
+        else:
+            linestyle = '--'
+        label = model.__class__.__name__
+        plt.plot(r, normalize(model, 0, max_radius)(r), color=color, linestyle=linestyle, label=label)
+    plt.xlim(0, max_radius)
+    plt.ylim(0, 0.28)
+    plt.xlabel('Galactocentric Distance [kpc]')
+    plt.ylabel('Normalized Surface Density [kpc^-2]')
+    plt.legend(prop={'size': 10})
+    plt.show()
+
+
+Velocity distributions
+======================
+Here is a comparison plot of all available velocity distribution functions:
+
+.. plot::
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    from gammapy.astro.population import velocity_distributions
+    from gammapy.utils.distributions import normalize
+
+    v_min, v_max = 10, 3000  # km / s
+    v = np.linspace(v_min, v_max, 200)
+    colors = ['b', 'k', 'g']
+
+    for color, key in zip(colors, velocity_distributions.keys()):
+    	model = velocity_distributions[key]()
+    	label = model.__class__.__name__
+    	plt.plot(v, normalize(model, v_min, v_max)(v), color=color, linestyle='-', label=label)
+
+    plt.xlim(v_min, v_max)
+    plt.ylim(0, 0.004)
+    plt.xlabel('Velocity [km/s]')
+    plt.ylabel('Probability Density [(km / s)^-1]')
+    plt.semilogx()
+    plt.legend(prop={'size': 10})
+    plt.show()
+
 
 Reference/API
 =============

docs/astro/source/index.rst

@@ -9,14 +9,22 @@ Astrophysical source models (`gammapy.astro.source`)
 Introduction
 ============
 
-`gammapy.astro.source` implements some common astrophysical source and population models.
+The `gammapy.astro.source` module contains classes of source models, which can be used for population synthesis of galactic gamma-ray sources.
+
+
+Source Models
+=============
+The following source models are available:
+
+.. toctree::
+   :maxdepth: 1
+
+   snr
+   pwn
+   pulsar
 
-TODO: describe
 
-Getting Started
-===============
 
-TODO: describe
 
 Reference/API
 =============

docs/astro/source/pulsar.rst

@@ -0,0 +1,23 @@
+Pulsar Source Models
+====================
+
+Plot spin frequency of the pulsar with time:
+
+
+ .. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.units import Quantity
+    from gammapy.astro.source import Pulsar
+
+    t = Quantity(np.logspace(0, 6, 100), 'yr')
+    pulsar = Pulsar(P_0=Quantity(0.01, 's'), logB=12)
+    plt.plot(t.value, 1 / pulsar.period(t).cgs.value, color='b')
+    plt.xlabel('time [years]')
+    plt.ylabel('frequency [1/s]')
+    plt.ylim(1E0, 1E2)
+    plt.loglog()
+    plt.show()
+

docs/astro/source/pwn.rst

@@ -0,0 +1,28 @@
+Pulsar Wind Nebula Source Models
+================================
+
+Plot the evolution of the radius of the PWN:
+
+
+ .. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.units import Quantity
+    from astropy.constants import M_sun
+    from gammapy.astro.source import PWN, SNRTrueloveMcKee
+
+    t = Quantity(np.logspace(1, 5, 100), 'yr')
+    n_ISM = Quantity(1, 'cm^-3')
+    snr = SNRTrueloveMcKee(m_ejecta=8 * M_sun, n_ISM=n_ISM)
+    pwn = PWN(snr=snr)
+    pwn.pulsar.L_0 = Quantity(1E40, 'erg/s')
+    plt.plot(t.value, pwn.radius(t).to('pc').value, label='Radius SNR')
+    plt.plot(t.value, pwn.snr.radius_reverse_shock(t).to('pc').value, label='Reverse Shock SNR')
+    plt.plot(t.value, pwn.snr.radius(t).to('pc').value, label='Radius PWN')
+    plt.xlabel('time [years]')
+    plt.ylabel('radius [pc]')#
+    plt.legend(loc=4)
+    plt.loglog()
+    plt.show()

docs/astro/source/snr.rst

@@ -0,0 +1,61 @@
+Supernova Remnant Models
+========================
+
+Plot the evolution of radius of the SNR:
+
+
+ .. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.units import Quantity
+    from gammapy.astro.source import SNR, SNRTrueloveMcKee
+
+    snr_models = [SNR, SNRTrueloveMcKee]
+    densities = Quantity([1, 0.1], 'cm^-3')
+    linestyles = ['-', '--']
+    colors = ['b', 'g']
+    t = Quantity(np.logspace(0, 5, 100), 'yr')
+
+    for color, density in zip(colors, densities):
+        for linestyle, snr_model in zip(linestyles, snr_models):
+            snr = snr_model(n_ISM=density)
+            plt.plot(t.value, snr.radius(t).to('pc').value,
+                     label=snr.__class__.__name__ + ' (n_ISM = {0})'.format(density.value),
+                     linestyle=linestyle, color=color)
+    plt.xlabel('time [years]')
+    plt.ylabel('radius [pc]')
+    plt.legend(loc=4)
+    plt.loglog()
+    plt.show()
+
+
+Plot the evolution of the flux of the SNR above 1 TeV and at 1 kpc distance:
+
+ .. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.units import Quantity
+    from gammapy.astro.source import SNR, SNRTrueloveMcKee
+
+    densities = Quantity([1, 0.1], 'cm^-3')
+    colors = ['b', 'g']
+
+    t = Quantity(np.logspace(0, 5, 100), 'yr')
+
+    for color, density in zip(colors, densities):
+        snr = SNR(n_ISM=density)
+        F = snr.luminosity_tev(t) / (4 * np.pi * Quantity(1, 'kpc') ** 2)
+        plt.plot(t.value, F.to('ph s^-1 cm^-2').value, color=color, label='n_ISM = {0}'.format(density.value))
+        plt.vlines(snr.sedov_taylor_begin.to('yr').value, 1E-13, 1E-10, linestyle='--', color=color)
+        plt.vlines(snr.sedov_taylor_end.to('yr').value, 1E-13, 1E-10, linestyle='--', color=color)
+    plt.xlim(1E2, 1E5)
+    plt.ylim(1E-13, 1E-10)
+    plt.xlabel('time [years]')
+    plt.ylabel('flux @ 1kpc [ph s^-1 cm ^-2]')
+    plt.legend(loc=4)
+    plt.loglog()
+    plt.show()