Merge 1a477d4 into fd106f8

REQUIRE

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

src/Cosmology.jl

@@ -2,20 +2,20 @@ __precompile__()
 
 module Cosmology
 
-using QuadGK
+using QuadGK, Unitful, UnitfulAstro
 
 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,64 +133,64 @@ 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) * u"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 * u"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 * u"Gyr"
+hubble_time(c::AbstractCosmology, z) = hubble_time0(c)/E(c,z)
 
 # distances
 
 Z(c::AbstractCosmology, z::Real) = ((q,_) = QuadGK.quadgk(a::Float64->1.0/a2E(c,a), scale_factor(z), 1); q)
 
-comoving_radial_dist_mpc(c::AbstractCosmology, z) = hubble_dist_mpc0(c)*Z(c, z)
+comoving_radial_dist(c::AbstractCosmology, z) = hubble_dist0(c)*Z(c, z)
 
-comoving_transverse_dist_mpc(c::AbstractFlatCosmology, z) =
-    comoving_radial_dist_mpc(c, z)
-function comoving_transverse_dist_mpc(c::AbstractOpenCosmology, z)
+comoving_transverse_dist(c::AbstractFlatCosmology, z) =
+    comoving_radial_dist(c, z)
+function comoving_transverse_dist(c::AbstractOpenCosmology, z)
     sqrtok = sqrt(c.Ω_k)
-    hubble_dist_mpc0(c)*sinh(sqrtok*Z(c,z))/sqrtok
+    hubble_dist0(c)*sinh(sqrtok*Z(c,z))/sqrtok
 end
-function comoving_transverse_dist_mpc(c::AbstractClosedCosmology, z)
+function comoving_transverse_dist(c::AbstractClosedCosmology, z)
     sqrtok = sqrt(abs(c.Ω_k))
-    hubble_dist_mpc0(c)*sin(sqrtok*Z(c,z))/sqrtok
+    hubble_dist0(c)*sin(sqrtok*Z(c,z))/sqrtok
 end
 
-angular_diameter_dist_mpc(c::AbstractCosmology, z) =
-    comoving_transverse_dist_mpc(c, z)/(1 + z)
+angular_diameter_dist(c::AbstractCosmology, z) =
+    comoving_transverse_dist(c, z)/(1 + z)
 
-luminosity_dist_mpc(c::AbstractCosmology, z) =
-    comoving_transverse_dist_mpc(c, z)*(1 + z)
+luminosity_dist(c::AbstractCosmology, z) =
+    comoving_transverse_dist(c, z)*(1 + z)
 
 distmod(c::AbstractCosmology, z) =
-    5.0 * log10(luminosity_dist_mpc(c, z)) + 25.0
+    5.0 * log10(luminosity_dist(c, z) / u"Mpc") + 25.0
 
 # volumes
 
-comoving_volume_gpc3(c::AbstractFlatCosmology, z) =
-    (4pi/3)*(comoving_radial_dist_mpc(c,z)*1e-3)^3
-function comoving_volume_gpc3(c::AbstractOpenCosmology, z)
-    DH = hubble_dist_mpc0(c)
-    x = comoving_transverse_dist_mpc(c,z)/DH
+comoving_volume(c::AbstractFlatCosmology, z) =
+    (4pi/3)*(comoving_radial_dist(c,z))^3
+function comoving_volume(c::AbstractOpenCosmology, z)
+    DH = hubble_dist0(c)
+    x = comoving_transverse_dist(c,z)/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)
-    DH = hubble_dist_mpc0(c)
-    x = comoving_transverse_dist_mpc(c,z)/DH
+function comoving_volume(c::AbstractClosedCosmology, z)
+    DH = hubble_dist0(c)
+    x = comoving_transverse_dist(c,z)/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) =
-    1e-9*hubble_dist_mpc0(c,z)*angular_diameter_dist_mpc(c,z)^2/a2E(c,scale_factor(z))
+comoving_volume_element(c::AbstractCosmology, z) =
+    hubble_dist0(c,z)*angular_diameter_dist(c,z)^2/a2E(c,scale_factor(z))
 
 # times
 
 T(c::AbstractCosmology, a0, a1) = ((q,_) = QuadGK.quadgk(x->x/a2E(c,x), a0, a1); q)
-age_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)*T(c, 0., scale_factor(z))
-lookback_time_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)*T(c, scale_factor(z), 1.)
+age(c::AbstractCosmology, z) = hubble_time0(c)*T(c, 0., scale_factor(z))
+lookback_time(c::AbstractCosmology, z) = hubble_time0(c)*T(c, scale_factor(z), 1.)
 
 end # module

test/runtests.jl

@@ -1,91 +1,88 @@
 using Cosmology
-using Test
-
-function test_approx_eq_rtol(va, vb, rtol, astr, bstr)
-    diff = maximum(abs(va - vb))
-    tol = rtol*max(maximum(abs(va)), maximum(abs(vb)))
-    if diff > tol
-        sdiff = string("|", astr, " - ", bstr, "| <= ", tol)
-        error("assertion failed: ", sdiff,
-              "\n  ", astr, " = ", va,
-              "\n  ", bstr, " = ", vb,
-              "\n  difference = ", diff, " > ", tol)
-    end
-end
-
-macro test_approx_eq_rtol(a, b, c)
-    :(test_approx_eq_rtol($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
-end
+using Test, Unitful, UnitfulAstro
 
 # values from http://icosmos.co.uk/
 
 dist_rtol = 1e-6
 age_rtol = 2e-4
 
-c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1651.9145 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3303.829 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 151.0571 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6607.6579 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.1002 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 13.4694 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.7527 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 13.4694-5.7527 age_rtol
-
-c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1619.9588 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3209.784 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 140.0856 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6479.8352 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.0578 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 13.064 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.5466 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 13.064-5.5466 age_rtol
+@testset "FlatLCDM" begin
+    c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
+    @test angular_diameter_dist(c,1) ≈ 1651.9145u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3303.829u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 151.0571u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6607.6579u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.1002 rtol = dist_rtol
+    @test age(c,0) ≈ 13.4694u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.7527u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (13.4694-5.7527)u"Gyr" rtol = age_rtol
+end
 
-c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1686.5272 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3408.937 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 163.8479 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6746.1088 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.1453 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 13.925 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.9868 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 13.925-5.9868 age_rtol
+@testset "OpenLCDM" begin
+    c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
+    @test angular_diameter_dist(c,1) ≈ 1619.9588u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3209.784u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 140.0856u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6479.8352u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.0578 rtol = dist_rtol
+    @test age(c,0) ≈ 13.064u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.5466u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (13.064-5.5466)u"Gyr" rtol = age_rtol
+end
 
-c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1612.0585 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3224.1169 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 140.3851 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6448.2338 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.0472 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 13.1915 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.6464 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 13.1915-5.6464 age_rtol
+@testset "ClosedLCDM" begin
+    c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
+    @test angular_diameter_dist(c,1) ≈ 1686.5272u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3408.937u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 163.8479u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6746.1088u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.1453 rtol = dist_rtol
+    @test age(c,0) ≈ 13.925u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.9868u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (13.925-5.9868)u"Gyr" rtol = age_rtol
+end
 
-c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1588.0181 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3147.6227 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 132.0466 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6352.0723 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.0146 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 12.8488 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.4659 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 12.8488-5.4659 age_rtol
+@testset "FlatWCDM" begin
+    c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
+    @test angular_diameter_dist(c,1) ≈ 1612.0585u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3224.1169u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 140.3851u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6448.2338u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.0472 rtol = dist_rtol
+    @test age(c,0) ≈ 13.1915u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.6464u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (13.1915-5.6464)u"Gyr" rtol = age_rtol
+end
 
-c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1637.5993 dist_rtol
-@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3307.9932 dist_rtol
-@test_approx_eq_rtol comoving_volume_gpc3(c,1) 149.8301 dist_rtol
-@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6550.3973 dist_rtol
-@test_approx_eq_rtol distmod(c,1) 44.0813 dist_rtol
-@test_approx_eq_rtol age_gyr(c,0) 13.5702 age_rtol
-@test_approx_eq_rtol age_gyr(c,1) 5.8482 age_rtol
-@test_approx_eq_rtol lookback_time_gyr(c,1) 13.5702-5.8482 age_rtol
+@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(c,1) ≈ 1588.0181u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3147.6227u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 132.0466u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6352.0723u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.0146 rtol = dist_rtol
+    @test age(c,0) ≈ 12.8488u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.4659u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (12.8488-5.4659)u"Gyr" rtol = age_rtol
+end
 
-# Test that FlatLCDM works with non-Float64 (BigFloat in this example)
-c = cosmology(h=0.7, OmegaM=big(0.3), OmegaR=0)
-@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1651.9145 dist_rtol
+@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(c,1) ≈ 1637.5993u"Mpc" rtol = dist_rtol
+    @test comoving_radial_dist(c,1) ≈ 3307.9932u"Mpc" rtol = dist_rtol
+    @test comoving_volume(c,1) ≈ 149.8301u"Gpc^3" rtol = dist_rtol
+    @test luminosity_dist(c,1) ≈ 6550.3973u"Mpc" rtol = dist_rtol
+    @test distmod(c,1) ≈ 44.0813 rtol = dist_rtol
+    @test age(c,0) ≈ 13.5702u"Gyr" rtol = age_rtol
+    @test age(c,1) ≈ 5.8482u"Gyr" rtol = age_rtol
+    @test lookback_time(c,1) ≈ (13.5702-5.8482)u"Gyr" rtol = age_rtol
+end
 
-# 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_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1612.0585 dist_rtol
+@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(c,1) ≈ 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(c,1) ≈ 1612.0585u"Mpc" rtol = dist_rtol
+end