File:Annual temperature changes of planet equator ecc 0p1 tilt 30 mvelp 0 1 1 1 1.png

1. -*- coding: utf-8 -*
4. annual temperature changes in specific latitude of planet
5. 1. simple 1d energy balance model, no temperature diffusion
6. no eccentricity effect to slongitude of sun
7. 1. python 3 source code
8. 1. 12.11.2023 v 0000.0002 beta only

import math import numpy as np import scipy

from scipy.ndimage import uniform_filter1d

import matplotlib.pyplot as plt from matplotlib import cm import PIL

def distans(e, t,period):

   p = 1 - e**2
   theta = 2 * math.pi * t / period
   distans = p / (1 + e * np.cos(theta))
   return (distans)

def energy_balance_temperature(simu_orbits, dt, speklen, Scoeff, orbitperiod_earthyrs, tilt_deg,ecc,velp,rotperiod_h,latitude_deg,albedo,tau,density,depth,cap): S1=1361.5*Scoeff #S2=S1*math.cos(math.radians(latitude)) Kelvin=273.15 #albedo=0.3 ## earth sb=5.670374419e-8 no_absorb_scatter=0.925 earthday=24*3600 yearlen=365.25*orbitperiod_earthyrs yearlen_s=yearlen*24*3600 dt=3600 ## s #simulen=int((yearlen/dt))*simu_orbits simulen=int((simu_orbits*orbitperiod_earthyrs*3600*24*365.25)/dt) #print (simulen) #quit(-1) ## teq of planet Teq1=math.pow(((S1*(1-albedo)) / (2*sb)), 1/4) #print(Teq1) rotperiod_s=3600*rotperiod_h Ts=Teq1 tim=0 tees=[] cosphases=[] temps=[] hoos=[]

for tee in range(0,rotperiod_s*simulen, dt): ## approx lam=((tee%yearlen_s)/yearlen_s)*math.pi*2 hoo=(((tee%rotperiod_s)/rotperiod_s)*math.pi*2)-math.pi #print(hoo)

latitude_rad=math.radians(latitude_deg)

## !!! approx, if ecvc small sundekl_rad=math.asin(math.sin(math.radians(tilt_deg))*math.sin(lam-math.radians(velp))) #sundekl=0

sundekl_deg=math.degrees(sundekl_rad)

#sunsethoo=-math.tan(latitude_rad)*math.tan(sundekl) #altangle=-0.83 ## atmpsphere refraction altangle=0

indaya=1 if((90-abs(latitude_deg))<=tilt_deg): indaya=0 inday=0

if(indaya==1): sunsethoo=math.acos((math.sin(altangle)-math.sin(latitude_rad)*math.sin(sundekl_rad))/(math.cos(latitude_rad)*math.cos(sundekl_rad))) inday=1 if(hoo<(-sunsethoo)): inday=0 if(hoo>(sunsethoo)):inday=0

#zenita=math.radians(latitude_deg-sundekl_deg) zenita_rad=math.acos(math.sin(latitude_rad)*math.sin(sundekl_rad)+math.cos(latitude_rad)*math.cos(sundekl_rad)*math.cos (hoo))

radisun=math.cos(zenita_rad)*inday

#sunhmax_lat=90-(latitude_deg-sundekl_deg) #season=sunhmax_lat

tee_velp=(velp/360)*yearlen_s dis=distans(ecc, tee+tee_velp, yearlen_s) sol=1/(dis*dis) #ret1=1+ecc*math.cos(lam-math.radians(velp)) #ret2=1-(ecc*ecc) #sol=math.pow((ret1/ret1),2)

S2=S1*sol*radisun

#phase=(tim%rotperiod_s)/rotperiod_s #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*S3*no_absorb_scatter solrad=S2*no_absorb_scatter asr=(1-albedo)*solrad olr=tau*sb*math.pow(Ts,4) dT=(asr-olr)/(cap*depth*density) Ts=Ts+dT temps.append(Ts) #cosphases.append(cosphase) tees.append(tim) #hoos.append(hoo) tim=tim+dt

tees1=np.array(tees) temps1=np.array(temps) #hoos1=np.array(hoos) datalen=len(tees) planetdaysimulen=int(rotperiod_s/dt) planetyearsimulen=int(yearlen_s/dt) #tees2=tees1[datalen-planetdaysimulen:] #print(rotperiod, dt) tees2=np.linspace(0, rotperiod_s,int(rotperiod_s/dt)) temps2=temps1[datalen-planetdaysimulen:] tees3=np.linspace(0, yearlen_s,int(yearlen_s/dt)) temps3=temps1[datalen-planetyearsimulen:] siz3=len(temps3) filtsiz=24 temps4 = uniform_filter1d(temps3, size=filtsiz) return(tees2, temps2, tees3, temps3)

simu_orbits=10

speklen=3400*24

orbitperiod_earthyrs=1

dt=3600 ## s rotperiod_h=24 ## rotperiod_h

Scoeff=1.0

1. tilt=23.44

1. tilt=23.44

tilt=30

1. ecc=0.0167

ecc=0.1

1. 1. ecc=0.5
2. ecc=0

velp=0

latitude=0

1. albedo=0.125

albedo=0.25

tau=0.612

1. depth=0.1

depth=1

cap=4.186 density=1000

S1=Scoeff*1361.5 Kelvin=273.15

tees2, temps2, tees3, temps3 =energy_balance_temperature(simu_orbits, dt, speklen, Scoeff, orbitperiod_earthyrs, tilt,ecc,velp,rotperiod_h,latitude,albedo,tau,density,depth,cap)

tees1=tees3/(3600*24) temps1=temps3

1. tees1=tees2
2. temps1=temps2

print(tees1) print(temps1)

mintemp1=round((np.min(temps1)-Kelvin),1) maxtemp1=round((np.max(temps1)-Kelvin),1) deltemp1=round((maxtemp1-mintemp1),1) meantemp1=round((np.mean(temps1)-Kelvin),1)

plt.plot(tees1, temps1-Kelvin, color="red")

1. plt.scatter(tees1,temps1-Kelvin, c=cm.viridis(np.abs(temps1-10)), edgecolor='none')
2. plt.scatter(tees1,temps1, c=cm.hot(np.abs(temps1)), edgecolor='none')
3. plt.scatter(tees1,temps1, edgecolor='none')

1. plt.plot(tees1/(3600), temps1)

plt.suptitle("Annual temperature changes of planet", fontsize=18)

plt.title("S="+ str(S1)+" tilt="+str(tilt)+" ecc="+str(ecc)+" velp="+str(velp)+ " latitude="+str(latitude)+" albedo="+str(albedo)+ " active depth="+str(depth)+" m", fontsize=10 ) plt.ylabel("Temperature degC",fontsize=16, color="#3f0000") plt.xlabel("Day of year", fontsize=16)

1. plt.xlabel("hourangle")

plt.xticks(fontsize=12) plt.yticks(fontsize=12)

plt.grid()

1. plt.axhline(y=0, linestyle=":", color="#7f7fff", label="water freezes")

plt.legend()

plt.show()