Merge d7921b4 into 3247019

README.md

@@ -11,7 +11,7 @@ Installation
 To install the package:
 
 ```jlcon
-julia> Pkg.add("Cosmology")
+pkg> add Cosmology
 ```
 
 Then, to load into your session:
@@ -20,6 +20,9 @@ Then, to load into your session:
 julia> using Cosmology
 ```
 
+`Cosmology.jl` uses [`Unitful.jl`](https://github.com/ajkeller34/Unitful.jl) and
+[`UnitfulAstro.jl`](https://github.com/JuliaAstro/UnitfulAstro.jl) to handle units.
+
 Cosmological Models
 -------------------
 
@@ -79,19 +82,19 @@ Distances
 
 <table>
   <tr>
-    <td>angular_diameter_dist_mpc(cosmo,&nbsp;z)</td>
+    <td>angular_diameter_dist(cosmo,&nbsp;z)</td>
     <td>Ratio of an object's proper transverse size (in Mpc) to its angular size (in radians)</td>
   </tr>
   <tr>
-    <td>comoving_radial_dist_mpc(cosmo,&nbsp;z)</td>
+    <td>comoving_radial_dist(cosmo,&nbsp;z)</td>
     <td>Comoving radial distance to redshift z, in Mpc</td>
   </tr>
   <tr>
-    <td>comoving_volume_gpc3(cosmo,&nbsp;z)</td>
+    <td>comoving_volume(cosmo,&nbsp;z)</td>
     <td>Comoving volume out to redshift z, in Gpc<sup>3</sup></td>
   </tr>
   <tr>
-    <td>luminosity_dist_mpc(cosmo, z)</td>
+    <td>luminosity_dist(cosmo, z)</td>
     <td>Bolometric luminosity distance, in Mpc</td>
   </tr>
   <tr>
@@ -106,20 +109,20 @@ julia> using Cosmology
 julia> c = cosmology(OmegaM=0.26)
 FlatLCDM(0.69,0.7399122024007928,0.26,8.779759920715362e-5)
 
-julia> angular_diameter_dist_mpc(c, 1.2)
-1784.0089227105118
+julia> angular_diameter_dist(c, 1.2)
+1784.0089227105113 Mpc
 ```
 
 Times
 -----
 
 <table>
   <tr>
-    <td>age_gyr(cosmo, z)</td>
+    <td>age(cosmo, z)</td>
     <td>Age of the universe at redshift z, in Gyr</td>
   </tr>
   <tr>
-    <td>lookback_time_gyr(cosmo, z)</td>
+    <td>lookback_time(cosmo, z)</td>
     <td>Difference between age at redshift 0 and age at redshift z, in Gyr</td>
   </tr>
 </table>
@@ -130,6 +133,6 @@ julia> using Cosmology
 julia> c = cosmology(OmegaM=0.26)
 FlatLCDM(0.69,0.7399122024007928,0.26,8.779759920715362e-5)
 
