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

Add new helper functions, IrrationalConstants #47

Merged
merged 1 commit into from
Sep 19, 2024
Merged

Add new helper functions, IrrationalConstants #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 3 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ version = "0.1.0"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
KissMCMC = "79d62d8d-4dfd-5781-bc85-ce78e0ac132a"
LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
Expand All @@ -30,8 +31,8 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"

[extensions]
TypedTablesExt = "TypedTables"
DataFramesExt = "DataFrames"
TypedTablesExt = "TypedTables"

[compat]
DataFrames = "1"
Expand All @@ -41,6 +42,7 @@ Documenter = "1"
DynamicHMC = "3.4.2" # Changed stack_posterior_matrices order https://github.com/tpapp/DynamicHMC.jl/pull/175
InitialMassFunctions = "0.1"
Interpolations = "0.13.4, 0.14, 0.15"
IrrationalConstants = "0.2"
KissMCMC = "0.2"
LBFGSB = "0.4"
LineSearches = "7.2"
Expand Down
9 changes: 9 additions & 0 deletions docs/src/helpers.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,15 @@ StarFormationHistories.distance_modulus_to_distance

## Magnitudes and Luminosities

```@docs
StarFormationHistories.mag2flux
StarFormationHistories.flux2mag
StarFormationHistories.magerr
StarFormationHistories.fluxerr
StarFormationHistories.snr_magerr
StarFormationHistories.magerr_snr
```

## [Metallicities](@id metallicity_helpers)

