Photo mags

Project.toml

@@ -27,9 +27,16 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TableOperations = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
 Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 
+[weakdeps]
+PhotometricFilters = "482c37c2-cde1-4201-a412-83b81535d120"
+
+[extensions]
+PhotometryExt = "PhotometricFilters"
+
 [compat]
 CSV = "0.8, 0.9, 0.10"
 Compat = "4.10 - 4"
+edFiles = "1.9.1"
 DSP = "0.7, 0.8"
 DataFrames = "1.7.0"
 DiffResults = "1"

ext/PhotometryExt.jl

@@ -0,0 +1,86 @@
+module PhotometryExt
+
+public compute_magnitude
+
+import Interpolations: linear_interpolation
+using Trapz
+import PhotometricFilters
+
+"""
+    compute_magnitude(spec, photometric_filter)
+    compute_magnitude(flux, wavelengths, photometric_filter)
+    compute_magnitude(spec, filter_trans, filter_wave)
+    compute_magnitude(flux, wavelengths, filter_trans, filter_wave)
+
+Compute a photometric magnitude from a synthetic spectrum. This should only be used to create colors or magnitude differences (e.g. relative to Vega).
+
+# Arguments
+- `spec`, an object returned by `Korg.synthesize`. You can pass `flux` and `wavelengths` (in Å) directly instead (which are in Korg's default units).
+- `photometric_filter`, a photometric filter object from `PhotometricFilters.jl`. You can pass `filter_trans` (the throughput) and `filter_wave` (in Å) directly instead.
+
+Adapted from Eddie Schlafly.
+"""
+function compute_magnitude(spec, photometric_filter::PhotometricFilters.PhotometricFilter)
+    compute_magnitude_core(spec.flux, spec.wavelengths, photometric_filter.throughput, photometric_filter.wave)
+end
+function compute_magnitude(spec_flux, spec_wave, photometric_filter::PhotometricFilters.PhotometricFilter)
+    compute_magnitude_core(spec_flux, spec_wave, photometric_filter.throughput, photometric_filter.wave)
+end
+function compute_magnitude(spec, filter_trans, filter_wave)
+    compute_magnitude_core(spec.flux, spec.wavelengths, filter_trans, filter_wave)
+end
+function compute_magnitude(spec_flux, spec_wave, filter_trans, filter_wave)
+    compute_magnitude_core(spec_flux, spec_wave, filter_trans, filter_wave)
+end
+
+"""
+Compute a magnitude in a given filter from a spectrum and filter transmission curve.
+"""
+function compute_magnitude_core(spec_flux, spec_wave, filter_trans, filter_wave)
+    # add atmospheric transmission separately later
+    # maybe I should pass a default resolution on which to compute the integral?
+    @assert length(spec_wave) == length(spec_flux)
+
+    # check spec_wave extends beyond filter_wave
+    if minimum(filter_wave) < minimum(spec_wave) || maximum(filter_wave) > maximum(spec_wave)
+        error("Spectrum does not cover the filter range")
+    end
+    # check spec_wave is longer than filter_wave
+    if length(spec_wave) < length(filter_wave)
+        throw(ArgumentError("Spectrum (length = $(length(spec_wave))) is lower resolution than the filter (length "*
+            "= $(length(filter_wave)). Please use a higher resolution spectrum."))
+    end
+
+    # interpolate filter_wave to spec_wave
+    interp_linear_extrap = linear_interpolation(filter_wave, filter_trans; extrapolation_bc=0)
+    filter_trans_interp = interp_linear_extrap(spec_wave)
+
+    # it is not clear to me that the units work out correctly here even if
+    # one uses erg/s/cm^2/Å (which requires using get_radius)
+    h = 6.62606957e-27  # erg * s
+    c = 2.99792458e10  # cm/s
+    flux_to_number = h * c ./ (spec_wave * 1e-8)  # erg
+
+    spec_ref = 3631e-23 * c ./ (spec_wave * 1e-8) ./ spec_wave # working Mgy
+    Iref = trapz(spec_wave, spec_ref .* filter_trans_interp .* flux_to_number)
+    Ispec = trapz(spec_wave, spec_flux .* filter_trans_interp .* flux_to_number)
+
+    return -2.5 * log10(Ispec / Iref) # AB magnitudes
+end
+
+# returns radius in cm given logg base 10 of cm/s^2
+# assumes 1 solar mass
+# conceptually use this as compute_magnitude(sol.flux.*(get_radius(4.43775056282).^2)./(1e8),sol.wavelengths,t.throughput,t.wave)
+function get_radius(logg)
+    # Constants
+    G = 6.67430e-8  # gravitational constant in cgs
+    M_sun = 1.989e33  # solar mass in g
+    # Surface gravity in cgs from log(g)
+    g = 10^logg # cm/s^2
+    # Using g = GM/R², solve for R
+    # R = sqrt(GM/g)
+    radius = sqrt((G * M_sun) / g)  # in cm
+    return radius
+end
+
+end # module

src/Korg.jl

@@ -26,6 +26,7 @@ include("synth.jl")
 include("prune_linelist.jl")           # select strong lines from a linelist
 include("Fit/Fit.jl")                  # routines to infer stellar params from data
 include("qfactors.jl")                 # formalism to compute theoretical RV precision
+# include("photometry.jl")               # synthetic photometry, moved to extension
 
 @compat public get_APOGEE_DR17_linelist, get_GALAH_DR3_linelist, get_GES_linelist, save_linelist,
                Fit, apply_LSF, compute_LSF_matrix, air_to_vacuum, vacuum_to_air, line_profile,

src/synthesize.jl

@@ -9,8 +9,8 @@ The result of a synthesis. Returned by [`synthesize`](@ref).
 
 # Fields
 
-  - `flux`: the output spectrum
-  - `cntm`: the continuum at each wavelength
+  - `flux`: the output spectrum (in units of erg/s/cm^5)
+  - `cntm`: the continuum at each wavelength (in units of erg/s/cm^5)
   - `intensity`: the intensity at each wavelength and mu value, and possibly each layer in the model
     atmosphere, depending on the radiative transfer scheme.
   - `alpha`: the linear absorption coefficient at each wavelength and atmospheric layer a Matrix of