ruthangus / stardate Goto Github PK

import stardate as sd

# Create a dictionary of observables
iso_params = {"teff": (4386, 50),     # Teff with uncertainty.
              "logg": (4.66, .05),    # logg with uncertainty.
              "feh": (0.0, .02),      # Metallicity with uncertainty.
              "parallax": (1.48, .1),  # Parallax in milliarcseconds.
              "maxAV": .1}            # Maximum extinction

prot, prot_err = 29, 3

# Set up the star object.
star = sd.Star(iso_params, prot=prot, prot_err=prot_err)  # Here's where you add a rotation period

# Run the MCMC
star.fit(max_n=1000)

# max_n is the maximum number of MCMC samples. I recommend setting this
# much higher when running for real, or using the default value of 100000.

# Print the median age with the 16th and 84th percentile uncertainties.
age, errp, errm, samples = star.age_results()
print("stellar age = {0:.2f} + {1:.2f} + {2:.2f}".format(age, errp, errm))

platform darwin -- Python 3.8.13, pytest-7.1.2, pluggy-1.0.0
rootdir: /Users/caroline/Research/GitHub/stardate
collected 13 items                                                                                    

test/isochrone_test.py F                                                                        [  7%]
test/lhf_test.py ..F..F......                                                                   [100%]

============================================== FAILURES ===============================================
___________________________________________ test_iso_lnlike ___________________________________________

    def test_iso_lnlike():
    
        iso_params = {"teff": (5770, 10),
                    "feh": (0, .01),
                    "logg": (4.44, .1),
                    "parallax": (1, 1)}
    
        mod = SingleStarModel(mist, **iso_params)  # StarModel isochrones obj
        params = [354, np.log10(4.56*1e9), 0., 1000, 0.]
        lnparams = [354, np.log10(4.56*1e9), 0., np.log(1000), 0.]
        lnpr = mod.lnprior(params)
        # lnparams = [350, 9, 0., 6, 0.]
        args1 = [mod, iso_params]
    
        # Calculate the lnprob above
        iso_lp = iso_lnprob(lnparams, *args1)
        start = time.time()
        for i in range(100):
            iso_lp = iso_lnprob(lnparams, *args1)
        end = time.time()
    
        # Calculate the stardate lnprob
        args2 = [mod, None, None, True, False, True, "praesepe"]
    
        start = time.time()
        for i in range(100):
            lp = lnprob(lnparams, *args2)[0]
        end = time.time()
        # print("time = ", end - start)
    
        ll = lnlike(lnparams, *args2)
        lnprior = lp - ll
    
        assert iso_lp == lp
        assert np.isclose(iso_lp - lnpr, ll)
>       assert lnpr == lnprior
E       assert -18.46720822861122 == -18.46720822861107

test/isochrone_test.py:62: AssertionError
____________________________________________ test_calc_bv _____________________________________________

    def test_calc_bv():
        sun = [355, np.log10(4.56*1e9), 0., np.log(1000), 0.]
        bv_sun = calc_bv(sun)
>       assert .6 < bv_sun, "Are you sure you're on the isochrones eep branch?"
E       AssertionError: Are you sure you're on the isochrones eep branch?
E       assert 0.6 < 0.42834431703148956

test/lhf_test.py:92: AssertionError
_______________________________________ test_gyro_model_rossby ________________________________________

    def test_gyro_model_rossby():
        """
        Make sure that the rossby model gives Solar values for the Sun.
        Make sure it gives a maximum rotation period of pmax.
        """
        age = np.log10(4.56*1e9)
        sun = [355, age, 0., np.log(1000), 0.]
    
        # iso_params = {"teff": (5777, 10),
        #               "logg": (4.44, .05),
        #               "feh": (0., .001),
        #               "parallax": (1., .01),  # milliarcseconds
        #               "B": (15.48, 0.02)}
    
        # # Set up the StarModel isochrones object.
        # mod = StarModel(mist, **iso_params)
        # # the lnprob arguments]
        # args = [mod, 26., 1., False, True, "praesepe"]
        # print(lnprob(sun, *args))
        # assert 0
    
        prot_sun = 10**gyro_model_rossby(sun, model="praesepe")[0]
>       assert 21 < prot_sun
E       assert 21 < 6.909888089400217

test/lhf_test.py:142: AssertionError