Skip to content

Commit

Permalink
Add new helper functions, IrrationalConstants
Browse files Browse the repository at this point in the history
  • Loading branch information
@cgarling
cgarling committed Sep 19, 2024
1 parent 5ca85f7 commit 28d7c33
Show file tree
Hide file tree
Showing 5 changed files with 84 additions and 12 deletions.
4 changes: 3 additions & 1 deletion Project.toml
View file
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
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
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
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
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

0 comments on commit 28d7c33

Please sign in to comment.