-julia> age_gyr(c, 1.2)
-5.445600787626434
+julia> age(c, 1.2)
+5.445600787626434 Gyr
 ```

REQUIRE

@@ -1,2 +1,4 @@
 julia 0.7
 QuadGK 0.1.1
+Unitful 0.9.0
+UnitfulAstro 0.1.1

src/Cosmology.jl

@@ -3,19 +3,21 @@ __precompile__()
 module Cosmology
 
 using QuadGK
+using Unitful: km, s
+using UnitfulAstro: Mpc, Gyr
 
 export cosmology,
-       age_gyr,
-       angular_diameter_dist_mpc,
-       comoving_radial_dist_mpc,
-       comoving_transverse_dist_mpc,
-       comoving_volume_gpc3,
+       age,
+       angular_diameter_dist,
+       comoving_radial_dist,
+       comoving_transverse_dist,
+       comoving_volume,
        distmod,
        H,
-       hubble_dist_mpc,
-       hubble_time_gyr,
-       luminosity_dist_mpc,
-       lookback_time_gyr,
+       hubble_dist,
+       hubble_time,
+       luminosity_dist,
+       lookback_time,
        scale_factor
 
 abstract type AbstractCosmology end
@@ -133,65 +135,65 @@ end
 
 scale_factor(z) = 1/(1 + z)
 E(c::AbstractCosmology, z) = (a = scale_factor(z); a2E(c,a)/a^2)
-H(c::AbstractCosmology, z) = 100. * c.h * E(c, z)
+H(c::AbstractCosmology, z) = 100. * c.h * E(c, z) * km / s / Mpc
 
-hubble_dist_mpc0(c::AbstractCosmology) = 2997.92458/c.h
-hubble_dist_mpc(c::AbstractCosmology, z) = hubble_dist_mpc0(c)/E(c,z)
+hubble_dist0(c::AbstractCosmology) = 2997.92458/c.h * Mpc
+hubble_dist(c::AbstractCosmology, z) = hubble_dist0(c)/E(c,z)
 
-hubble_time_gyr0(c::AbstractCosmology) = 9.77814/c.h
-hubble_time_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)/E(c,z)
+hubble_time0(c::AbstractCosmology) = 9.77814/c.h * Gyr
+hubble_time(c::AbstractCosmology, z) = hubble_time0(c)/E(c,z)
 
 # distances
 
 Z(c::AbstractCosmology, z::Real; kws...) =
     QuadGK.quadgk(a::Float64->1.0/a2E(c,a), scale_factor(z), 1; kws...)[1]
 
-comoving_radial_dist_mpc(c::AbstractCosmology, z; kws...) = hubble_dist_mpc0(c)*Z(c, z; kws...)
+comoving_radial_dist(c::AbstractCosmology, z; kws...) = hubble_dist0(c)*Z(c, z; kws...)
 
-comoving_transverse_dist_mpc(c::AbstractFlatCosmology, z; kws...) =
-    comoving_radial_dist_mpc(c, z; kws...)
-function comoving_transverse_dist_mpc(c::AbstractOpenCosmology, z; kws...)
+comoving_transverse_dist(c::AbstractFlatCosmology, z; kws...) =
+    comoving_radial_dist(c, z; kws...)
+function comoving_transverse_dist(c::AbstractOpenCosmology, z; kws...)
     sqrtok = sqrt(c.Ω_k)
-    hubble_dist_mpc0(c)*sinh(sqrtok*Z(c,z; kws...))/sqrtok
+    hubble_dist0(c)*sinh(sqrtok*Z(c,z; kws...))/sqrtok
 end
-function comoving_transverse_dist_mpc(c::AbstractClosedCosmology, z; kws...)
+function comoving_transverse_dist(c::AbstractClosedCosmology, z; kws...)
     sqrtok = sqrt(abs(c.Ω_k))
-    hubble_dist_mpc0(c)*sin(sqrtok*Z(c,z; kws...))/sqrtok
+    hubble_dist0(c)*sin(sqrtok*Z(c,z; kws...))/sqrtok
 end
 
-angular_diameter_dist_mpc(c::AbstractCosmology, z; kws...) =
-    comoving_transverse_dist_mpc(c, z; kws...)/(1 + z)
+angular_diameter_dist(c::AbstractCosmology, z; kws...) =
+    comoving_transverse_dist(c, z; kws...)/(1 + z)
 
-luminosity_dist_mpc(c::AbstractCosmology, z; kws...) =
-    comoving_transverse_dist_mpc(c, z; kws...)*(1 + z)
+luminosity_dist(c::AbstractCosmology, z; kws...) =
+    comoving_transverse_dist(c, z; kws...)*(1 + z)
 
 distmod(c::AbstractCosmology, z; kws...) =
-    5.0 * log10(luminosity_dist_mpc(c, z; kws...)) + 25.0
+    5.0 * log10(luminosity_dist(c, z; kws...) / Mpc) + 25.0
 
 # volumes
 
-comoving_volume_gpc3(c::AbstractFlatCosmology, z; kws...) =
-    (4pi/3)*(comoving_radial_dist_mpc(c,z; kws...)*1e-3)^3
-function comoving_volume_gpc3(c::AbstractOpenCosmology, z; kws...)
-    DH = hubble_dist_mpc0(c)
-    x = comoving_transverse_dist_mpc(c,z; kws...)/DH
+comoving_volume(c::AbstractFlatCosmology, z; kws...) =
+    (4pi/3)*(comoving_radial_dist(c,z; kws...))^3
+function comoving_volume(c::AbstractOpenCosmology, z; kws...)
+    DH = hubble_dist0(c)
+    x = comoving_transverse_dist(c,z; kws...)/DH
     sqrtok = sqrt(c.Ω_k)
-    2pi*(DH*1e-3)^3*(x*sqrt(1. + c.Ω_k*x^2) - asinh(sqrtok*x)/sqrtok)/c.Ω_k
+    2pi*(DH)^3*(x*sqrt(1. + c.Ω_k*x^2) - asinh(sqrtok*x)/sqrtok)/c.Ω_k
 end
-function comoving_volume_gpc3(c::AbstractClosedCosmology, z; kws...)
-    DH = hubble_dist_mpc0(c)
-    x = comoving_transverse_dist_mpc(c,z; kws...)/DH
+function comoving_volume(c::AbstractClosedCosmology, z; kws...)
+    DH = hubble_dist0(c)
+    x = comoving_transverse_dist(c,z; kws...)/DH
     sqrtok = sqrt(abs(c.Ω_k))
-    2pi*(DH*1e-3)^3*(x*sqrt(1. + c.Ω_k*x^2) - asin(sqrtok*x)/sqrtok)/c.Ω_k
+    2pi*(DH)^3*(x*sqrt(1. + c.Ω_k*x^2) - asin(sqrtok*x)/sqrtok)/c.Ω_k
 end
 
-comoving_volume_element_gpc3(c::AbstractCosmology, z; kws...) =
-    1e-9*hubble_dist_mpc0(c,z)*angular_diameter_dist_mpc(c,z; kws...)^2/a2E(c,scale_factor(z))
+comoving_volume_element(c::AbstractCosmology, z; kws...) =
+    hubble_dist0(c,z)*angular_diameter_dist(c,z; kws...)^2/a2E(c,scale_factor(z))
 
 # times
 
 T(c::AbstractCosmology, a0, a1; kws...) = QuadGK.quadgk(x->x/a2E(c,x), a0, a1; kws...)[1]
-age_gyr(c::AbstractCosmology, z; kws...) = hubble_time_gyr0(c)*T(c, 0., scale_factor(z); kws...)
-lookback_time_gyr(c::AbstractCosmology, z; kws...) = hubble_time_gyr0(c)*T(c, scale_factor(z), 1.; kws...)
+age(c::AbstractCosmology, z; kws...) = hubble_time0(c)*T(c, 0.0, scale_factor(z); kws...)
+lookback_time(c::AbstractCosmology, z; kws...) = hubble_time0(c)*T(c, scale_factor(z), 1.0; kws...)
 
 end # module

test/runtests.jl

@@ -1,5 +1,5 @@
 using Cosmology
-using Test
+using Test, Unitful, UnitfulAstro
 
 # values from http://icosmos.co.uk/
 
@@ -8,81 +8,81 @@ age_rtol = 2e-4
 
 @testset "FlatLCDM" begin
     c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1651.9145 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3303.829 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 151.0571 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6607.6579 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1651.9145u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1,rtol=dist_rtol) ≈ 3303.829u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 151.0571u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6607.6579u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.1002 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.4694 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.7527 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.4694-5.7527 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 13.4694u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.7527u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (13.4694-5.7527)u"Gyr" rtol = age_rtol
 end
 
 @testset "OpenLCDM" begin
     c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1619.9588 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3209.784 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 140.0856 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6479.8352 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1619.9588u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1,rtol=dist_rtol) ≈ 3209.784u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 140.0856u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6479.8352u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.0578 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.064 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.5466 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.064-5.5466 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 13.064u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.5466u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (13.064-5.5466)u"Gyr" rtol = age_rtol
 end
 
 @testset "ClosedLCDM" begin
     c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1686.5272 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3408.937 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 163.8479 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6746.1088 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1686.5272u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1,rtol=dist_rtol) ≈ 3408.937u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 163.8479u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6746.1088u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.1453 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.925 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.9868 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.925-5.9868 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 13.925u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.9868u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (13.925-5.9868)u"Gyr" rtol = age_rtol
 end
 
 @testset "FlatWCDM" begin
     c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1612.0585 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3224.1169 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 140.3851 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6448.2338 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1612.0585u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1,rtol=dist_rtol) ≈ 3224.1169u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 140.3851u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6448.2338u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.0472 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.1915 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.6464 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.1915-5.6464 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 13.1915u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.6464u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (13.1915-5.6464)u"Gyr" rtol = age_rtol
 end
 
 @testset "OpenWCDM" begin
     c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1588.0181 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3147.6227 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 132.0466 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6352.0723 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1588.0181u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,rtol=dist_rtol,1) ≈ 3147.6227u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 132.0466u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6352.0723u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.0146 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 12.8488 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.4659 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 12.8488-5.4659 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 12.8488u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.4659u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (12.8488-5.4659)u"Gyr" rtol = age_rtol
 end
 
 @testset "ClosedWCDM" begin
     c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1637.5993 rtol = dist_rtol
-    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3307.9932 rtol = dist_rtol
-    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 149.8301 rtol = dist_rtol
-    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6550.3973 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1637.5993u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1,rtol=dist_rtol) ≈ 3307.9932u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1,rtol=dist_rtol) ≈ 149.8301u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1,rtol=dist_rtol) ≈ 6550.3973u"Mpc" rtol = dist_rtol
     @test distmod(c,1,rtol=dist_rtol) ≈ 44.0813 rtol = dist_rtol
-    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.5702 rtol = age_rtol
-    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.8482 rtol = age_rtol
-    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.5702-5.8482 rtol = age_rtol
+    @test age(c,0,rtol=age_rtol) ≈ 13.5702u"Gyr" rtol = age_rtol
+    @test age(c,1,rtol=age_rtol) ≈ 5.8482u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1,rtol=age_rtol) ≈ (13.5702-5.8482)u"Gyr" rtol = age_rtol
 end
 
 @testset "Non-Float64" begin
     # Test that FlatLCDM works with non-Float64 (BigFloat in this example)
     c = cosmology(h=0.7, OmegaM=big(0.3), OmegaR=0)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1651.9145 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1651.9145u"Mpc" rtol = dist_rtol
     # Test that FlatWCDM works with non-Float64 (BigFloat in this example)
     c = cosmology(h=big(0.7), OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1612.0585 rtol = dist_rtol
+    @test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1612.0585u"Mpc" rtol = dist_rtol
 end