Commit adf7b220 authored by Francesco Montanari's avatar Francesco Montanari

Merge branch 'devel'

* devel:
  Remove cosmo-pedia
  Test comoving volume
  Add comoving-volume command
  Add parallax distance and comoving volume
parents 47370502 48eef127
......@@ -48,7 +48,6 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
expansion of the Universe.
Definitions follow Hogg (1999) [[[https://arxiv.org/abs/astro-ph/9905116][arXiv:astro-ph/9905116]]]. The
=cosmo-pedia= command is based on [[http://fiteoweb.unige.ch/~durrer/Book.html][Durrer (2008)]]. The
=cosmo-calculator= command is based on Gnuastro =astcosmiccal=
program (v0.2).
......@@ -116,6 +115,9 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
<...>
#+END_EXAMPLE
- cosmo-comoving-volume :: Display comoving volume in
mini-buffer.
- cosmo-hubble :: Display Hubble parameter in mini-buffer.
- cosmo-hubble-distance :: Display Hubble distance c/H0 [Mpc]
......@@ -130,8 +132,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
- cosmo-luminosity-distance :: Display luminosity distance in
mini-buffer.
- cosmo-pedia :: Display a reference to basic cosmological
definitions.
- cosmo-parallax-distance :: Display parallax distance in mini-buffer.
- cosmo-set-params :: Change the values of cosmological parameters.
......
......@@ -3,9 +3,9 @@
;; Copyright (C) 2017 Francesco Montanari
;;
;; Author: Francesco Montanari <fmnt@fmnt.info>
;; Maintainer: Francesco Montanari <fmnt@fmnt.info>
;; Created: 22 April 2017
;; Version: 0.1
;; Package-Requires: ((emacs "24.4"))
;; Keywords: tools
;; Homepage: https://gitlab.com/montanari/cosmo-el
;;
......@@ -60,18 +60,9 @@
;;; Todo:
;; In priority order:
;;
;; - Refactor tests, now they are too cumbersome and repetitive.
;;
;; - Add all quantities from Hogg 1999.
;;
;; - Suggest default parameters when reading them with the related
;; command; set the to default values if none is entered.
;;
;; - Allow users to load a customized file with the cosmo-pedia
;; command (or describe how to do it from ~/.emacs).
;;
;; - (Consider the following only if performance becomes critical.) At
;; a fixed redshift, only cosmo-get-los-comoving-distance perform
;; the integral. Other distances are defined in terms of this
......@@ -95,42 +86,6 @@
"Maximum number of steps in the numerical integral algorithm."
:group 'cosmo)
(defvar cosmo-pedia-text
;; Text appearing in the *Cosmopedia* buffer.
"(`q` to quite)
Units system: hbar = c = k_Boltzmann = 1.
Distances relations
-------------------
- Comoving distance (transverse): D_M
- Angular diameter distance: D_A = D_M / (1+z)
- Luminosity distance: D_L = (1+z) D_M = (1+z)^2 D_A
Conversion factors, units
-------------------------
- 1GeV = 1.6022e-3 erg
= 1.1605e13 K
= 1.7827e-24 g
= 5.0684e13 1/cm
= 1.5192 1/s
- 1 pc = 3.2612 light years
= 3.0856e18 cm
- 1 Mpc = 1e6pc ~ 3e24 cm ~ 1e14 s
- 1 AU = 1.4960e13 cm
- 1 Jy = 1e-23 erg/cm^2/s/Hz
= 2.4730e-48 GeV^3
- 1 yr ~ pi * 1e17 s
Important constants
-------------------
- Hubble constant: H0 = 100h km/s/Mpc
= 2.1332e-42 h GeV
- Hubble time, distance: 1/H0 = 3.0856e17/h s
= 9.7776e9/h yr
= 2997.9/h Mpc
= 9.2503e27/h cm")
;; Hash table containing all independent cosmological parameters.
(defvar cosmo--params
(let ((table (make-hash-table :test #'equal)))
......@@ -182,6 +137,10 @@ Important constants
"Hyperbolic sine of real arguments X."
(* 0.5 (- (exp x) (exp (- x)))))
(defun cosmo-asinh (x)
"Inverse hyperbolic sine of real arguments X."
(log (+ x (sqrt (1+ (* x x))))))
(defun cosmo--trapzd (func lo hi sum niter)
"Extended trapezoidal rule to integrate FUNC from LO to HI.
Compute the refinement to an old SUM result, previous to the
......@@ -350,6 +309,55 @@ Optional argument JMAX maximum number of steps."
(message (format "%s Mpc"
(cosmo-get-luminosity-distance z)))))
(defun cosmo-get-parallax-distance (redshift)
"Parallax distance [Mpc] for Lambda-CDM at a given REDSHIFT."
(let* ((DH (cosmo-get-hubble-distance))
(DM (cosmo-get-transverse-comoving-distance redshift))
(dM (/ DM DH))
(ocurvature (cosmo-get-ocurvature)))
(/ DM (+ dM (sqrt (1+ (* ocurvature dM dM)))))))
(defun cosmo-parallax-distance ()
"Display parallax distance in mini-buffer."
(interactive)
(let ((z (cosmo--read-param "redshift")))
(message (format "%s Mpc"
(cosmo-get-parallax-distance z)))))
(defun cosmo--get-comoving-volume-nonflat (DM DH ocurvature)
"Return the comoving volume for non-vanishing curvature.
Argument DM comoving distance (transverse).
Argument DH Hubble distance.
Argument OCURVATURE curvature density parameter."
(let* ((DM-over-DH (/ DM DH))
(sqrt-ok (sqrt (abs (cosmo-get-ocurvature))))
(pref (* 2.0 float-pi (/ (expt DH 3.0) ocurvature)))
(func (cond ((> ocurvature 0.0)
#'cosmo-asinh)
((< ocurvature 0.0)
#'asin)
(t (error "Error: wrong curvature parameter value")))))
(setq term1 (* DM-over-DH
(sqrt (1+ (* ocurvature (expt DM-over-DH 2.0))))))
(setq term2 (/ (funcall func (* sqrt-ok DM-over-DH)) sqrt-ok))
(* pref (- term1 term2))))
(defun cosmo-get-comoving-volume (redshift)
"Comoving volume [Mpc^3] for Lambda-CDM at a given REDSHIFT."
(let* ((DH (cosmo-get-hubble-distance))
(DM (cosmo-get-transverse-comoving-distance redshift))
(ocurvature (cosmo-get-ocurvature)))
(if (= ocurvature 0.0)
(* (/ 4.0 3.0) float-pi (expt DM 3.0))
(cosmo--get-comoving-volume-nonflat DM DH ocurvature))))
(defun cosmo-comoving-volume ()
"Display comoving volume in mini-buffer."
(interactive)
(let ((z (cosmo--read-param "redshift")))
(message (format "%s Mpc^3"
(cosmo-get-comoving-volume z)))))
;;; Handle output.
(defun cosmo--write-calc-header ()
......@@ -433,14 +441,6 @@ Argument ANGULAR-DIST angular diameter distance at given redshift."
los-dist transverse-dist luminosity-dist
angular-dist))))
(defun cosmo-pedia ()
"Display a reference to basic cosmological definitions."
(interactive)
(let* ((cosmo-buffer "*Cosmopedia*"))
(with-output-to-temp-buffer cosmo-buffer
(pop-to-buffer cosmo-buffer)
(insert cosmo-pedia-text))))
(provide 'cosmo)
;;; cosmo.el ends here
......@@ -50,6 +50,8 @@ def print_curves(reds, cosmo):
name='luminosity-distance [Mpc]')
sexp_print(reds, cosmo.angular_diameter_distance(reds).value,
name='angular-diameter-distance [Mpc]')
sexp_print(reds, cosmo.comoving_volume(reds).value,
name='comoving-volume [Mpc^3]')
def print_open_cosmo():
......@@ -59,7 +61,9 @@ def print_open_cosmo():
Ode0 = 0.7
cosmo = LambdaCDM(H0=H0, Om0=Om0, Ode0=Ode0)
print('* Open cosmology:')
print('===============')
print('Open cosmology:')
print('===============')
print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
.format(H0, Om0, Ode0, cosmo.Ogamma(0.)+cosmo.Onu(0.),
cosmo.Ok(0.)))
......@@ -75,7 +79,9 @@ def print_close_cosmo():
Ode0 = 0.7
cosmo = LambdaCDM(H0=H0, Om0=Om0, Ode0=Ode0)
print('* Close cosmology:')
print('================')
print('Close cosmology:')
print('================')
print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
.format(H0, Om0, Ode0, cosmo.Ogamma(0.)+cosmo.Onu(0.),
cosmo.Ok(0.)))
......@@ -90,7 +96,9 @@ def print_flat_cosmo():
Om0 = 0.3
cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Tcmb0=0.)
print('* Flat cosmology:')
print('===============')
print('Flat cosmology:')
print('===============')
print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
.format(H0, Om0, cosmo.Ode(0.), cosmo.Ogamma(0.)+cosmo.Onu(0.),
cosmo.Ok(0.)))
......
......@@ -28,6 +28,10 @@
;; This package provides an unit test for cosmo.el, a cosmological
;; calculator.
;;; Todo:
;; - Refactor tests, now they are too cumbersome and repetitive.
;;; Code:
......@@ -119,6 +123,10 @@ Argument OREL relativistic density parameter today."
(assert (cosmo-almost-eq
(cosmo-get-angular-diameter-distance 1000.0) 13.2194016018 1e-3)))
(defun cosmo-test-comoving-volume ()
(assert (cosmo-almost-eq
(cosmo-get-comoving-volume 1000.0) 1.00014515316e+13 3e-3)))
(cosmo-test-efunc)
(cosmo-test-inv-efunc)
(cosmo-test-hubble)
......@@ -128,6 +136,7 @@ Argument OREL relativistic density parameter today."
(cosmo-test-transverse-comoving-distance-open)
(cosmo-test-luminosity-distance)
(cosmo-test-angular-diameter-distance)
(cosmo-test-comoving-volume)
;;; Close cosmology.
......@@ -144,7 +153,12 @@ Argument OREL relativistic density parameter today."
(cosmo-get-transverse-comoving-distance 1000.0)
14832.933576 1e-3)))
(defun cosmo-test-comoving-volume ()
(assert (cosmo-almost-eq
(cosmo-get-comoving-volume 1000.0) 1.24288502093e+13 3e-3)))
(cosmo-test-transverse-comoving-distance-close)
(cosmo-test-comoving-volume)
;;; Flat cosmology.
......@@ -158,10 +172,15 @@ Argument OREL relativistic density parameter today."
(defun cosmo-test-transverse-comoving-distance-flat ()
(assert (cosmo-almost-eq
(cosmo-get-transverse-comoving-distance 1000.0)
13660.5292969 1e-3)))
(cosmo-get-transverse-comoving-distance 1000.0) 13660.5292969 1e-3)))
(defun cosmo-test-comoving-volume ()
(assert (cosmo-almost-eq
(cosmo-get-comoving-volume 1000.0) 1.06780313213e+13 3e-3)))
(cosmo-test-transverse-comoving-distance-flat)
(cosmo-test-comoving-volume)
;;; Benchmarks
......
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