```@docs
Expand Down
1 change: 1 addition & 0 deletions src/StarFormationHistories.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ using Distributions: Distribution, Sampleable, Univariate, Continuous, pdf, logp
import Distributions: _rand! # Extending
import DynamicHMC # For random uncertainties in SFH fits
using Interpolations: interpolate, Gridded, Linear, deduplicate_knots! # extrapolate, Throw
using IrrationalConstants: logten
import LBFGSB # Used for one method in fitting.jl
import LineSearches # For configuration of Optim.jl
# Need mul! for composite!, ∇loglikelihood!;
Expand Down
62 changes: 60 additions & 2 deletions src/utilities.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,66 @@ function angular_transformation_distance(angle, distance0, distance1)
end

#### Luminosity Utilities
"""
mag2flux(m, zpt=0)
Convert a magnitude `m` to a flux assuming a photometric zeropoint of `zpt`, defined as
the magnitude of an object that produces one count (or data number, DN) per second.
```jldoctest; setup = :(import StarFormationHistories: mag2flux)
julia> mag2flux(15.0, 25.0) ≈ exp10(4 * (25.0 - 15.0) / 10)
true
```
"""
mag2flux(m, zpt=0) = exp10(4 * (zpt - m) / 10)
"""
flux2mag(f, zpt=0)
Convert a flux `f` to a magnitude assuming a photometric zeropoint of `zpt`, defined as
the magnitude of an object that produces one count (or data number, DN) per second.
```jldoctest; setup = :(import StarFormationHistories: flux2mag)
julia> flux2mag(10000.0, 25.0) ≈ 25.0 - 5 * log10(10000.0) / 2
true
```
"""
flux2mag(f, zpt=0) = zpt - 5 * log10(f) / 2
"""
magerr(f, σf)
Returns an error in magnitudes given a flux and a flux uncertainty.

```jldoctest; setup = :(import StarFormationHistories: magerr)
julia> magerr(100.0, 1.0) ≈ 2.5 / log(10) * (1.0 / 100.0)
true
```
"""
magerr(f, σf) = 5//2 * σf / f / logten
"""
fluxerr(f, σm)
Returns an error in flux given a flux and a magnitude uncertainty.

```jldoctest; setup = :(import StarFormationHistories: fluxerr)
julia> fluxerr(100.0, 0.01) ≈ (0.01 * 100.0) / 2.5 * log(10)
true
```
"""
fluxerr(f, σm) = σm * f * logten / 5//2
"""
snr_magerr(σm)
Returns a signal-to-noise ratio ``(f/σf)`` given an uncertainty in magnitudes.

```jldoctest; setup = :(import StarFormationHistories: snr_magerr)
julia> snr_magerr(0.01) ≈ 2.5 / log(10) / 0.01
true
```
"""
snr_magerr(σm) = 5//2 / σm / logten
"""
magerr_snr(snr)
Returns a magnitude uncertainty given a signal-to-noise ratio ``(f/σf)``.

```jldoctest; setup = :(import StarFormationHistories: magerr_snr)
julia> magerr_snr(100.0) ≈ 2.5 / log(10) / 100.0
true
```
"""
magerr_snr(snr) = 5//2 / snr / logten
# Absolute magnitude of Sun in V-band is 4.83 = 483//100
L_from_MV(absmagv) = mag2flux(absmagv, 483//100)
MV_from_L(lum) = flux2mag(lum, 483//100)
Expand Down Expand Up @@ -120,7 +178,7 @@ Completeness model of [Martin et al. 2016](https://ui.adsabs.harvard.edu/abs/201

`m` is the magnitude of interest, `A` is the maximum completeness, `m50` is the magnitude at which the data are 50% complete, and `ρ` is an effective slope modifier.
"""
Martin2016_complete(m,A,m50,ρ) = A / (1 + exp((m-m50) / ρ))
Martin2016_complete(m,A,m50,ρ) = A / (1 + exp((m - m50) / ρ))

"""
exp_photerr(m, a, b, c, d)
Expand All @@ -133,7 +191,7 @@ Exponential model for photometric errors of the form

Reported values for some HST data were `a=1.05, b=10.0, c=32.0, d=0.01`.
"""
exp_photerr(m, a, b, c, d) = a^(b * (m-c)) + d
exp_photerr(m, a, b, c, d) = a^(b * (m - c)) + d

"""
process_ASTs(ASTs::Union{DataFrames.DataFrame,
Expand Down
20 changes: 11 additions & 9 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,7 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa

##################################
####### Testing generate_stars_mag
absmaglim = T(-7)
absmaglim = T(-8)
# Test with NoBinaries() model
result = SFH.generate_stars_mag(m_ini, mags, mag_names, absmaglim, mag_names[2], imf; dist_mod=dmod, rng=rng, mag_lim=T(Inf), mag_lim_name=mag_names[2], binary_model=SFH.NoBinaries())
@test result[1] isa Vector{SVector{1,S}}
Expand All @@ -187,7 +187,7 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
# brighter than the requested absmaglim
apparent_mag = SFH.flux2mag(sum(map(x->SFH.mag2flux(x[2]), result[2])))
abs_mag = apparent_mag - dmod
@test T(-0.06) <= (abs_mag - absmaglim) <= T(0)
@test T(-0.05) <= (abs_mag - absmaglim) <= T(0)

# Test with RandomBinaryPairs() model
result = SFH.generate_stars_mag(m_ini, mags, mag_names, absmaglim, mag_names[2], imf; dist_mod=dmod, rng=rng, mag_lim=T(Inf), mag_lim_name=mag_names[2], binary_model=SFH.RandomBinaryPairs(T(4//10)))
Expand All @@ -198,7 +198,7 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
# brighter than the requested absmaglim
apparent_mag = SFH.flux2mag(sum(map(x->SFH.mag2flux(x[2]), result[2])))
abs_mag = apparent_mag - dmod
@test T(-0.06) <= (abs_mag - absmaglim) <= T(0)
@test T(-0.05) <= (abs_mag - absmaglim) <= T(0)

# Test with BinaryMassRatio() model
result = SFH.generate_stars_mag(m_ini, mags, mag_names, absmaglim, mag_names[2], imf; dist_mod=dmod, rng=rng, mag_lim=T(Inf), mag_lim_name=mag_names[2], binary_model=SFH.BinaryMassRatio(T(4//10), Uniform(T(1//10),T(1))))
Expand All @@ -209,7 +209,7 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
# brighter than the requested absmaglim
apparent_mag = SFH.flux2mag(sum(map(x->SFH.mag2flux(x[2]), result[2])))
abs_mag = apparent_mag - dmod
@test T(-0.06) <= (abs_mag - absmaglim) <= T(0)
@test T(-0.05) <= (abs_mag - absmaglim) <= T(0)

# Test errors
@test_throws ArgumentError SFH.generate_stars_mag(m_ini, mags, mag_names, absmaglim, "V", imf; dist_mod=dmod, rng=rng, mag_lim=T(Inf), mag_lim_name=mag_names[2], binary_model=SFH.BinaryMassRatio(T(4//10), Uniform(T(1//10),T(1))))
Expand All @@ -219,7 +219,7 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
# Test generate_stars_mass_composite
# We'll spoof a second isochrone by just shifting
# the F814W mags slightly lower and slightly altering m_ini
composite_masses = [m_ini,m_ini .+ T(0.01)]
composite_masses = [m_ini, m_ini .+ T(0.01)]
composite_mags = [mags, [mags[1], mags[2] .- T(0.02)]]
@test length(composite_masses) == length(composite_mags)
# nisochrones = length(composite_masses)
Expand Down Expand Up @@ -250,11 +250,13 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
@test result[1][i] isa Vector{SVector{1,S}} # Masses
@test result[2][i] isa Vector{SVector{2,T}} # Magnitudes
end
# Test that total magnitude of sampled population is (slightly)
# Test that total magnitude of sampled population is only slightly
# brighter than the requested absmaglim
apparent_mag = SFH.flux2mag( sum( sum(map(x->SFH.mag2flux(x[2]), i)) for i in result[2]) )
abs_mag = apparent_mag - dmod
@test T(-0.06) <= (abs_mag - absmaglim) <= T(0)
# absmag = SFH.flux2mag( sum( sum(map(x->SFH.mag2flux(x[2] - dmod), i)) for i in result[2]) )
# @test T(-0.05) <= (abs_mag - absmaglim) <= T(0)
# Test total flux instead of magnitude
flux_total = sum( sum(map(x->SFH.mag2flux(x[2] - dmod), i)) for i in result[2])
@test 1 ≤ flux_total / SFH.mag2flux(absmaglim) ≤ 1.05

# Test errors
@test_throws ArgumentError SFH.generate_stars_mag_composite(composite_masses, composite_mags, mag_names, absmaglim, "V", T[1//2, 1//2], imf; frac_type="lum", dist_mod=dmod, rng=rng, mag_lim=T(Inf), mag_lim_name=mag_names[2], binary_model=SFH.NoBinaries()) # Test bad absmag_name
Expand Down