Revision 61a4b1bc63133e04f774b8ab6325f39c05bf7693 authored by Hanno Rein on 20 July 2023, 22:04:30 UTC, committed by Hanno Rein on 20 July 2023, 22:04:30 UTC
1 parent b24da72
Raw File
units.py
# Unit constants
import math
from . import clibrebound
from ctypes import c_char_p, c_uint32

def hash_to_unit(hash):
    clibrebound.reb_hash.restype = c_uint32
    for u in times_SI.keys():
        uhash = clibrebound.reb_hash(c_char_p(u.encode("ascii")))
        if uhash == hash:
            return u
    for u in masses_SI.keys():
        uhash = clibrebound.reb_hash(c_char_p(u.encode("ascii")))
        if uhash == hash:
            return u
    for u in lengths_SI.keys():
        uhash = clibrebound.reb_hash(c_char_p(u.encode("ascii")))
        if uhash == hash:
            return u
    return None

# All units entered in SI (kg, m, s)
G_SI = 6.67408e-11
times_SI = {'s':1.,
    'hr':3600.,
    'day': 86400.,
    'days': 86400.,
    'd': 86400.,
    'yr':31557600., # Julian year (exact)
    'year':31557600., 
    'years':31557600., 
    'yrs':31557600., 
    'jyr':31557600.,
    'sidereal_yr':31558149.7635,
    'yr2pi':math.sqrt(149597870700.**3/1.3271244004193938e20), # chosen to make G=1
    'kyr':31557600.*1.e3,
    'myr':31557600.*1.e6,
    'gyr':31557600.*1.e9}
lengths_SI =  {'m':1.,
    'cm':0.01,
    'km':1000.,
    'au':149597870700.,
    'aus':149597870700.,
    'pc':3.085677581e16,
    'parsec':3.085677581e16
    }

    #What we measure accurately is GM, so set mass units such that G*M gives the value of GM in horizons.py (in the list at the end of horizons.py, the NAIF codes ending in 99 refer to the planets, single digits to the total mass of the planet plus its moons).  Have to multiply by 10**9 since that list has G in kg^-1km^3/s^2 and we use SI.

masses_SI = {'kg':1.,
    'g':1.0e-3,
    'gram':1.0e-3,
    'msun':1.3271244004193938E+11/G_SI*10**9,
    'solarmass':1.3271244004193938E+11/G_SI*10**9,
    'sunmass':1.3271244004193938E+11/G_SI*10**9,
    'msolar':1.3271244004193938E+11/G_SI*10**9,
    'mmercury':2.2031780000000021E+04/G_SI*10**9,
    'mvenus':3.2485859200000006E+05/G_SI*10**9,
    'mearth':3.9860043543609598E+05/G_SI*10**9,
    'mmars':4.282837362069909E+04/G_SI*10**9,
    'mjupiter':1.266865349218008E+08/G_SI*10**9,
    'msaturn':3.793120749865224E+07/G_SI*10**9,
    'muranus':5.793951322279009E+06/G_SI*10**9,
    'mneptune':6.835099502439672E+06/G_SI*10**9,
    'mpluto':8.696138177608748E+02/G_SI*10**9,
    'massist':4.48485856027459e+14/G_SI*10**9, # Sun has mass 0.00029591220828412 in these units. Used to keep G=1 while length=AU and time=day
             }


def units_convert_particle(p, old_l, old_t, old_m, new_l, new_t, new_m):
    p.m = convert_mass(p.m, old_m, new_m)
    p.x = convert_length(p.x, old_l, new_l) 
    p.y = convert_length(p.y, old_l, new_l)
    p.z = convert_length(p.z, old_l, new_l)
    p.vx = convert_vel(p.vx, old_l, old_t, new_l, new_t)
    p.vy = convert_vel(p.vy, old_l, old_t, new_l, new_t)
    p.vz = convert_vel(p.vz, old_l, old_t, new_l, new_t)
    p.ax = convert_acc(p.ax, old_l, old_t, new_l, new_t)
    p.ay = convert_acc(p.ay, old_l, old_t, new_l, new_t)
    p.az = convert_acc(p.az, old_l, old_t, new_l, new_t)
    return p

def convert_mass(mass, old_m, new_m):
    return mass*masses_SI[old_m]/masses_SI[new_m]

def convert_length(length, old_l, new_l):
    return length*lengths_SI[old_l]/lengths_SI[new_l]

def convert_vel(vel, old_l, old_t, new_l, new_t):
    in_SI=vel*lengths_SI[old_l]/times_SI[old_t]
    return in_SI*times_SI[new_t]/lengths_SI[new_l]

def convert_acc(acc, old_l, old_t, new_l, new_t):
    in_SI=acc*lengths_SI[old_l]/times_SI[old_t]**2
    return in_SI*times_SI[new_t]**2/lengths_SI[new_l]

def convert_G(newunits):
    new_l, new_t, new_m = newunits
    return G_SI*masses_SI[new_m]*times_SI[new_t]**2/lengths_SI[new_l]**3
       
def check_units(newunits):   
    if len(newunits) != 3:
        raise Exception("Error: Need to pass exactly 3 units for length, time, and mass (any order), see ipython_examples/Units.ipynb")

    if isinstance(newunits, dict): 
        # keys are not important as they are inferred from the values anyway
        newunits = newunits.values() 
    
    l_unit = t_unit = m_unit = None
    for unit in newunits:
        unit = unit.lower()
        if unit in lengths_SI:
            l_unit = unit
        if unit in times_SI:
            t_unit = unit
        if unit in masses_SI:
            m_unit = unit

    if l_unit is None or t_unit is None or m_unit is None:
        raise Exception("Error: Need to assign rebound.units a tuple consisting of 3 units for length, time, and mass (any order).  See ipython/examples/Units.ipynb.  If you passed such a tuple, at least one of your units isn't in our list.  Please update the dictionaries at the top of rebound/rebound/units.py and send a pull request!")

    return (l_unit, t_unit, m_unit)
back to top