File:Planet rotation rate effect to daily temperatures simple energy balance model 1.png

1. -*- coding: utf-8 -*
4. planet rotation effect to daily temperature change
   1. simple energy balance model
5. 1. python 3 source code
6. 1. 04.11.2023 v 0000.0001

import math import numpy as np import scipy

import matplotlib.pyplot as plt import PIL

def rotating_planet_energy_balance_temperature(latitude, S1,Ts,albedo, cap, depth, rotperiod, simulen): #Ts=Teq1 tim=0 tees=[] cosphases=[] temps=[]

for tee in range(0,rotperiod*simulen, dt): phase=(tim%rotperiod)/rotperiod phase2=-math.pi+phase*math.pi*2 cosphase=math.cos(phase2) if(phase2<-math.pi/2):cosphase=0 if(phase2>math.pi/2):cosphase=0 solrad=cosphase*S2*no_absorb_scatter asr=(1-albedo)*solrad olr=tau*sb*math.pow(Ts,4) dT=(asr-olr)/(cap*depth) Ts=Ts+dT temps.append(Ts) cosphases.append(cosphase) tees.append(tim) tim=tim+dt

tees1=np.array(tees) temps1=np.array(temps) datalen=len(tees) planetdaysimulen=int(rotperiod/dt) #tees2=tees1[datalen-planetdaysimulen:] #print(rotperiod, dt) tees2=np.linspace(0, rotperiod,int(rotperiod/dt)) #print("LT",len(tees2)) temps2=temps1[datalen-planetdaysimulen:] return(tees2, temps2)

fii=0

hourangle=0

S1=1361.5

lat=0

S2=S1*math.cos(math.radians(lat))

Kelvin=273.15

albedo=0.3

1. albedo=0.3 ## earth

sb=5.670374419e-8 no_absorb_scatter=0.925

tau=0.612

simulen=36

1. 1. seconds

dt=60*1

1. Ts=310

depth=5 cap=4186 albedo=0.3

1. depth=10
2. albedo=0.3
3. cap=800 700, 1200, 1500?

latitude=0

1. 1. teq of planet

Teq1=math.pow(((S1*(1-albedo)) / (2*sb)), 1/4)

hours1=24

roths0=[] mintemps0=[] maxtemps0=[] deltemps0=[] meantemps0=[]

beginh=12 diffh=12

1. maxh=24*60
2. maxh=1000

maxh=24*30

for hours1 in range(beginh,maxh, diffh): rotperiod=3600*hours1*1 tees2, temps2=rotating_planet_energy_balance_temperature(latitude, S1,Teq1,albedo, cap, depth, rotperiod, simulen) mintemp1=round((np.min(temps2)-Kelvin),1) maxtemp1=round((np.max(temps2)-Kelvin),1) deltemp1=round((maxtemp1-mintemp1),1) meantemp1=round((np.mean(temps2)-Kelvin),1) roths0.append(hours1) mintemps0.append(mintemp1) maxtemps0.append(maxtemp1) deltemps0.append(deltemp1) meantemps0.append(meantemp1)

roths0=np.array(roths0) mintemps0=np.array(mintemps0) maxtemps0=np.array(maxtemps0) deltemps0=np.array(deltemps0) meantemps0=np.array(meantemps0)

plt.suptitle("Planet rotation period vs equator daily temperatures", fontsize=18)

plt.title("S="+ str(S1)+" tilt=0 albedo="+str(albedo)+ " active depth="+str(depth)+" m", fontsize=16 ) plt.ylabel("Temperature degC",fontsize=15) plt.xlabel("Rotation period days",fontsize=15) plt.xticks(fontsize=12) plt.yticks(fontsize=12)

plt.plot(roths0/24, mintemps0, lw=2,color="blue", label="T min") plt.plot(roths0/24, maxtemps0, lw=2,color="red", label="T max") plt.plot(roths0/24, meantemps0, lw=2, color="green", label="T mean")

plt.plot(roths0/24, deltemps0, lw=2,color="brown", linestyle=":", label="T max-min") plt.axhline(y=100, linestyle=":", color="#ff3f3f", label="water boils") plt.axhline(y=0, linestyle=":", color="#7f7fff", label="water freezes")

plt.legend()

plt.show()