natashabatalha / virga

A Cloud Model for Exoplanets and Brown Dwarfs
https://natashabatalha.github.io/virga/
27 stars 13 forks source link

Integrate Peter Gao's Code #37

Open natashabatalha opened 8 months ago

natashabatalha commented 8 months ago

From @sirpetergao

First, I added a few aerosol species to gas_properties and pvap, including CO2, H2SO4, and S8. They don't have metallicity dependencies and the gas_mmr is set to that of water, but the idea is that the user needs to give a below-cloud mixing ratio for these to work properly (though the gas_mmr input is still set as optional for now). For H2SO4 I've hard coded an 80% H2SO4 weight percentage - this should be a function of the water vapor mixing ratio and temperature but I think this is enough for now. Second, for justdoit, I added an extra input argument into the init function of the Atmosphere class, bkg_atm, which is a string that determines what the background gas is; right now there is H2, N2, and CO2; if any other string is input it defaults to H2. Then I split the "true" constants in the constants function (R_GAS, AVOGADRO, K_BOLTZ) into its own function and made 3 other functions, constants_h2, constants_n2, and constants_co2; in these I put eps_k, d_molecule, and c_p_factor, which are different for each of H2, N2, and CO2 background atmospheres. I got the eps_k and d_molecule values for N2 and CO2 from the same source as you, Rosner 2000; for the c_p_factor of CO2 I'm using wikipedia for now, which shows that it's roughly 13/3 (or rather, gamma ~1.3 at room temperature, and c_p_factor = gamma/(gamma-1)). And that's it! I probably missed something but I went through the whole code and these are the parts that stuck out to me. I would love to know if there are other parts that I need to change! Now I just need to put in the refractive indices for those new condensates and make the .mieff files. gaocode.zip

natashabatalha commented 8 months ago

@logan-pearce is going to help integrate this

natashabatalha commented 8 months ago

Tholin:

      pvap_tholin_bars = 10**(7.6106 - 11382./T) !this is actually the KCl vapor pressure, not a real tholin one. 
c     Then convert from bars to dynes/cm^2    
      pvap_tholin = pvap_tholin_bars*1e6 

Soot:

      pvap_soot_bars = 10.d0**(7.6106 - 11382./t) !this is actually the KCl vapor pressure, not a real soot one. 
c     Then convert from bars to dynes/cm^2    
      pvap_soot = pvap_soot_bars*1e6   

Gas Properties from @semoran


def tholin(mw_atmos, mh = 1,gas_mmr=None):
    """Defines properties for Haze with Khare 1984 tholin refractive indices as if it were solid in atmosphere.
    Parameters
    ----------
    mw_atmos : float
        Mean molecular weight of the atmosphere amu
    gas_mmr : float , optional
        Gas mass mixing ratio
        None points to the default value of : 8.80e-7
    mh : float , optional
        Metallicity, Default is 1=1xSolar
    Returns
    -------
    Mean molecular weight of gas,
    Gas mass mixing ratio
    Density of gas cgs
    """
    if mh != 1: raise Exception("Alert: No M/H Dependence in FakeHaze Routine. Consult your local theorist to determine next steps.")
    if isinstance(gas_mmr, type(None)):
        gas_mmr = 1e-7
    gas_mw = 500. #this actually quite low but for testing.
    gas_mmr = gas_mmr * (gas_mw/mw_atmos)
    rho_p =  1.38 #average in g/cm3 from He+2017 --- can't find Khare report, so we use this for now
    return gas_mw, gas_mmr, rho_p