# Author: Taylor Bell
# Last Update: 2018-11-28
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import astropy.constants as const
from .KeplerOrbit import KeplerOrbit
from .Map import Map
from . import H2_Dissociation_Routines as h2
[docs]class Planet(object):
"""A planet.
Attributes:
a (float): The planet's semi-major axis in m.
absorptivity (float): The absorptivity of the emitting layer (between 0 and 1).
albedo (float): The planet's Bond albedo.
argobliq (float): The reference orbital angle used for the obliquity (in degrees from inferior conjunction).
argp (float): The planet's argument of periastron (in degrees CCW from Omega).
C (float, optional): The planet's heat capacity in J/m^2/K.
cp (float or callable): The planet's isobaric specific heat capacity in J/kg/K.
cpParams (iterable, optional): Any parameters to be passed to cp if using the bell2018 LTE H2+H mix cp
e (float): The planet's orbital eccentricity.
emissivity (float): The emissivity of the emitting layer (between 0 and 1).
g (float): The planet's surface gravity in m/s^2.
inc (float): The planet's orbial inclination (in degrees above face-on)
map (Bell_EBM.Map): The planet's temperature map.
mass (float): The planet's mass in kg.
mlDensity (float): The density of the planet's mixed layer.
mlDepth (float): The depth of the planet's mixed layer (in units of m for rock/water or Pa for gas/bell2018).
obliq (float): The planet's obliquity (axial tilt) (in degrees toward star).
Omega (float): The planet's longitude of ascending node (in degrees CCW from line-of-sight).
orbit (Bell_EBM.KeplerOrbit): The planet's orbit.
plType (str): The planet's composition.
Porb (float): The planet's orbital period in days.
Prot (float): The planet's rotational period in days.
rad (float): The planet's radius in m.
t0 (float): The planet's linear ephemeris in days.
T_exponent (float): The exponent which determinges the rate at which the planet cools (4 for blackbody cooling,
1 for Newtonian cooling) when calculating Fout with bolo=True.
trasmissivity (float): The trasmissivity of the emitting layer (between 0 and 1).
useHealpix (bool): Whether the planet's map uses a healpix grid.
vWind (float): The planet's wind velocity in m/s.
"""
def __init__(self, plType='gas', rad=1, mass=1,
a=0.03, Porb=None, Prot=None, inc=90, t0=0, e=0, Omega=270, argp=90, obliq=0, argobliq=0,
vWind=0, albedo=0, cp=None, mlDepth=None, mlDensity=None, T_exponent=4, emissivity=1, trasmissivity=0,
nlat=16, nlon=None, useHealpix=False, nside=7):
"""Initialization function.
Args:
plType (str, optional): The planet's composition.
rad (float, optional): The planet's radius in m.
mass (float, optional): The planet's mass in kg.
a (float, optional): The planet's semi-major axis in m.
Porb (float, optional): The planet's orbital period in days.
Prot (float, optional): The planet's rotational period in days.
inc (float, optional): The planet's orbial inclination (in degrees above face-on)
t0 (float, optional): The planet's linear ephemeris in days.
e (float, optional): The planet's orbital eccentricity.
Omega (float, optional): The planet's longitude of ascending node (in degrees CCW from line-of-sight).
argp (float, optional): The planet's argument of periastron (in degrees CCW from Omega).
obliq (float, optional): The planet's obliquity (axial tilt) (in degrees toward star).
argobliq (float, optional): The reference orbital angle used for obliq (in degrees from inferior conjunction).
vWind (float, optional): The planet's wind velocity in m/s.
albedo (float, optional): The planet's Bond albedo.
cp (float or callable, optional): The planet's isobaric specific heat capacity in J/kg/K.
mlDepth (float, optional): The depth of the planet's mixed layer (in units of m for rock/water models,
or Pa for gas/bell2018 models).
mlDensity (float, optional): The density of the planet's mixed layer (in kg/m^3) if rock/water models,
or the inverse of the planet's surface gravity (in s^2/m) for gas/bell2018 models.
T_exponent (float): The exponent which determinges the rate at which the body cools (4 for blackbody cooling,
1 for Newtonian cooling).
emissivity (float, optional): The emissivity of the emitting layer (between 0 and 1).
trasmissivity (float, optional): The trasmissivity of the emitting layer (between 0 and 1).
nlat (int, optional): The number of latitudinal cells to use for rectangular maps.
nlon (int, optional): The number of longitudinal cells to use for rectangular maps.
If nlon==None, uses 2*nlat.
useHealpix (bool, optional): Whether the planet's map uses a healpix grid.
nside (int, optional): A parameter that sets the resolution of healpix maps.
"""
#Planet Attributes
self.plType = plType
self.rad = rad*const.R_jup.value # m
self.mass = mass*const.M_jup.value # kg
self.g = const.G.value*self.mass/self.rad**2 # m/s^2
self.albedo = albedo # None
self.cp = cp
self.mlDepth = mlDepth
self.mlDensity = mlDensity
self.T_exponent = T_exponent
if emissivity > 1:
self.emissivity = 1
elif emissivity < 0:
self.emissivity = 0
else:
self.emissivity = emissivity
if trasmissivity > 1:
self.emissivity = 1
elif trasmissivity < 0:
self.trasmissivity = 0
else:
self.trasmissivity = trasmissivity
self.absorptivity = 1-self.albedo-self.trasmissivity
# Planet's Thermal Attributes
if self.plType.lower()=='water':
#water
if self.cp == None:
self.cp = 4.1813e3 # J/kg/K
if self.mlDepth == None:
self.mlDepth = 50 # m
if self.mlDensity == None:
self.mlDensity = 1e3 # kg/m^3
self.C = self.mlDepth*self.mlDensity*self.cp
elif self.plType.lower() == 'rock':
#basalt
if self.cp == None:
self.cp = 0.84e3 # J/kg/K
if self.mlDepth == None:
self.mlDepth = 0.5 # m
if self.mlDensity == None:
self.mlDensity = 3e3 # kg/m^3
self.C = self.mlDepth*self.mlDensity*self.cp
elif self.plType.lower() == 'gas':
# H2 atmo
if self.cp == None:
self.cp = 14.31e3 # J/kg/K
if self.mlDepth == None:
self.mlDepth = 0.1e5 # Pa
if self.mlDensity == None:
self.mlDensity = 1/self.g # s^2/m
self.C = self.mlDepth*self.mlDensity*self.cp
elif self.plType.lower() == 'bell2018':
# LTE H2+H atmo
if self.cp == None:
self.cp = h2.true_cp
if self.mlDepth == None:
self.cpParams = None
if self.mlDensity == None:
self.mlDepth = 0.1e5 # Pa
self.mlDensity = 1/self.g # s^2/m
else:
print('Planet type not accepted!')
return False
#Map Attributes
self.map = Map(nlat=nlat, nlon=nlon, useHealpix=useHealpix, nside=nside)
#Orbital Attributes
self.a = a*const.au.value # m
self.Porb = Porb # days (if None, will be set to Kepler expectation when loaded into system)
self.e = e # None
self.Omega = Omega # degrees ccw from line-of-sight
self.inc = inc # degrees above face-on
self.argp = argp # degrees ccw from Omega
self.obliq = obliq # degrees toward star
if self.obliq <= 90:
self.ProtSign = 1
else:
self.ProtSign = -1
self.argobliq = argobliq # degrees from t0
self.t0 = t0 # days
if self.Porb is not None:
self.orbit = KeplerOrbit(a=self.a, Porb=self.Porb, inc=self.inc, t0=self.t0,
e=self.e, Omega=self.Omega, argp=self.argp, m2=self.mass)
else:
self.orbit = None
#Rotation Rate Attributes
if Prot is None:
self.Prot_input = self.Porb
else:
self.Prot_input = Prot # days
if self.Prot_input is None:
self.vRot = 0
else:
self.vRot = 2*np.pi*self.rad/(self.Prot_input*24*3600) # m/s
self.vWind = vWind # m/s
if self.vWind == 0:
self.Prot = self.Prot_input
else:
self.Prot = 2*np.pi*self.rad/((self.vWind+self.vRot)*(24*3600))
# Used to propogate any changes to the planet's attributes through the other attributes
[docs] def update(self):
"""Update the planet's properties
Used to propogate any manual changes to the planet's attributes through the other, dependent attributes.
"""
g_old = self.g
self.g = const.G.value*self.mass/self.rad**2
if self.plType.lower() == 'gas':
# H2 atmo
if self.mlDensity == 1/g_old:
self.mlDensity = 1/self.g # s^2/m
self.C = self.mlDepth*self.mlDensity*self.cp
elif self.plType.lower() == 'bell2018':
# LTE H2+H atmo
if self.mlDensity == 1/g_old:
self.mlDensity = 1/self.g # s^2/m
else:
self.C = self.mlDepth*self.mlDensity*self.cp
if self.Porb is not None and self.Prot is None:
self.Prot_input = self.Porb
if self.Prot_input is not None:
self.vRot = 2*np.pi*self.rad/(self.Prot_input*24*3600) # m/s
if self.vWind == 0:
self.Prot = self.Prot_input
else:
self.Prot = 2*np.pi*self.rad/((self.vWind+self.vRot)*(24*3600))
# if self.Porb is not None:
# self.orbit = KeplerOrbit(a=self.a, Porb=self.Porb, inc=self.inc, t0=self.t0,
# e=self.e, Omega=self.Omega, argp=self.argp, m2=self.mass)
[docs] def ssp(self, t):
"""Calculate the sub-stellar longitude and latitude.
Args:
t (ndarray): The time in days.
Returns:
list: A list of 2 ndarrays containing the sub-stellar longitude and latitude.
Each ndarray is in the same shape as t.
"""
if type(t)!=np.ndarray:
t = np.array([t])
tshape = t.shape
elif len(t.shape)!=1:
tshape = t.shape
t = t.reshape(-1)
else:
tshape = t.shape
trueAnom = (self.orbit.true_anomaly(t)*180/np.pi)
trueAnom[trueAnom<0] += 360
sspLon = (trueAnom-(t-self.orbit.peri_time())/self.Prot*360)
sspLon = sspLon%180+(-180)*(np.rint(np.floor(sspLon%360/180) > 0))
sspLat = self.obliq*np.cos(t/self.Porb*2*np.pi-self.argobliq*np.pi/180)
return sspLon.reshape(tshape), sspLat.reshape(tshape)
[docs] def sop(self, t):
"""Calculate the sub-observer longitude and latitude.
Args:
t (ndarray): The time in days.
Returns:
list: A list of 2 ndarrays containing the sub-observer longitude and latitude.
Each ndarray is in the same shape as t.
"""
if type(t)!=np.ndarray:
t = np.array([t])
sopLon = 180-((t-self.orbit.t0)/self.Prot)*360
sopLon = sopLon%180+(-180)*(np.rint(np.floor(sopLon%360/180) > 0))
sopLat = 90-self.inc-self.obliq
return sopLon, sopLat
[docs] def Fout(self, T=None, bolo=True, wav=1e-6):
"""Calculate the instantaneous total outgoing flux.
Args:
T (ndarray): The temperature (if None, use self.map.values).
bolo (bool, optional): Determines whether computed flux is bolometric
(True, default) or wavelength dependent (False).
wav (float, optional): The wavelength to use if bolo==False.
Returns:
ndarray: The emitted flux in the same shape as T.
"""
if T is None:
T = self.map.values
elif type(T)!=np.ndarray:
T = np.array([T])
if bolo:
return self.emissivity*const.sigma_sb.value*T**self.T_exponent
else:
a = (2*const.h.value*const.c.value**2/wav**5)
b = (const.h.value*const.c.value)/(wav*const.k_B.value)
return self.emissivity*a/np.expm1(b/T)
[docs] def weight(self, t, refPos='SSP'):
"""Calculate the weighting of map pixels.
Weight flux by visibility/illumination kernel, assuming the star/observer are infinitely far away for now.
Args:
t (ndarray): The time in days.
refPos (str, optional): The reference position; SSP (sub-stellar point) or SOP (sub-observer point).
Returns:
ndarray: The weighting of map mixels at time t. Has shape (t.size, self.map.npix).
"""
if type(t)!=np.ndarray or len(t.shape)==1:
t = np.array([t]).reshape(-1,1)
if refPos.lower() == 'ssp':
refLon, refLat = self.ssp(t)
elif refPos.lower() == 'sop':
refLon, refLat = self.sop(t)
else:
print('Reference point "'+str(refPos)+'" not understood!')
return False
weight = (np.cos(self.map.latGrid*np.pi/180)*np.cos(refLat*np.pi/180)*np.cos((self.map.lonGrid-refLon)*np.pi/180)+
np.sin(self.map.latGrid*np.pi/180)*np.sin(refLat*np.pi/180))
weight = np.max(np.append(np.zeros_like(weight[np.newaxis,:]), weight[np.newaxis,:], axis=0), axis=0)
return weight
[docs] def Fp_vis(self, t, T=None, bolo=True, wav=4.5e-6):
"""Calculate apparent outgoing planetary flux (used for making phasecurves).
Weight flux by visibility/illumination kernel, assuming the star/observer are infinitely far away for now.
Args:
t (ndarray): The time in days.
T (ndarray): The temperature (if None, use self.map.values).
bolo (bool, optional): Determines whether computed flux is bolometric
(True, default) or wavelength dependent (False).
wav (float, optional): The wavelength to use if bolo==False
Returns:
ndarray: The apparent emitted flux. Has shape (t.size, self.map.npix).
"""
if T is None:
T = self.map.values
if type(t)!=np.ndarray or len(t.shape)==1:
t = np.array(t).reshape(-1,1)
weights = self.weight(t, 'SOP')
flux = self.Fout(T, bolo, wav)*self.map.pixArea*self.rad**2
# used to try to remove wiggles from finite number of pixels coming in and out of view
weightsNormed = weights*(4*np.pi/self.map.npix)/np.pi
return np.sum(flux*weights, axis=1)#/np.sum(weightsNormed, axis=1)
[docs] def plot_map(self, tempMap=None, time=None):
"""A convenience routine to plot the planet's temperature map.
Args:
tempMap (ndarray): The temperature map (if None, use self.map.values).
time (float, optional): The time corresponding to the map used to de-rotate the map.
Returns:
figure: The figure containing the plot.
"""
if tempMap is None:
if time is None:
subStellarLon = self.ssp(self.map.time)[0].flatten()
else:
subStellarLon = self.ssp(time)[0].flatten()
else:
self.map.set_values(tempMap, time)
if time is not None:
subStellarLon = self.ssp(time)[0].flatten()
else:
subStellarLon = None
return self.map.plot_map(subStellarLon)
[docs] def plot_dissociation(self, tempMap=None, time=None):
"""A convenience routine to plot the planet's H2 dissociation map.
Args:
tempMap (ndarray, optional): The temperature map (if None, use self.map.values).
time (float, optional): The time corresponding to the map used to de-rotate the map.
Returns:
figure: The figure containing the plot.
"""
if tempMap is None:
if time is None:
subStellarLon = self.ssp(self.map.time)[0].flatten()
else:
subStellarLon = self.ssp(time)[0].flatten()
else:
self.map.set_values(tempMap, time)
if time is not None:
subStellarLon = self.ssp(time)[0].flatten()
else:
subStellarLon = None
self.map.plot_dissociation(subStellarLon)
return plt.gcf()