File:Habitable zone over time sun mass metallicity 002 1 1.svg

2. 1. visualize starpasta output
3. to habitable zone vs time diagram
5. python 3.8 source code
6. https://github.com/hersfeldtn/Starpasta
9. 24.4.2023 0000.0000

import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib from scipy import interpolate import math

def calc_earthline_au(lstar1, teff1): ## original snowline in planet system! earthline_au=1*np.power(lstar1,1/2 )*np.power((teff1/5780), -0.799) #print(" L Teff snowline ", lstar1, teff1, snowline_au) return(earthline_au)

1. 1. in datafile
2. step, stage, t, mt, Mc, McCO, ML, L, R, Rcr, Teff, hzoptin, hzconin, hzconout, hzoptout

1. 1. starpasta file name

1. 1. here Mstar=1 msun, Z=0.02

sfilen1="Z0.02_M1.0.csv"

1. 1. sep=","

data1=pd.read_csv(sfilen1, header=None)

1. data2 = data1.loc[data1[1] == 5 ]

time1=data1[2] luminosity1=data1[7] radius1=data1[8] teff1=data1[10] hzin=data1[12] hzout=data1[13] hzin2=data1[11] hzout2=data1[14]

print("t ", time1[1]) print("L ", luminosity1[1]) print("Teff ",teff1[1]) print("R ",radius1[1]) print("Optimistic HZ in middle out ", round(hzin[1],3), round(((hzin[1]+hzout[1])/2), 3), round(hzout[1],4) )

1. time1=data2[2]
2. luminosity1=data2[7]
3. radois1=data2[8]
4. teff1=data2[10]

time1=data1[2] logluminosity1=np.log10(data1[7]) logradius1=np.log10(data1[8]) logteff1=np.log10(data1[10])

1. print(data1)

1. print(time1)

1. quit(-1)

1. print (len(time1))
2. print(len(luminosity1))

f1 = interpolate.interp1d(time1, luminosity1) f2 = interpolate.interp1d(time1, teff1)

1. 1. first 40 Ma

billion=1e9 million=1e6

timenew3 = np.arange(10*billion, 11*billion, 100*million) lumnew3 = f1(luminosity1) teffnew3 = f2(teff1)

1. plt.xlim([1000,10000])
2. plt.plot(teff1, luminosity1)

font = {'family' : 'DejaVu Sans',

       'weight' : 'bold',
       'size'   : 16}

diagram=2

1. 1. H-R

if(diagram==1): plt.xlim([3.4,3.8]) plt.ylim([-1,4]) plt.gca().invert_xaxis() #print (timenew3) #print (lumnew3) len2=int(len(timenew3)/10) for i in range(0,10,1): x=math.log10(teffnew3[i]) y=math.log10(lumnew3[i]) print(x,y) timema1=timenew3[i] timema2=round(timema1,1) timema3=str(timema2) #plt.plot(x,y,color="blue",marker='x') #plt.text(x,y,timema3, color="green",size=15, alpha=0.7)

plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.suptitle("1 Sun mass star, Z=0.02 ", size=17) plt.title("Evolution track", color="green", alpha=0.7, size=15) plt.xlabel("Log10 temperature Teff K", size=17) plt.ylabel("Log10 luminosity L/L_sun", size=17) plt.grid(color='r', linestyle='-', linewidth=1, alpha=0.2) plt.plot(logteff1 , logluminosity1) plt.show()

1. 1. habitable zone vs time

if(diagram==2): earthline1=calc_earthline_au(luminosity1, teff1) plt.suptitle("1 Sun mass star, Z=0.02 ", size=17) plt.title("Habitable zone", color="green", alpha=0.7, size=15) plt.xlabel("Time million years ", size=17) plt.ylabel("Distance AU" , size=17) plt.xticks(fontsize=15) plt.yticks(fontsize=15) #plt.xlim([2.0,4.5]) #plt.ylim([-1,2.5]) plt.ylim(0, 3) plt.xlim(0, 1.2e4) logtime1=np.log10(time1) loghzin1=np.log10(hzin) loghzout1=np.log10(hzout) plt.scatter(4600, 1, s=80, marker="o") plt.axhline(y = 1, color = 'g', linestyle = ":", alpha=0.2, label = "Earth distance") plt.axhline(y = 1.523, color = 'r', linestyle = ":", alpha=0.2, label = "Mars distance") plt.axhline(y = 0.723, color = '#7f0000', linestyle = ":", alpha=0.2, label = "Venus distance") plt.axhline(y = 0.3871, color = '#7f0000', linestyle = ":", alpha=0.2, label = "Mercury distance") #plt.plot(logtime1 ,loghzin1, color="red") #plt.plot(logtime1 ,loghzout1, color="blue") #plt.plot(time1 ,hzin, color="red") #plt.plot(time1 ,hzout, color="blue") plt.fill_between(time1, hzin2, hzout2, color="green", alpha=0.05) plt.fill_between(time1, hzin, hzout, color="green", alpha=0.2)

#plt.plot(time1 ,hzin2, color="red", linestyle = "--", alpha=0.5) #plt.plot(time1 ,hzout2, color="blue", linestyle = "--",alpha=0.5) plt.plot(time1 ,(hzin+hzout)/2, color="green", linestyle = "--", alpha=0.5) plt.plot(time1 ,earthline1, color="green", linestyle = "-", alpha=0.75) ## distance corresponds Earth today #plt.plot(time1 ,np.sqrt(luminosity1)*1.25, color="#7f7fff", linestyle = ":") ## possible lot of CO2 freeze limit plt.show()