-
Notifications
You must be signed in to change notification settings - Fork 68
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #9 from MRCWirtz/master
Add example for confirming Liouville's theorem
- Loading branch information
Showing
5 changed files
with
258 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,99 @@ | ||
|
|
||
| # coding: utf-8 | ||
|
|
||
| # # Galactic Lensing - confirm liouville | ||
| # | ||
| # Galactic lensing can be applied with the requirement that Liouville's theorem | ||
| # holds, thus in this context: from an isotropic distribution outside the area | ||
| # of influence of the Galactic magnetic field follows an isotropic arrival | ||
| # distribution at any point within our Galaxy. | ||
| # First, we are setting up the oberver which we will place further outside of | ||
| # the Galactic center than Earth to exaggerate the observed effects: | ||
|
|
||
| # In[1]: | ||
|
|
||
| import crpropa | ||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
|
|
||
| n = 10000000 | ||
|
|
||
| # simulation setup | ||
| sim = crpropa.ModuleList() | ||
| # we just need propagation in straight lines here to demonstrate the effect | ||
| sim.add(crpropa.SimplePropagation()) | ||
|
|
||
| # collect arriving cosmic rays at observer 19 kpc outside of the Galactic center | ||
| # to exaggerate effects | ||
| obs = crpropa.Observer() | ||
| pos_earth = crpropa.Vector3d(-19, 0, 0) * crpropa.kpc | ||
| # observer with radius 500 pc to collect fast reasonable statistics | ||
| obs.add(crpropa.ObserverSurface(crpropa.Sphere(pos_earth, 0.5 * crpropa.kpc))) | ||
| # Use crpropa's particle collector to only collect the cosmic rays at Earth | ||
| output = crpropa.ParticleCollector() | ||
| obs.onDetection(output) | ||
| sim.add(obs) | ||
|
|
||
| # discard outwards going cosmic rays, that missed Earth and leave the Galaxy | ||
| obs_trash = crpropa.Observer() | ||
| obs_trash.add(crpropa.ObserverSurface(crpropa.Sphere(crpropa.Vector3d(0), 21 * crpropa.kpc))) | ||
| sim.add(obs_trash) | ||
|
|
||
|
|
||
| # ## Lambert's distribution | ||
| # For the source setup we have to consider that from an isotropic propagation | ||
| # in the Extragalactic universe, the directions on any surface element follows | ||
| # the Lambert's distribution (https://en.wikipedia.org/wiki/Lambert%27s_cosine_law). | ||
| # You could also phrase: vertical incident angles are more frequent due to the | ||
| # larger visible size of the area of the surface element than flat angles. | ||
|
|
||
|
|
||
| # In[2]: | ||
|
|
||
| # source setup | ||
| source = crpropa.Source() | ||
| # inward=True for inwards directed emission and False for outwards directed emission | ||
| center, radius, inward = crpropa.Vector3d(0, 0, 0) * crpropa.kpc, 20 * crpropa.kpc, True | ||
| source.add(crpropa.SourceLambertDistributionOnSphere(center, radius, inward)) | ||
| source.add(crpropa.SourceParticleType(-crpropa.nucleusId(1, 1))) | ||
| source.add(crpropa.SourceEnergy(100 * crpropa.EeV)) | ||
|
|
||
| sim.run(source, n) | ||
|
|
||
| print("Number of hits: %i" % len(output)) | ||
| lons = [] | ||
| for c in output: | ||
| v = c.current.getDirection() | ||
| lons.append(v.getPhi()) | ||
|
|
||
| plt.hist(np.array(lons), bins=30, color='k', histtype='step') | ||
| plt.xlabel('lon', fontsize=20) | ||
| plt.ylabel('counts', fontsize=20) | ||
| plt.savefig('lon_distribution_lamberts.png', bbox_inches='tight') | ||
| plt.close() | ||
|
|
||
| # One can see this results in an isotropic arrival distribution. | ||
| # Note, that one instead obtains anisotropies if one assumes an isotropic | ||
| # emission from sources that are distributed uniformly on the sphere shell by e.g. | ||
|
|
||
| # In[3]: | ||
|
|
||
| source = crpropa.Source() | ||
| source.add(crpropa.SourceUniformShell(center, radius)) | ||
| source.add(crpropa.SourceIsotropicEmission()) | ||
| source.add(crpropa.SourceParticleType(-crpropa.nucleusId(1, 1))) | ||
| source.add(crpropa.SourceEnergy(100 * crpropa.EeV)) | ||
|
|
||
| sim.run(source, n) | ||
|
|
||
| print("Number of hits: %i" % len(output)) | ||
| lons = [] | ||
| for c in output: | ||
| v = c.current.getDirection() | ||
| lons.append(v.getPhi()) | ||
|
|
||
| plt.hist(np.array(lons), bins=30, color='k', histtype='step') | ||
| plt.xlabel('lon', fontsize=20) | ||
| plt.ylabel('counts', fontsize=20) | ||
| plt.savefig('lon_distribution_double_isotropic.png', bbox_inches='tight') | ||
| plt.close() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,148 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Galactic Lensing - confirm liouville\n", | ||
| "\n", | ||
| "Galactic lensing can be applied with the requirement that Liouville's theorem\n", | ||
| "holds, thus in this context: from an isotropic distribution outside the area\n", | ||
| "of influence of the Galactic magnetic field follows an isotropic arrival\n", | ||
| "distribution at any point within our Galaxy.\n", | ||
| "First, we are setting up the oberver which we will place further outside of\n", | ||
| "the Galactic center than Earth to exaggerate the observed effects:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import crpropa\n", | ||
| "import matplotlib.pyplot as plt\n", | ||
| "import numpy as np\n", | ||
| "\n", | ||
| "n = 10000000\n", | ||
| "\n", | ||
| "# Simulation setup\n", | ||
| "sim = crpropa.ModuleList()\n", | ||
| "# We just need propagation in straight lines here to demonstrate the effect\n", | ||
| "sim.add(crpropa.SimplePropagation())\n", | ||
| "\n", | ||
| "# collect arriving cosmic rays at Observer 19 kpc outside of the Galactic center\n", | ||
| "# to exaggerate effects\n", | ||
| "obs = crpropa.Observer()\n", | ||
| "pos_earth = crpropa.Vector3d(-19, 0, 0) * crpropa.kpc\n", | ||
| "# observer with radius 500 pc to collect fast reasonable statistics\n", | ||
| "obs.add(crpropa.ObserverSurface(crpropa.Sphere(pos_earth, 0.5 * crpropa.kpc)))\n", | ||
| "# Use crpropa's particle collector to only collect the cosmic rays at Earth\n", | ||
| "output = crpropa.ParticleCollector()\n", | ||
| "obs.onDetection(output)\n", | ||
| "sim.add(obs)\n", | ||
| "\n", | ||
| "# Discard outwards going cosmic rays, that missed Earth and leave the Galaxy\n", | ||
| "obs_trash = crpropa.Observer()\n", | ||
| "obs_trash.add(crpropa.ObserverSurface(crpropa.Sphere(crpropa.Vector3d(0), 21 * crpropa.kpc)))\n", | ||
| "sim.add(obs_trash)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Lambert's distribution\n", | ||
| "For the source setup we have to consider that from an isotropic propagation\n", | ||
| "in the Extragalactic universe, the directions on any surface element follows\n", | ||
| "the Lambert's distribution (https://en.wikipedia.org/wiki/Lambert%27s_cosine_law).\n", | ||
| "You could also phrase: vertical incident angles are more frequent due to the\n", | ||
| "larger visible size of the area of the surface element than flat angles." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# source setup\n", | ||
| "source = crpropa.Source()\n", | ||
| "# inward=True for inwards directed emission, and False for outwards directed emission\n", | ||
| "center, radius, inward = crpropa.Vector3d(0, 0, 0) * crpropa.kpc, 20 * crpropa.kpc, True\n", | ||
| "source.add(crpropa.SourceLambertDistributionOnSphere(center, radius, inward))\n", | ||
| "source.add(crpropa.SourceParticleType(-crpropa.nucleusId(1, 1)))\n", | ||
| "source.add(crpropa.SourceEnergy(100 * crpropa.EeV))\n", | ||
| "\n", | ||
| "sim.run(source, n)\n", | ||
| "\n", | ||
| "print(\"Number of hits: %i\" % len(output))\n", | ||
| "lons = []\n", | ||
| "for c in output:\n", | ||
| " v = c.current.getDirection()\n", | ||
| " lons.append(v.getPhi())\n", | ||
| "\n", | ||
| "plt.hist(np.array(lons), bins=30, color='k', histtype='step')\n", | ||
| "plt.xlabel('lon', fontsize=20)\n", | ||
| "plt.ylabel('counts', fontsize=20)\n", | ||
| "plt.savefig('lon_distribution_lamberts.png', bbox_inches='tight')" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "One can see, this results in an isotropic arrival distribution. Note, that one\n", | ||
| "instead obtain anisotropies if one assumes an isotropic emission from sources\n", | ||
| "that are distributed uniformly on the sphere shell by e.g.." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "source = crpropa.Source()\n", | ||
| "source.add(crpropa.SourceUniformShell(center, radius))\n", | ||
| "source.add(crpropa.SourceIsotropicEmission())\n", | ||
| "source.add(crpropa.SourceParticleType(-crpropa.nucleusId(1, 1)))\n", | ||
| "source.add(crpropa.SourceEnergy(100 * crpropa.EeV))\n", | ||
| "\n", | ||
| "sim.run(source, n)\n", | ||
| "\n", | ||
| "print(\"Number of hits: %i\" % len(output))\n", | ||
| "lons = []\n", | ||
| "for c in output:\n", | ||
| " v = c.current.getDirection()\n", | ||
| " lons.append(v.getPhi())\n", | ||
| "\n", | ||
| "plt.hist(np.array(lons), bins=30, color='k', histtype='step')\n", | ||
| "plt.xlabel('lon', fontsize=20)\n", | ||
| "plt.ylabel('counts', fontsize=20)\n", | ||
| "plt.savefig('lon_distribution_double_isotropic.png', bbox_inches='tight')" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 2 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython2", | ||
| "version": "2.7.16" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 2 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters