Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Computing distances from z1 to z2 #24

Merged
merged 2 commits into from
Dec 16, 2019
Merged

Computing distances from z1 to z2 #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6dc5a6b
Computing distances from z1 to z2
astrozot Dec 10, 2019
e4b7038
Readme and tests
astrozot Dec 16, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,11 @@ Distances
<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>angular_diameter_dist(cosmo,&nbsp;z<sub>1</sub>,&nbsp;z<sub>2</sub>)</td>
<td>Ratio of the proper transverse size (in Mpc) of an object at redshift z<sub>2</sub> to its
angular size (in radians), as seen by an observer at z<sub>1</sub></td>
</tr>
<tr>
<td>comoving_radial_dist(cosmo,&nbsp;z)</td>
<td>Comoving radial distance to redshift z, in Mpc</td>
Expand All @@ -107,6 +112,9 @@ FlatLCDM(0.69,0.7399122024007928,0.26,8.779759920715362e-5)

julia> angular_diameter_dist(c, 1.2)
1784.0089227105113 Mpc

julia> angular_diameter_dist(c, 0.7, 1.2)
606.6521737365097 Mpc
```

For each function returning a unitful number, you can specify a different unit
Expand Down
17 changes: 16 additions & 1 deletion src/Cosmology.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,22 +148,37 @@ hubble_time(c::AbstractCosmology, z) = hubble_time0(c)/E(c,z)

Z(c::AbstractCosmology, z::Real; kws...) =
QuadGK.quadgk(a->1/a2E(c,a), scale_factor(z), 1; kws...)[1]
Z(c::AbstractCosmology, z₁::Real, z₂::Real; kws...) =
giordano marked this conversation as resolved.
Show resolved Hide resolved
QuadGK.quadgk(a->1/a2E(c,a), scale_factor(z₂), scale_factor(z₁); kws...)[1]

comoving_radial_dist(c::AbstractCosmology, z; kws...) = hubble_dist0(c)*Z(c, z; kws...)
comoving_radial_dist(c::AbstractCosmology, z₁, z₂; kws...) = hubble_dist0(c)*Z(c, z₁, z₂; kws...)
giordano marked this conversation as resolved.
Show resolved Hide resolved

comoving_transverse_dist(c::AbstractFlatCosmology, z; kws...) =
comoving_radial_dist(c, z; kws...)
comoving_transverse_dist(c::AbstractFlatCosmology, z₁, z₂; kws...) =
giordano marked this conversation as resolved.
Show resolved Hide resolved
comoving_radial_dist(c, z₁, z₂; kws...)
function comoving_transverse_dist(c::AbstractOpenCosmology, z; kws...)
sqrtok = sqrt(c.Ω_k)
hubble_dist0(c)*sinh(sqrtok*Z(c,z; kws...))/sqrtok
hubble_dist0(c)*sinh(sqrtok*Z(c, z; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractOpenCosmology, z₁, z₂; kws...)
giordano marked this conversation as resolved.
Show resolved Hide resolved
sqrtok = sqrt(c.Ω_k)
hubble_dist0(c)*sinh(sqrtok*Z(c, z₁, z₂; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractClosedCosmology, z; kws...)
sqrtok = sqrt(abs(c.Ω_k))
hubble_dist0(c)*sin(sqrtok*Z(c,z; kws...))/sqrtok
end
function comoving_transverse_dist(c::AbstractClosedCosmology, z₁, z₂; kws...)
giordano marked this conversation as resolved.
Show resolved Hide resolved
sqrtok = sqrt(abs(c.Ω_k))
hubble_dist0(c)*sin(sqrtok*Z(c, z₁, z₂; kws...))/sqrtok
end

angular_diameter_dist(c::AbstractCosmology, z; kws...) =
comoving_transverse_dist(c, z; kws...)/(1 + z)
angular_diameter_dist(c::AbstractCosmology, z₁, z₂; kws...) =
giordano marked this conversation as resolved.
Show resolved Hide resolved
comoving_transverse_dist(c, z₁, z₂; kws...)/(1 + z₂)

luminosity_dist(c::AbstractCosmology, z; kws...) =
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should luminosity_dist and distmod have the new method?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably not: I never saw a luminosity distance between two different redshifts, and I cannot easily imagine a situation where such a quantity would be used. But, for completeness, one could also include it. FYI, the Python astrolib package does not include it.

Anyway, let me know your decision, and I will modify the code (and the Readme file) accordingly.

P.S. Not sure how it works: if I make a new commit to the code and push it to the repository, is the pull request automatically updated or should I make a new pull request?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, no need to add these methods, but please add tests. You can just keep pushing to the same branch, the pull request will be updated

comoving_transverse_dist(c, z; kws...)*(1 + z)
Expand Down
8 changes: 8 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ integrand(c, z) = 4pi*ustrip(comoving_volume_element(c, z))
@testset "FlatLCDM" begin
c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
@test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1651.9145u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 625.3444u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -26,6 +27,7 @@ end
@testset "OpenLCDM" begin
c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
@test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1619.9588u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 598.9118u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -40,6 +42,7 @@ end
@testset "ClosedLCDM" begin
c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
@test angular_diameter_dist(c,1,rtol=dist_rtol) ≈ 1686.5272u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 655.6019u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -54,6 +57,7 @@ end
@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,rtol=dist_rtol) ≈ 1612.0585u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 607.6802u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -68,6 +72,7 @@ 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(c,1,rtol=dist_rtol) ≈ 1588.0181u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 585.4929u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -82,6 +87,7 @@ 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(c,1,rtol=dist_rtol) ≈ 1637.5993u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 632.5829u"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 quadgk(z -> integrand(c, z), 0, 2.5)[1] ≈ ustrip(comoving_volume(c, 2.5))
Expand All @@ -97,10 +103,12 @@ end
# 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,rtol=dist_rtol) ≈ 1651.9145u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 625.3444u"Mpc" rtol = dist_rtol
@test comoving_volume_element(c, big(1.41)) ≈ 3.4030879e10u"Mpc^3" 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,rtol=dist_rtol) ≈ 1612.0585u"Mpc" rtol = dist_rtol
@test angular_diameter_dist(c,1,2,rtol=dist_rtol) ≈ 607.6802u"Mpc" rtol = dist_rtol
@test comoving_volume_element(c, big(1.41)) ≈ 3.1378625e10u"Mpc^3" rtol = dist_rtol
end

Expand Down