Fix bug in Z_from_MH

Sign error in derivation led to correct answers at low MH (<-1) but divergences from the true answer at higher metallicities. Errors were minor until super-solar Z which is where I noticed them. This fixes the bug and updates the relevant tests.

src/utilities.jl

@@ -148,18 +148,18 @@ function Z_from_MH(MH, solZ=0.01524; Y_p = 0.2485, γ = 1.78)
     # Z/X = exp10( [M/H] + log(Z/X)☉ )
     # X = 1 - Y - Z
     # Y ≈ Y_p + γ * Z for parsec (see Y_from_Z above)
-    # so X ≈ 1 - (Y_p + γ * Z) - Z = 1 - Y_p + (1 + γ) * Z
+    # so X ≈ 1 - (Y_p + γ * Z) - Z = 1 - Y_p - (1 + γ) * Z
     # Substitute into line 2,
-    # Z / (1 - Y_p + (1 + γ) * Z) = exp10( [M/H] + log(Z/X)☉ )
-    # Z = (1 - Y_p + (1 + γ) * Z) * exp10( [M/H] + log(Z/X)☉ )
+    # Z / (1 - Y_p - (1 + γ) * Z) = exp10( [M/H] + log(Z/X)☉ )
+    # Z = (1 - Y_p - (1 + γ) * Z) * exp10( [M/H] + log(Z/X)☉ )
     # let A = exp10( [M/H] + log(Z/X)☉ )
-    # Z = (1 - Y_p) * A + (1 + γ) * Z * A
-    # Z - (1 + γ) * Z * A = (1 - Y_p) * A
-    # Z (1 - (1 + γ) * A) = (1 - Y_p) * A
-    # Z = (1 - Y_p) * A * (1 - (1 + γ) * A)
+    # Z = (1 - Y_p) * A - (1 + γ) * Z * A
+    # Z + (1 + γ) * Z * A = (1 - Y_p) * A
+    # Z (1 + (1 + γ) * A) = (1 - Y_p) * A
+    # Z = (1 - Y_p) * A / (1 + (1 + γ) * A)
     # Originally had X_from_Z(solZ) without passing through the Y_p. Don't remember why
-    zoverx = exp10(MH + log10(solZ / X_from_Z(solZ, Y_p, γ))) 
-    return (1 - Y_p) * zoverx * (1 - (1 + γ) * zoverx)
+    zoverx = exp10(MH + log10(solZ / X_from_Z(solZ, Y_p, γ)))
+    return (1 - Y_p) * zoverx / (1 + (1 + γ) * zoverx)
 end
 
 # PARSEC says that the solar Z is 0.0152 and Z/X = 0.0207, but they don't quite agree

test/fitting/log_amr_test.jl

@@ -17,20 +17,20 @@ using Test
     high_constraint = (-1.0, 0.0)
     T_max = last(low_constraint)
     α, β = SFH.calculate_αβ_logamr( low_constraint, high_constraint, T_max, SFH.Z_from_MH )
-    @test α ≈ 0.00011345544581771879
-    @test β ≈ 5.106276398722378e-5
+    @test α ≈ 0.00011345962863988529
+    @test β ≈ 5.106276580989658e-5
     @test SFH.Z_from_MH( min(first(low_constraint),first(high_constraint)) ) ≈ β
     @test α ≈ ( SFH.Z_from_MH( first(high_constraint) ) - # dZ / dt
         SFH.Z_from_MH( first(low_constraint) ) ) /
         (last(low_constraint) - last(high_constraint))
     # Test that passing different max_age works
     @test all(SFH.calculate_αβ_logamr( low_constraint,
-                                       high_constraint, 14.0) .≈ (0.00011345544581771879, 1.702613024190806e-5))
+                                       high_constraint, 14.0) .≈ (0.00011345962863988529, 1.7024877217930916e-5))
     # Test that passing different function to calculate Z from MH works
     @test all(SFH.calculate_αβ_logamr( low_constraint,
                                        high_constraint, T_max,
                                        x -> SFH.Z_from_MH(x, 0.017; Y_p = 0.25) ) .≈
-                                           (0.00012735499210578944, 5.736184884707544e-5) )
+                                           (0.0001273609414391884, 5.736185144138171e-5) )
 end
 
 @testset "calculate_coeffs_logamr" begin

test/runtests.jl

@@ -722,9 +722,12 @@ const rtols = (1e-3, 1e-7) # Relative tolerance levels to use for the above floa
                 @test SFH.X_from_Z(convert(T,1e-3), convert(T,0.25)) ≈ 0.74722 rtol=rtols[i] # Return type not guaranteed
                 @test SFH.X_from_Z(convert(T,1e-3), convert(T,0.25), convert(T,1.78)) ≈ 0.74722 rtol=rtols[i] # Return type not guaranteed
                 @test SFH.MH_from_Z(convert(T,1e-3), convert(T,0.01524)) ≈ -1.206576807011171 rtol=rtols[i] # Return type not guaranteed
-                @test SFH.Z_from_MH(convert(T,-2), convert(T,0.01524); Y_p=convert(T,0.2485)) ≈ 0.00016140865968917453 rtol=rtols[i] # Return type not guaranteed
+                @test SFH.Z_from_MH(convert(T,-2), convert(T,0.01524); Y_p=convert(T,0.2485)) ≈ 0.00016140871730361718 rtol=rtols[i] # Return type not guaranteed
                 # These two functions should be inverses; test that they are
                 @test SFH.MH_from_Z(SFH.Z_from_MH(convert(T,-2), convert(T,0.01524); Y_p=convert(T,0.2485), γ=convert(T,1.78)), convert(T,0.01524); Y_p=convert(T,0.2485)) ≈ -2 rtol=rtols[i] # Return type not guaranteed
+                # Due to a previous bug, Z_from_MH passed the above test but diverged at higher Z. Test at positive [M/H] here.
+                @test SFH.MH_from_Z(SFH.Z_from_MH(convert(T,1.0), convert(T,0.01524); Y_p=convert(T,0.2485), γ=convert(T,1.78)), convert(T,0.01524); Y_p=convert(T,0.2485)) ≈ 1.0 rtol=rtols[i] # Return type not guaranteed
+
                 @test SFH.dMH_dZ(convert(T,1e-3), convert(T,0.01524); Y_p = convert(T,0.2485), γ = convert(T,1.78)) ≈ 435.9070188458886 rtol=rtols[i] # Return type not guaranteed
                 @test SFH.Martin2016_complete(T[20.0, 1.0, 25.0, 1.0]...) ≈ big"0.9933071490757151444406380196186748196062559910927034697307877569401159160854199" rtol=rtols[i]
                 @test SFH.Martin2016_complete(T[20.0, 1.0, 25.0, 1.0]...) isa T