<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Activity Description
see
Step2:  Test equilibrium_temperature
Step4:  Mass-luminosity relation
Step5:  test luminosity_from_mass
Step6:  Equilibrium temperature as a function of mass
Step7:  test equilibrium_temperature_from_mass
Step8:  Results for activities
Step9:  Get equilibrium temperatures of all the planets
Step10:  Calculate T_eq for all exoplanets
<ASSISTANT_TASK:>
Python Code:
def equilibrium_temperature(luminosity, distance, albedo=0):
    
    Calculates the equilibrium temperature of a planet, assuming blackbody radiation
    and thermodynamic equilibrium.
    
    Parameters
    ----------
    luminosity : float
        luminosity of the host star [ergs s**-1]
        
    distance : float
        distance between the star and the planet [cm]
        (assumes circular orbits)
    
    albedo : float
        albedo of the planet (0 for perfect absorber, 1 for perfect reflector)
        
    
    Returns
    -------
    T_eq : float
        Equilibrium temperature of the planet
    
    
    T_eq = (luminosity / distance**2)**.25  \
            * ((1-albedo) / (16 * np.pi * const.sigma_sb.cgs.value))**.25
    return T_eq
Explanation: Activity Description
see: https://docs.google.com/document/d/1ALSBqY4wv7BfCtdi1uFhnFNNlMcIHSg1VooC6rgxvds/edit?usp=sharing
Define and test functions
End of explanation
T_eq_earth = equilibrium_temperature(const.L_sun.cgs.value, units.AU.to(units.cm))
print("Earth equilibrium temperature: ", T_eq_earth, " K")
Explanation: Test equilibrium_temperature:
(there are better ways to do testing using libraries like pytest, nose or unittest)
End of explanation
def luminosity_from_mass(mass):
    
    Convert mass to luminosity using a standard mass-luminosity relation
    
    Parameters
    ----------
    mass : float
        stellar mass [g]
        
    Returns
    -------
    luminosity : float
        stellar luminosity [ergs s**-1]
        
    Notes
    -----
    not vectorized.
        `mass` cannot be an array
    
    
    Sources
    -------
    Maurizio & Cassisi (2005) Evolution of stars and stellar populations
        http://books.google.com/books?id=r1dNzr8viRYC&lpg=PA138&dq=Mass-Luminosity%20relation&lr=&client=firefox-a&pg=PA138#v=onepage&q=&f=false
    
    Nebojsa (2004) Advanced astrophysics
        http://books.google.com/books?id=-ljdYMmI0EIC&lpg=PA19&ots=VdMUIiCdP_&dq=Mass-luminosity%20relation&pg=PA19#v=onepage&q=&f=false
        
    
    # convert to solar units
    mass /= const.M_sun.cgs.value
    
    if mass < .43:
        luminosity = .23 * (mass)**2.3
    elif mass < 2:
        luminosity = (mass)**4
    elif mass < 20:
        luminosity = 1.5 * (mass)**3.5
    else:
        luminosity = 3200
        
    # convert back into cgs units
    luminosity *= const.L_sun.cgs.value
    
    return luminosity
    
Explanation: Mass-luminosity relation
End of explanation
if np.isclose(luminosity_from_mass(const.M_sun.cgs.value) / const.L_sun.cgs.value, 1 ):
    print("OK")
else:
    print("Error: luminosity_from_mass gave wrong answer")
Explanation: test luminosity_from_mass:
End of explanation
def equilibrium_temperature_from_mass(mass, distance, albedo=0):
    return equilibrium_temperature(luminosity_from_mass(mass),
                                   distance,
                                   albedo = albedo)
Explanation: Equilibrium temperature as a function of mass
End of explanation
if np.isclose(equilibrium_temperature_from_mass(const.M_sun.cgs.value, units.AU.to(units.cm)),
           equilibrium_temperature(const.L_sun.cgs.value, units.AU.to(units.cm))):
    print("OK")
else:
    print("Error: equilibrium_temperature_from_mass gave wrong answer")
Explanation: test equilibrium_temperature_from_mass
End of explanation
print("albedo = 0")
T_eq_earth = equilibrium_temperature(const.L_sun.cgs.value, units.AU.to(units.cm))
print("Earth equilibrium temperature: ", T_eq_earth, " K")
print()
Earth_albedo = .3
print("albedo = ", Earth_albedo)
T_eq_earth = equilibrium_temperature(const.L_sun.cgs.value,
                                     units.AU.to(units.cm),
                                     albedo=Earth_albedo)
print("Earth equilibrium temperature: ", T_eq_earth, " K")
mass = .96 * const.M_sun.cgs.value
luminosity = luminosity_from_mass(mass)
T_eq_earth = equilibrium_temperature(luminosity, units.AU.to(units.cm))
print("Earth equilibrium temperature from lower mass sun: ", T_eq_earth)
print("For a 10 K change in T_eq, the sun would need to be ", 
      round(100*(mass / const.M_sun.cgs.value)),
      "% of its current mass")
Explanation: Results for activities
End of explanation
# the better way to do this is with objects,
# but we're not going to teach object-oriented programming in this workshop
# tuples should be ("name", distance [in cm])
planets = [
    ("Mercury", 5.79e12),
    ("Venus",   1.08e13),
    ("Earth", 1.50e13),
    ("Mars", 2.28e13),
    ("Jupiter", 7.78e13),
    ("Saturn", 1.43e14),
    ("Uranus", 2.87e14),
    ("Neptune", 4.50e14),
    ("Pluto", 5.90e14)
]
solar_system_T_eqs = dict()
for planet in planets:
    solar_system_T_eqs[planet[0]] = equilibrium_temperature(const.L_sun.cgs.value,
                                                            planet[1])
solar_system_T_eqs
# if you wanted to keep the solar system in order, 
# you could have used collections.OrderedDict instead of dict
Explanation: Get equilibrium temperatures of all the planets
End of explanation
import exoplanets
exoplanets.download_data()
data = exoplanets.parse_data()
print("column names of `data`: ")
data.dtype.names
luminosities = const.L_sun.cgs.value * 10**data["st_lum"]
distances    = units.AU.to(units.cm) * data["pl_orbsmax"]
plt.scatter(distances / units.AU.to(units.cm), 
            luminosities / const.L_sun.cgs.value)
plt.xscale("log")
plt.yscale("log")
plt.ylim(10**-3, 10**3)
plt.xlabel("Star-Planet distance [AU]")
plt.ylabel("Star luminosity [L_sun]")
equilibrium_temperatures = equilibrium_temperature(luminosities, distances)
plt.hist(equilibrium_temperatures)
plt.xlabel("Equilibrium Temperature [K]")
plt.ylabel("Number of planets")
Explanation: Calculate T_eq for all exoplanets
End of explanation
<END_TASK>