https://en.wikipedia.org/wiki/TRAPPIST-1

https://github.com/taylorbell57/Bell_EBM

1. coding: utf-8

1. 1. EBM for trappist-1e
   2. based on Bell's energy balance model

1. 1. 18.10.2023 v 0000.0001b

1. 1. you must pip install . from directory or so on Bell EBM from github
   2. https://github.com/taylorbell57/Bell_EBM
   3. you must have Python 3
   4. copy directory to your computer, go the directoty, then "pip install ."
2. 1. https://bell-ebm.readthedocs.io/en/v1.3/

import numpy as np import scipy.ndimage

import matplotlib.pyplot as plt import matplotlib as mpl

import astropy.constants as const

import Bell_EBM as ebm

1. 1. params of planet

albedo=0.32

1. albedo=0.70

1. 1. Trappist-1e

rad=0.920 mass=0.692 a=0.02925

1. 1. Trappist-1c

1. rad=1.097
2. mass=1.308
3. a=0.0158
4. S=2.214
   1. Teq 339.7 K 66.6 C

1. 1. Trappist-1f

1. rad=1.045
2. mass=1.039
3. a=0.03949
4. S=0.242
5. Teq=217.7, -55.5

e=0.0 argp=0 obliquity=0

star_teff=2550. star_rad=0.114 star_mass=0.089

1. star_teff=2511.
2. star_rad=0.121
3. star_mass=0.08

simulen=10 simustepj=1000

cp_N2 = 1.039e3 # J/(kg K) cp_H2O = 4.2e6 # J/(kg K) cp_rock = .800e6 # J/(kg K)

1. We'll pretend the whole atmosphere absorbs and radiates

P0 = const.atm.value

1. 1. quite flat, maybe graanhouse model!

1. planet = ebm.Planet('water', rad=const.R_earth.value*rad, mass=const.M_earth.value*mass,a=a*const.au.value, e=e, argp=argp, albedo=albedo)

planet = ebm.Planet('rock', rad=const.R_earth.value*rad, mass=const.M_earth.value*mass,a=a*const.au.value, e=e, argp=argp, albedo=albedo)

1. planet = ebm.Planet('gas', rad=const.R_earth.value*rad, mass=const.M_earth.value*mass,a=a*const.au.value, e=e, argp=argp,cp=cp_N2, mlDepth=P0, albedo=albedo)

1. star = ebm.Star()

star = ebm.Star(teff=star_teff, rad=star_rad, mass=star_mass)

system = ebm.System(star, planet)

1. system.planet.Prot = system.planet.Porb*(2./3.)

system.planet.Prot = system.planet.Porb*(1./1.)

Teq = np.median(system.get_teq(np.linspace(0.,system.planet.Porb,simustepj, endpoint=False))) T0 = Teq*np.ones_like(system.planet.map.values) t0 = 0 t1 = t0+system.planet.Porb*simulen dt = system.planet.Porb/simustepj times, maps = system.run_model(T0, t0, t1, dt, verbose=False)

T0 = maps[-1] t0 = times[-1] t1 = t0+system.planet.Porb dt = system.planet.Porb/simustepj times, maps = system.run_model(T0, t0, t1, dt, verbose=False, intermediates=True)

phases = system.get_phase(times) phasePeri = system.get_phase_periastron() indexPeri = np.argmin(np.abs(phases-phasePeri)) subStellLon = system.planet.orbit.get_ssp(system.planet.orbit.t_peri)[0]

1. mapp1=maps[indexPeri]-273.15

1. mapps1=maps[indexPeri]

1. 1. take last orbit to average

mapps1=maps[-simustepj:]

1. print(np.shape(mapps1))

mappa1=np.mean(mapps1, axis=0)

mapp1=mappa1-273.15

tmean1=round(np.mean(mapp1),2) tmin1=round(np.min(mapp1),2) tmax1=round(np.max(mapp1),2)

print(" Temperature degC mean(grid) min max " , tmean1, tmin1, tmax1)

1. quit(-1)

def fmt(x):

   s = f"{x:.2f}"
   return rf"{s} "

plt.title("TRAPPIST-1e EBM simulation: temperature deg C", fontsize=16) plt.xticks(fontsize=14) plt.yticks(fontsize=14)

mapp3 = scipy.ndimage.zoom(mapp1, 4)

levels1=[-300,-200,-180,-150,-120,-110,-100,-90,-80,-70,-60,-50,-40,-20,-10,-5,0,5,10,15,20,25,30,35,40,45,50,60,70,80,90,100,110,120,130,140,150,200,250,300,400,500,600,700,800,1000,2000,3000,4000]

contourcolours1=["#00003f"]

vmin1=-80 vmax1=40

cmap1="coolwarm" cmap1="jet" cmap1="bwr" cmap1="RdBu_r" cmap1="seismic" cmap1="RdYlBu_r" cmap1="Spectral_r" cmap1="hsv_r" cmap1="turbo"

1. cmap1="rainbow"
2. cmap1="gist_rainbow_r"

plt.imshow(mapp3,interpolation= "bicubic", cmap=cmap1, origin="lower", vmin=vmin1, vmax=vmax1, extent=[-180, 180, -90, 90]) CS=plt.contour(mapp3,alpha=0.5, extent=[-180, 180, -90, 90], levels=levels1 ,origin='lower', colors=contourcolours1) plt.clabel(CS, fontsize=14)

plt.show()